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

Added scatterplot to 01 elevation, comparing measured and estimated elevation

parent 48314be8
Branches
No related tags found
No related merge requests found
......@@ -124,7 +124,7 @@ clusterEvalQ(cl, {
})
foreach(i = todo[length(todo):1]) %do% {
print("create sp and project data")
#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]) %>%
......@@ -135,11 +135,11 @@ foreach(i = todo[length(todo):1]) %do% {
proj4string = CRS("+init=epsg:4326"))
sp.tile <- spTransform(sp.tile, CRSobj = CRS("+init=epsg:3857")) #project to mercator
print("retrieve dem raster")
#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"))}
)
print("extract and summarize elevation data")
#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)))
......@@ -149,7 +149,7 @@ foreach(i = todo[length(todo):1]) %do% {
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
save(output.tile, file = paste("../_derived/elevatr/elevation_tile_", i, ".RData", sep=""))
print(i)
}
......@@ -158,11 +158,13 @@ foreach(i = todo[length(todo):1]) %do% {
For those tiles that failed, extract elevation of remaining plots one by one
```{r, eval=F}
#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),
......@@ -172,6 +174,7 @@ sp.tile0 <- SpatialPointsDataFrame(coords=header.sp %>%
proj4string = CRS("+init=epsg:4326"))
sp.tile0 <- spTransform(sp.tile, 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,
......@@ -185,7 +188,7 @@ for(i in 1:nrow(sp.tile0)){
DEM.res=NA))
print(paste("could not retrieve DEM for", sp.tile$PlotObservationID))}
)
print("extract and summarize elevation data")
#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)
......@@ -198,7 +201,6 @@ for(i in 1:nrow(sp.tile0)){
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=""))
```
......@@ -224,8 +226,8 @@ for(i in 1:length(myfiles)){
write_csv(elevation.out, path ="../_derived/elevatr/Elevation_out.csv")
```
Check output
```{r}
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") %>%
......@@ -233,10 +235,29 @@ knitr::kable(head(elevation.out,10),
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")
```
Check strange values (elevation < -100 m.s.l.)
......@@ -275,7 +296,7 @@ for(i in 1:nlevels(strange$tilenam)){
sp.i <- SpatialPoints(coords=strange.i %>%
dplyr::select(Longitude, Latitude),
proj4string = CRS("+init=epsg:4326"))
dem.i <- get_elev_raster(sp.i, z=10)
dem.i <- get_elev_raster(sp.i, z=10, expand = 1)
dem.df <- data.frame(xyFromCell(dem.i, cell=1:ncell(dem.i)), z=getValues(dem.i))
ggstrange <- ggplot(countries, aes(x=long, y=lat, group = group)) +
......@@ -291,7 +312,7 @@ for(i in 1:nlevels(strange$tilenam)){
box.padding = 0.35,
point.padding = 0.5,
segment.color = 'grey50') +
ggtitle(strange.i$Dataset[[i]]) +
scale_color_viridis() +
theme_bw() +
theme(axis.title = element_blank())
......@@ -300,6 +321,7 @@ for(i in 1:nlevels(strange$tilenam)){
```
```{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