A97_ElevationExtract.R

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 #should take around 40'
  indices <- 1:length(todo)
  todo.chunks <- split(indices, sort(indices%%99))
  
  print(paste("start main foreach loop on chunk n=", chunk.i))
  foreach(i = todo.cunks[[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
      save(raster.tile, file = paste(output, "elevation_tile_", i, "failed.RData", sep=""))
      message(paste("tile", i, "doesn't work!, skip to next"))
      rm(raster.tile)
      next
    }
    # 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)
  }