diff --git a/code/04_buildHeader.Rmd b/code/04_buildHeader.Rmd
index f481745e13a33db2d6a98f2ba4d399923ffe954e..1f4bb818c9658be16cf081a546ebad8dd66302e5 100644
--- a/code/04_buildHeader.Rmd
+++ b/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")
+
+```
+
+
+
+