Started integrating extrac elevation in buildHeader

code/04_buildHeader.Rmd

@@ -443,6 +443,208 @@ header <- header %>%
 ```
 
 
+## 4.3 Extract elevation
+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)
+```{r create tiles}
+header.tiles <- 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.tiles)` 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](https://cran.r-project.org/web/packages/elevatr/vignettes/introduction_to_elevatr.html#get_raster_elevation_data), which uses the [Terrain Tiles on Amazon Web Services](https://registry.opendata.aws/terrain-tiles/). 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](https://github.com/tilezen/joerd/blob/master/docs/data-sources.md)
+
+```{r }
+library(elevatr)
+
+continent.high.merc <- spTransform(continent.high, CRS( "+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m
++nadgrids=@null +no_defs"))
+
+require(parallel)
+require(doParallel)
+cl <- makeForkCluster(10, outfile="")
+registerDoParallel(cl)
+
+clusterEvalQ(cl, {
+  library(raster)
+  library(sp)
+  library(elevatr)
+  library(dplyr)
+  })
+
+#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.tiles$tilenam)
+todo <- todo[-which(todo %in% done)]
+#foreach(i = 1:nlevels(header.sp$tilenam)) %do% {
+foreach(i = todo) %dopar% {  
+  #create sp and project data
+  if(nrow(header.tiles %>% 
+          filter(tilenam %in% levels(header.tiles$tilenam)[i])) ==0 ) next()
+  sp.tile <- SpatialPointsDataFrame(coords=header.tiles %>% 
+                                    filter(tilenam %in% levels(header.tiles$tilenam)[i]) %>%
+                                    dplyr::select(Longitude, Latitude),
+                                  data=header.tiles %>% 
+                                    filter(tilenam %in% levels(header.tiles$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"))}
+  )
+  # clip dem tile with continent shape
+  raster.tile <- mask(raster.tile, continent.high.merc)
+  
+  #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)`
+                                      buffer=sp.tile@data$`Location uncertainty (m)`, 
+                                      small=T)
+  tmp <- round(mapply( quantile, 
+                            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)`, 
+                       na.rm=T))
+  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)
+}
+stopCluster(cl)
+
+```
+
+
+###### TO CHECK BELOW HERE
+
+For those tiles that failed, extract elevation of remaining plots one by one
+```{r}
+#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
+```{r }
+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
+```{r message=F}
+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  
+```{r scatterplot, cache=T}
+#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")
+
+```
+
+
+
+