-
Francesco Sabatini authoredFrancesco Sabatini authored
title: "sPlot 3.0 - Construction"
subtitle: "Step 2 - Extract elevation"
author: "Francesco Maria Sabatini"
date: "6/21/2019"
output: html_document
Timestamp: r date()
Drafted: Francesco Maria Sabatini
Version: 1.2
Changes to v1.2 - based on dataset sPlot_3.0.1, received on 29/06/2019 from SH. Calculates elevation also for plots without location uncertainty (arbitrarily set to 100 m).
This report describes how elevation data was matched to plot locations.
knitr::opts_chunk$set(echo = TRUE)
library(reshape2)
library(tidyverse)
library(readr)
library(dplyr)
library(data.table)
library(elevatr)
library(sp)
library(raster)
library(knitr)
library(kableExtra)
library(viridis)
library(grid)
library(gridExtra)
library(ggforce)
library(viridis)
library(ggrepel)
Import header data
header <- readr::read_delim("../sPlot_data_export/sPlot 3.0.1_header.csv", locale = locale(encoding = 'UTF-8'),
delim="\t", col_types=cols(
PlotObservationID = col_double(),
PlotID = col_double(),
`TV2 relevé number` = col_double(),
Country = col_factor(),
`Cover abundance scale` = col_factor(),
`Date of recording` = col_date(format="%d-%m-%Y"),
`Relevé area (m²)` = col_double(),
`Altitude (m)` = col_double(),
`Aspect (°)` = col_double(),
`Slope (°)` = col_double(),
`Cover total (%)` = col_double(),
`Cover tree layer (%)` = col_double(),
`Cover shrub layer (%)` = col_double(),
`Cover herb layer (%)` = col_double(),
`Cover moss layer (%)` = col_double(),
`Cover lichen layer (%)` = col_double(),
`Cover algae layer (%)` = col_double(),
`Cover litter layer (%)` = col_double(),
`Cover open water (%)` = col_double(),
`Cover bare rock (%)` = col_double(),
`Height (highest) trees (m)` = col_double(),
`Height lowest trees (m)` = col_double(),
`Height (highest) shrubs (m)` = col_double(),
`Height lowest shrubs (m)` = col_double(),
`Aver. height (high) herbs (cm)` = col_double(),
`Aver. height lowest herbs (cm)` = col_double(),
`Maximum height herbs (cm)` = col_double(),
`Maximum height cryptogams (mm)` = col_double(),
`Mosses identified (y/n)` = col_factor(),
`Lichens identified (y/n)` = col_factor(),
COMMUNITY = col_character(),
SUBSTRATE = col_character(),
Locality = col_character(),
ORIG_NUM = col_double(),
ALLIAN_REV = col_character(),
REV_AUTHOR = col_character(),
Forest = col_logical(),
Grassland = col_logical(),
Wetland = col_logical(),
`Sparse vegetation` = col_logical(),
Shrubland = col_logical(),
`Plants recorded` = col_factor(),
`Herbs identified (y/n)` = col_factor(),
Naturalness = col_factor(),
EUNIS = col_factor(),
Longitude = col_double(),
Latitude = col_double(),
`Location uncertainty (m)` = col_double(),
Dataset = col_factor(),
GUID = col_character()
))
Split data into tiles of 1 x 1 degrees, and create sp::SpatialPointsDataFrame files. Only for plots having a location uncertainty < 50 km. (Include also plots without location uncertainty, arbitrarily set to 100 m)
header.sp <- header %>%
dplyr::select(PlotObservationID, Dataset, Longitude, Latitude, `Location uncertainty (m)`) %>%
mutate(`Location uncertainty (m)`,
list=is.na(`Location uncertainty (m)`),
value=-100) %>%
filter(`Location uncertainty (m)`<= 50000) %>%
mutate_at(.vars=vars(Longitude, Latitude),
.funs=list(tile=~cut(., breaks = seq(-180,180, by=1)))) %>%
filter(!is.na(Longitude_tile) & !is.na(Latitude_tile) ) %>%
mutate(tilenam=factor(paste(Longitude_tile, Latitude_tile))) %>%
mutate(`Location uncertainty (m)`=abs(`Location uncertainty (m)`))
There are r nrow(header.sp)
plots out of r nrow(header)
plots with Location uncertainty <= 50km (or absent).
Extract elevation for each plot. Loops over tiles of 1 x 1°, projects to mercator, and extract elevation for plot coordinates, as well as 2.5, 50, and 97.5 quantiles for a buffer area having a radius equal to the location uncertainty of each plot (but only if location uncertainty < 50 km). DEM derive from package elevatr, which uses the Terrain Tiles on Amazon Web Services. Resolutions of DEM rasters vary by region. I set a zoom factor z=10, which corresponds to ~ 75-150 m. Sources are: SRTM, data.gov.at in Austria, NRCAN in Canada, SRTM, NED/3DEP 1/3 arcsec, data.gov.uk in United Kingdom, INEGI in Mexico, ArcticDEM in latitudes above 60°, LINZ in New Zealand, Kartverket in Norway, as described here
require(parallel)
require(doParallel)
cl <- makeForkCluster(6, outfile="")
registerDoParallel(cl)
clusterEvalQ(cl, {
library(raster)
library(sp)
library(elevatr)
library(dplyr)
})
#Define robustQuantile function
#if the center of the plot is not at elevation <0,
#and not all values are <0, calculate only quantiles for values >0
RobustQuantile <- function(x, probs, center, loc.uncert){
x <- x[!is.na(x)]
if(length(x)==0) return(NA) else {
if(all(x)<0 | center<0) return(stats::quantile(x, probs, na.rm=T)) else {
x2 <- x[which(x>=0)]
return(stats::quantile(x2, probs, na.rm=T))
}
}
}
#create list of tiles for which dem could not be extracted
myfiles <- list.files("../_derived/elevatr/")
done <- as.numeric(unlist(regmatches(myfiles, gregexpr("[[:digit:]]+", myfiles))))
todo <- 1:nlevels(header.sp$tilenam)
todo <- todo[-which(todo %in% done)]
#foreach(i = 1:nlevels(header.sp$tilenam)) %do% {
foreach(i = todo) %do% {
#create sp and project data
if(nrow(header.sp %>% filter(tilenam %in% levels(header.sp$tilenam)[i])) ==0 ) next()
sp.tile <- SpatialPointsDataFrame(coords=header.sp %>%
filter(tilenam %in% levels(header.sp$tilenam)[i]) %>%
dplyr::select(Longitude, Latitude),
data=header.sp %>%
filter(tilenam %in% levels(header.sp$tilenam)[i]) %>%
dplyr::select(-Longitude, -Latitude),
proj4string = CRS("+init=epsg:4326"))
sp.tile <- spTransform(sp.tile, CRSobj = CRS("+init=epsg:3857")) #project to mercator
#retrieve dem raster
tryCatch(raster.tile <- get_elev_raster(sp.tile, z=10, expand=max(sp.tile$`Location uncertainty (m)`), clip="bbox"),
error = function(e){next(paste("tile", i, "doesn't work!, skip to next"))}
)
#extract and summarize elevation data
elev.tile <- raster::extract(raster.tile, sp.tile, small=T)
elev.tile.buffer <- raster::extract(raster.tile, sp.tile, buffer=sp.tile$`Location uncertainty (m)`, small=T)
tmp <- round(mapply( RobustQuantile,
x=elev.tile.buffer,
center=elev.tile,
probs=rep(c(0.025, 0.5, 0.975), each=length(elev.tile)),
loc.uncert=sp.tile$`Location uncertainty (m)`))
elev.q95 <- setNames(data.frame(matrix(tmp, ncol = 3, nrow = length(elev.tile.buffer))),
c("Elevation_q2.5", "Elevation_median", "Elevation_q97.5"))
output.tile <- data.frame(PlotObservationID=sp.tile$PlotObservationID,
elevation=round(elev.tile),
elev.q95,
DEM.res=res(raster.tile)[1])
#save output
save(output.tile, file = paste("../_derived/elevatr/elevation_tile_", i, ".RData", sep=""))
print(i)
}
For those tiles that failed, extract elevation of remaining plots one by one
#create list of tiles for which dem could not be extracted
myfiles <- list.files("../_derived/elevatr/")
done <- as.numeric(unlist(regmatches(myfiles, gregexpr("[[:digit:]]+", myfiles))))
todo <- 1:nlevels(header.sp$tilenam)
todo <- todo[-which(todo %in% done)]
#create SpatialPointsDataFrame
sp.tile0 <- SpatialPointsDataFrame(coords=header.sp %>%
filter(tilenam %in% levels(header.sp$tilenam)[todo]) %>%
dplyr::select(Longitude, Latitude),
data=header.sp %>%
filter(tilenam %in% levels(header.sp$tilenam)[todo]) %>%
dplyr::select(-Longitude, -Latitude),
proj4string = CRS("+init=epsg:4326"))
sp.tile0 <- spTransform(sp.tile0, CRSobj = CRS("+init=epsg:3857")) #project to mercator
output.tile <- data.frame(NULL)
#Loop over all plots
for(i in 1:nrow(sp.tile0)){
sp.tile <- sp.tile0[i,]
tryCatch(raster.tile <- get_elev_raster(sp.tile, z=10,
expand=max(sp.tile$`Location uncertainty (m)`)),
error = function(e){bind_rows(output.tile,
data.frame(PlotObservationID=sp.tile$PlotObservationID,
elevation=NA,
Elevation_q2.5=NA,
Elevation_median=NA,
Elevation_q97.5=NA,
DEM.res=NA))
print(paste("could not retrieve DEM for", sp.tile$PlotObservationID))}
)
#extract and summarize elevation data
elev.tile <- raster::extract(raster.tile, sp.tile, small=T)
elev.tile.buffer <- raster::extract(raster.tile, sp.tile,
buffer=sp.tile$`Location uncertainty (m)`, small=T)
#elev.q95 <- t(round(sapply(elev.tile.buffer,stats::quantile, probs=c(0.025, 0.5, 0.975), na.rm=T)))
elev.q95 <- t(round(mapply( RobustQuantile,
x=elev.tile.buffer,
center=elev.tile,
probs=rep(c(0.025, 0.5, 0.975), each=length(elev.tile)),
loc.uncert=sp.tile$`Location uncertainty (m)`)))
output.tile <- bind_rows(output.tile,
data.frame(PlotObservationID=sp.tile$PlotObservationID,
elevation=round(elev.tile),
elev.q95,
DEM.res=res(raster.tile)[1]) %>%
rename(Elevation_q2.5=X2.5., Elevation_median=X50., Elevation_q97.5=X97.5.))
}
save(output.tile, file = paste("../_derived/elevatr/elevation_tile_", 0, ".RData", sep=""))
Compose tiles into a single output, and export
myfiles <- list.files(path)
myfiles <- myfiles[grep(pattern="*.RData$", myfiles)]
#create empty data.frame
elevation.out <- data.frame(PlotObservationID=header$PlotObservationID,
elevation=NA,
Elevation_q2.5=NA,
Elevation_median=NA,
Elevation_q97.5=NA,
DEM.res=NA)
for(i in 1:length(myfiles)){
load(paste(path, myfiles[i], sep=""))
#attach results to empty data.frame
mymatch <- match(output.tile$PlotObservationID, elevation.out$PlotObservationID)
elevation.out[mymatch,] <- output.tile
if(i %in% seq(1,length(myfiles), by=25)){print(i)}
}
write_csv(elevation.out, path ="../_derived/elevatr/Elevation_out_v301.csv")
Reimport output and check
elevation.out <- read_csv("../_derived/elevatr/Elevation_out.csv")
knitr::kable(head(elevation.out,10),
caption="Example of elevation output") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = F, position = "center")
summary(elevation.out)
There are r sum(elevation.out$elevation < -100, na.rm=T)
plots with elevation < -100 !
Create Scatterplot between measured elevation in the field, and elevation derived from DEM
#join measured and derived elevation
mydata <- header %>%
dplyr::select(PlotObservationID, `Altitude (m)`) %>%
filter(!is.na(`Altitude (m)`)) %>%
rename(elevation_measured=`Altitude (m)`) %>%
left_join(elevation.out %>%
dplyr::select(PlotObservationID, elevation) %>%
rename(elevation_dem=elevation),
by="PlotObservationID")
ggplot(data=mydata) +
geom_point(aes(x=elevation_measured, y=elevation_dem), alpha=1/10, cex=0.8) +
theme_bw() +
geom_abline(slope=0, intercept=0, col=2, lty=2) +
geom_abline(slope=1, intercept=1, col="Dark green")
Check strange values (elevation < -100 m.s.l.)
strange <- header %>%
filter(PlotObservationID %in% (elevation.out %>%
filter(elevation< -100))$PlotObservationID) %>%
left_join(elevation.out %>%
dplyr::select(PlotObservationID, elevation),
by="PlotObservationID") %>%
rename(elevation_DEM=elevation)
knitr::kable(strange %>%
dplyr::select(PlotObservationID, `TV2 relevé number`,
"Altitude (m)", elevation_DEM, Longitude:Dataset),
caption="Example of elevation output") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = F, position = "center")
Create graph for plots with elevation < -100
countries <- map_data("world")
robust.range <- function(x){
return(c(floor((min(x, na.rm=T)-0.01)/5)*5,
ceiling((max(x, na.rm=T)+0.01)/5)*5))
}
strange <- strange %>%
mutate_at(.vars=vars(Longitude, Latitude),
.funs=list(tile=~cut(., breaks = seq(-180,180, by=10)))) %>%
mutate(tilenam=factor(paste(Longitude_tile, Latitude_tile)))
for(i in 1:nlevels(strange$tilenam)){
strange.i <- strange %>%
filter(tilenam==levels(strange$tilenam)[i])
sp.i <- SpatialPoints(coords=strange.i %>%
dplyr::select(Longitude, Latitude),
proj4string = CRS("+init=epsg:4326"))
dem.i <- get_elev_raster(sp.i, z=ifelse(i!=15,5,8), expand=ifelse(i!=15,5,0.1)) #API doesn't work for tile i
dem.df <- data.frame(xyFromCell(dem.i, cell=1:ncell(dem.i)), z=getValues(dem.i))
ggstrange <- ggplot(countries, aes(x=long, y=lat, group = group)) +
geom_raster(data=dem.df, aes(x=x, y=y, fill=z, group=1)) +
geom_polygon(col=gray(0.3), lwd=0.3, fill = gray(0.9), alpha=1/5) +
geom_point(data=strange.i,
aes(x=Longitude, y=Latitude, group=1), col="red")+
coord_equal(xlim= robust.range(strange.i$Longitude),
ylim= robust.range(strange.i$Latitude)) +
geom_text_repel(data=strange.i,
aes(x=Longitude, y=Latitude, group=1,
label=strange.i$PlotObservationID)
) +
ggtitle(strange.i$Dataset[1], subtitle = strange.i$Country[1]) +
scale_fill_viridis() +
theme_bw() +
theme(axis.title = element_blank())
print(ggstrange)
}
sessionInfo()