From cc8378b2de37c0ec8bc0f8deeba5c2c01f019319 Mon Sep 17 00:00:00 2001
From: Francesco Sabatini <francesco.sabatini@idiv.de>
Date: Fri, 21 Jun 2019 18:44:31 +0200
Subject: [PATCH] Expanded 01_Extract_elevation to run in tiles. Added chunk to
recompose output
---
.gitignore | 3 +-
code/01_Extract_elevation.Rmd | 78 ++++++++++++++++++++++++++++++-----
2 files changed, 70 insertions(+), 11 deletions(-)
diff --git a/.gitignore b/.gitignore
index 07bb462..34e0cc9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
*.html
code/*_cache/*
code/*_files/*
-sPlot_data_export/*
\ No newline at end of file
+sPlot_data_export/*
+_derived/*
\ No newline at end of file
diff --git a/code/01_Extract_elevation.Rmd b/code/01_Extract_elevation.Rmd
index f9f0d51..fc37717 100644
--- a/code/01_Extract_elevation.Rmd
+++ b/code/01_Extract_elevation.Rmd
@@ -91,25 +91,83 @@ header <- readr::read_csv("../sPlot_data_export/sPlot_data_header_fix1.csv",
))
```
-Split data into tiles of 5 x 5 degrees, and create sp files
+Split data into tiles of 5 x 5 degrees, and create sp files. Only for plots having a location uncertainty < 50 km
```{r create tiles}
header.sp <- header %>%
- dplyr::select(PlotObservationID, Longitude, Latitude) %>%
+ dplyr::select(PlotObservationID, Dataset, Longitude, Latitude, `Location uncertainty (m)`) %>%
mutate_at(.vars=vars(Longitude, Latitude),
- .funs=list(tile=~cut(., breaks = seq(-180,180, by=5)))) %>%
+ .funs=list(tile=~cut(., breaks = seq(-180,180, by=3)))) %>%
filter(!is.na(Longitude_tile) & !is.na(Latitude_tile) ) %>%
- mutate(tilenam=factor(paste(Longitude_tile, Latitude_tile)))
+ mutate(tilenam=factor(paste(Longitude_tile, Latitude_tile))) %>%
+ filter(`Location uncertainty (m)`<= 50000) %>%
+ mutate(`Location uncertainty (m)`=abs(`Location uncertainty (m)`))
+
#group_by(Longitude_tile, Latitude_tile) %>%
#summarize(n())
+```
+
+Extract elevation for each plot. Loops over tiles of 5 x 5°, 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 eval=F}
+require(parallel)
+require(doParallel)
+cl <- makeForkCluster(8, outfile="")
+registerDoParallel(cl)
-#create sp
-sp.tile <- SpatialPoints(coords=header.sp %>%
- filter(tilenam %in% levels(header.sp$tilenam)[1]) %>%
- dplyr::select(Longitude, Latitude),
- proj4string = CRS("+init=epsg:4326"))
+clusterEvalQ(cl, {
+ library(raster)
+ library(sp)
+ library(elevatr)
+ library(dplyr)
+ })
+foreach(i = 1:nlevels(header.sp$tilenam)) %dopar% {
+ print("create sp and project data")
+ 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")
+ raster.tile <- get_elev_raster(sp.tile, z=10, expand=max(sp.tile$`Location uncertainty (m)`), clip="bbox")
+ 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)
+}
+
+```
+Compose output
+```{r}
+path <- "../_derived/elevatr/"
+#create empty data.frame
+elevation.out <- data.frame(PlotObservationID=header.sp$PlotObservationID,
+ elevation=NA,
+ Elevation_q2.5=NA,
+ Elevation_median=NA,
+ Elevation_q97.5=NA,
+ DEM.res=NA)
+for(i in 1:nlevels(header.sp$tilenam)){
+ load(paste(path, list.files(path)[i], sep=""))
+ #attach results to empty data.frame
+ mymatch <- match(output.tile$PlotObservationID, elevation.out$PlotObservationID)
+ elevation.out[mymatch,] <- output.tile
+}
-raster.tile <- get_elev_point(sp.tile )
+write_csv(elevation.out, file="../_derived/elevatr/Elevation_out.csv")
```
--
GitLab