Skip to content
Snippets Groups Projects
Commit cc8378b2 authored by Francesco Sabatini's avatar Francesco Sabatini
Browse files

Expanded 01_Extract_elevation to run in tiles. Added chunk to recompose output

parent 4dbc7d92
Branches
No related tags found
No related merge requests found
*.html *.html
code/*_cache/* code/*_cache/*
code/*_files/* code/*_files/*
sPlot_data_export/* sPlot_data_export/*
\ No newline at end of file _derived/*
\ No newline at end of file
...@@ -91,25 +91,83 @@ header <- readr::read_csv("../sPlot_data_export/sPlot_data_header_fix1.csv", ...@@ -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} ```{r create tiles}
header.sp <- header %>% header.sp <- header %>%
dplyr::select(PlotObservationID, Longitude, Latitude) %>% dplyr::select(PlotObservationID, Dataset, Longitude, Latitude, `Location uncertainty (m)`) %>%
mutate_at(.vars=vars(Longitude, Latitude), 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) ) %>% 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) %>% #group_by(Longitude_tile, Latitude_tile) %>%
#summarize(n()) #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 clusterEvalQ(cl, {
sp.tile <- SpatialPoints(coords=header.sp %>% library(raster)
filter(tilenam %in% levels(header.sp$tilenam)[1]) %>% library(sp)
dplyr::select(Longitude, Latitude), library(elevatr)
proj4string = CRS("+init=epsg:4326")) 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")
``` ```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment