-
Francesco Sabatini authoredFrancesco Sabatini authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
A97_ElevationExtract.R 4.75 KiB
ElevationExtract <- function(header, output, ncores, chunk.i){
print("load packages")
require(tidyverse)
require(rgdal)
require(sp)
require(rgeos)
require(raster)
require(rworldmap)
require(elevatr)
require(parallel)
require(doParallel)
print("Import header and divide in tiles")
header.shp <- readOGR(header)
header.tiles <- header.shp@data %>%
bind_cols(as.data.frame(header.shp@coords)) %>%
rename(PlotObservationID=PltObID, Longitude=coords.x1, Latitude=coords.x2) %>%
mutate(lc_ncrt=abs(lc_ncrt)) %>%
filter(lc_ncrt <= 50000) %>%
mutate_at(.vars=vars(Longitude, Latitude),
.funs=list(tile=~cut(., breaks = seq(-180,180, by=.2)))) %>%
mutate(tilenam=factor(paste(Longitude_tile, Latitude_tile)))
print("Get continent")
sPDF <- rworldmap::getMap(resolution="high")
continent.high <- sPDF[,"continent"]
crs(continent.high) <- CRS("+init=epsg:4326")
continent.high@data$continent <- fct_recode(continent.high@data$continent, "South America"="South America and the Caribbean")
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"))
print("go parallel")
cl <- makeForkCluster(ncores, outfile="")
registerDoParallel(ncores)
clusterEvalQ(cl, {
library(rgdal)
library(raster)
library(sp)
library(elevatr)
library(dplyr)
})
print("create list of tiles still to do")
myfiles <- list.files(path = output, pattern = "[A-Za-z]*_[0-9]+\\.RData$")
done <- NULL
done <- as.numeric(unlist(regmatches(myfiles, gregexpr("[[:digit:]]+", myfiles))))
todo <- 1:nlevels(header.tiles$tilenam)
if(length(done)>0) {todo <- todo[-which(todo %in% done)]}
todo <- sample(todo, replace=F)
print(paste(length(todo), "tiles to do"))
print("divide in chunks")
#divide in chunks
todo.chunks <- split(todo, sort(todo%%10))
print(paste("start main foreach loop on chunk n=", chunk.i))
print(paste(length(todo.chunks[[chunk.i]]), "to do"))
foreach(i = todo.chunks[[chunk.i]]) %dopar% {
print(paste("doing", i))
#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 +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 ")) #project to mercator
#retrieve dem raster
tryCatch(raster.tile <- get_elev_raster(sp.tile, z=9, expand=max(sp.tile$lc_ncrt), clip="bbox"),
error = function(e){e}
)
if(!exists("raster.tile")) {
raster.tile <- NA
print(paste("tile", i, "doesn't work!, skip to next"))
save(raster.tile, file = paste(output, "elevation_tile_", i, "failed.RData", sep=""))
rm(raster.tile)
} else {
# 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@data$lc_ncrt,
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)),
#lc_ncrt=sp.tile$lc_ncrt,
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(output, "elevation_tile_", i, ".RData", sep=""))
print(paste(i, "done"))
}}
stopCluster(cl)
}