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

Updated 01 to extract elevation from all plots, and check those with elev < -100

parent cc8378b2
Branches
No related tags found
No related merge requests found
---
title: "sPlot 3.0 - Construction"
subtitle: "Extract elevation"
subtitle: "Step 2 - Extract elevation"
author: "Francesco Maria Sabatini"
date: "6/21/2019"
output: html_document
......@@ -41,6 +41,7 @@ library(ggforce)
```
Import header data
```{r import}
header <- readr::read_csv("../sPlot_data_export/sPlot_data_header_fix1.csv",
......@@ -91,22 +92,22 @@ 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. Only for plots having a location uncertainty < 50 km
Split data into tiles of 1 x 1 degrees, and create sp::SpatialPointsDataFrame files. Only for plots having a location uncertainty < 50 km
```{r create tiles}
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=3)))) %>%
.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))) %>%
filter(`Location uncertainty (m)`<= 50000) %>%
mutate(`Location uncertainty (m)`=abs(`Location uncertainty (m)`))
#group_by(Longitude_tile, Latitude_tile) %>%
#summarize(n())
```
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 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)
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 eval=F}
require(parallel)
......@@ -121,8 +122,9 @@ clusterEvalQ(cl, {
library(dplyr)
})
foreach(i = 1:nlevels(header.sp$tilenam)) %dopar% {
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),
......@@ -133,7 +135,9 @@ foreach(i = 1:nlevels(header.sp$tilenam)) %dopar% {
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")
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)
......@@ -150,25 +154,140 @@ foreach(i = 1:nlevels(header.sp$tilenam)) %dopar% {
}
```
Compose output
```{r}
path <- "../_derived/elevatr/"
For those tiles that failed, extract elevation of remaining plots one by one
```{r, eval=F}
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
```{r, eval=F}
myfiles <- list.files(path)
myfiles <- myfiles[grep(pattern="*.RData$", myfiles)]
#create empty data.frame
elevation.out <- data.frame(PlotObservationID=header.sp$PlotObservationID,
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:nlevels(header.sp$tilenam)){
load(paste(path, list.files(path)[i], sep=""))
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, file="../_derived/elevatr/Elevation_out.csv")
write_csv(elevation.out, path ="../_derived/elevatr/Elevation_out.csv")
```
Check output
```{r}
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.)
```{r}
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
```{r, message=F, warning=F}
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)
}
```
```{r}
sessionInfo()
```
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment