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.0

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)

Import header data

header <- readr::read_csv("../sPlot_data_export/sPlot_data_header_fix1.csv",
                            col_types=cols(
  PlotObservationID = col_double(),
  PlotID = col_double(),
  `TV2 relevé number` = col_double(),
  Country = col_factor(),
  `Cover abundance scale` = col_factor(),
  Author = col_character(),
  `Date of recording` = col_date(format="%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_logical(),
  `Lichens identified (y/n)` = col_logical(),
  Locality = col_character(),
  Naturalness = col_factor(),
  Forest = col_logical(),
  Shrubland = col_logical(),
  Wetland = col_logical(),
  Grassland = col_logical(),
  `Sparse vegetation` = col_logical(),
  EUNIS = col_factor(),
  Longitude = col_double(),
  Latitude = col_double(),
  `Location uncertainty (m)` = col_double(),
  Dataset = col_factor()
))

Split data into tiles of 1 x 1 degrees, and create sp::SpatialPointsDataFrame files. Only for plots having a location uncertainty < 50 km

header.sp <- header %>%
  dplyr::select(PlotObservationID, Dataset, Longitude, Latitude, `Location uncertainty (m)`) %>%
  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.

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(8, outfile="")
registerDoParallel(cl)

clusterEvalQ(cl, {
  library(raster)
  library(sp)
  library(elevatr)
  library(dplyr)
  })

foreach(i = todo[length(todo):1]) %do% {
  print("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

  print("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"))}
  )
  print("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)))
  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.)
  
  print("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

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

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.tile, CRSobj = CRS("+init=epsg:3857")) #project to mercator
output.tile <- data.frame(NULL)
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))}
          )
  print("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)))
  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.csv")

Check output

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) plots with elevation < 100 !

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])
  ggstrange <- ggplot(countries, aes(x=long, y=lat, group = group)) +
    geom_polygon(col=gray(0.3), lwd=0.3, fill = gray(0.9)) + 
    geom_point(data=strange.i, 
              aes(x=Longitude, y=Latitude, group=1), col=2)+ 
    coord_equal(xlim= robust.range(strange.i$Longitude),
                    ylim= robust.range(strange.i$Latitude)) + 
    geom_label(data=strange.i, 
              aes(x=Longitude, y=Latitude, group=1,
                  label=strange.i$PlotObservationID), 
              nudge_y = 0.5)
    theme_bw()  + 
    theme(axis.title = element_blank())
  print(ggstrange)
}

sessionInfo()