01_Extract_elevation.Rmd

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()