R/mod08_add_DTM.R

In ecoservR/ecoserv_tool: Natural Capital Mapping Toolkit (UK)

#######################################
### Worflow functions               ###
### for EcoservR                    ###
### Sandra Angers-Blondin           ###
### 26 October 2020                 ###
#######################################

#' Add terrain data
#'
#' This function adds terrain data (elevation) to the basemap. Data input can be one or several rasters of tif or asc format.

#' @param mm The mm object loaded in the environment, can be at various stages of updating.
#' @param studyAreaBuffer The buffered study area generated during mod01 or reloaded when resuming a session.
#' @param projectLog The RDS project log file generated by the wizard app and containing all file paths to data inputs and model parameters
#' @return Saves a project_title_MM_08.RDS file to project folder
#' @export

add_DTM <- function(mm = parent.frame()$mm,
                    studyAreaBuffer = parent.frame()$studyAreaBuffer,
                    projectLog = parent.frame()$projectLog){

   timeA <- Sys.time() # start time


   ## Extract the file paths and other info from project log ----------------------

   output_temp <- projectLog$output_temp
   title <- projectLog$title
   scratch_path <- file.path(output_temp, "ecoservR_scratch")
   if (!dir.exists(scratch_path)) dir.create(scratch_path)

   # Get path
   terrainpath <- projectLog$df[projectLog$df$dataset == "terrain", ][["path"]]  # path to corine data, if available


   if (!is.na(terrainpath) & !is.null(terrainpath)){
      message("Preparing to update baseline with DTM data...")

      # DATA IMPORT ---------------------------------------------------------------------------------

      # Check all files at specified path
      dtm <- list.files(terrainpath, # folder with DTM tiles
                        pattern = paste(c(".asc$", ".tif$"),collapse="|"),  # find all tif or asc files
                        ignore.case = TRUE, # sometimes extensions are all caps
                        all.files=TRUE, full.names=TRUE,
                        recursive = TRUE)

      # Read in the dtm tiles
      dtm <- lapply(dtm, function(x) raster::raster(x))


      ## If no projection set, flash a warning and assign OSGB1936

      if (any(unlist(lapply(dtm, function(x) is.na(raster::crs(x)))))){
      message("Warning! No projection set for your DTM so we will assume it is OSGB 1936 (British National Grid).")

         dtm <- lapply(dtm, function(x){
            suppressWarnings({
               raster::crs(x) <- sp::CRS(SRS_string = "EPSG:27700")
               return(x)
            })
         })

      }

      # if (is.na(raster::crs(dtm[[1]]))){
      #
      #    message("Warning! No projection set for your DTM so we will assume it is OSGB 1936 (British National Grid).")
      #
      #
      #    dtm <- lapply(dtm, function(x){
      #       suppressWarnings({
      #          raster::crs(x) <- sp::CRS(SRS_string = "EPSG:27700")
      #          return(x)
      #       })
      #    })
      #
      # }

      # Calculate side of square (assuming projection is in meters)
      span <- dim(dtm[[1]])[[1]]*raster::res(dtm[[1]])[[1]]

      # TILE IF NEEDED ------------------------------------------------------------------------------

      ## If there is only one raster, we tile using OS grid to match grid ref of mm:

      if (length(dtm) == 1){

         # we can call the object "grid" directly

         gridSA <- suppressWarnings(sf::st_intersection(sf::st_set_crs(ecoservR:::grid, 27700), studyAreaBuffer))  # create gridded study area

         gridSA$TILE_NAME <- droplevels(gridSA$TILE_NAME)  # drop squares that are not in the study area

         gridSA <- gridSA[gridSA$TILE_NAME %in% names(mm),]  # make sure that tiles will be matched to existing mm tiles

         dtm_tiles <- vector("list", length = nrow(gridSA))  # create empty list that we will fill with tiles
         names(dtm_tiles) <- gridSA$TILE_NAME

         # Create the dtm tiles for each 10x10km square in study area, and save to temporary folder defined previously
         for (i in 1:length(dtm_tiles)){

            suppressWarnings({  # make sure extents overlap before forcing a clip, otherwise errors
               if (lengths(
                  sf::st_intersects(gridSA[gridSA$TILE_NAME == names(dtm_tiles)[[i]], ],
                                sf::st_set_crs(sf::st_as_sf(sf::st_as_sfc(sf::st_bbox(dtm[[1]]))), 27700)
                  )) == 0) {next}

            })

            ## wrapping in writeraster rather than using filename= in crop due to bug with dataType FLT4S
            dtm_tiles[[i]] <- raster::writeRaster(
               raster::crop(dtm[[1]], gridSA[gridSA$TILE_NAME == names(dtm_tiles)[[i]],]),
               filename = file.path(scratch_path, paste("dtm_tile", i, sep = "")),
               overwrite = TRUE)

         }

         dtm <- dtm_tiles
         rm(dtm_tiles)

      } else if (span < 10000){

         # If we imported tiles and they are smaller than 10x10 km, we take their center point and place it within a 10km grid ref, and assign the grid ref as name

         rcentro <- lapply(dtm, function(x) raster_centroid(x))  # find centroid (spatial point)
         rcentro <- do.call(rbind, rcentro)  # make into sf df

         # index which grid tiles each raster belongs to (several rasters will belong to one tile)
         rtiles <- sf::st_intersects(rcentro, sf::st_set_crs(ecoservR:::grid, 27700))

         # add the 10km grid ref as names to the rasters
         # This could crash if a rcentro point is exactly across two tiles, but shouldn't be if dtm from any reputable source
         names(dtm) <- ecoservR:::grid[unlist(rtiles),]$TILE_NAME

      } else if (span > 10000){
         # if imported tiles are larger than a 10km tile (unlikely), we just use the old extraction
         # setting names to null ensures that the tile matching workflow fails and reverts to slower extraction

         names(dtm) <- NULL
      }



      # EXTRACTION ----------------------------------------------------------------------------------

      message("Extracting DTM data...")

      ## If we managed to assign names to the dtm tiles, extraction will be one-to-one (merging smaller dtm tiles first if needed)

      mm <- mapply(function(x, n) extractRaster(x, dtm, fun = "mean", tile = n, newcol = "elev"),
                   x = mm,
                   n = names(mm), # passing the names of the tiles will allow to select corresponding raster, making function faster. If user is not working with named tiles, will be read as null and the old function will kick in (slower but works)
                   SIMPLIFY = FALSE)  # absolutely necessary



      # SLOPES ------------------------------------------------------------------

      message("Calculating slopes...")

      # Calculate slopes from the DTM; they automatically keep their name
      slopes <- lapply(dtm, function(x) raster::terrain(x, opt = "slope", unit = "degrees"))

      rm(dtm)

      message("Extracting slopes...")
      # Extract into map
      mm <- mapply(function(x, n) extractRaster(x, slopes, fun = "mean", tile = n, newcol = "slope"),
                   x = mm,
                   n = names(mm), # passing the names of the tiles will allow to select corresponding raster, making function faster. If user is not working with named tiles, will be read as null and the old function will kick in (slower but works)
                   SIMPLIFY = FALSE)  # absolutely necessary

      rm(slopes)

      # SAVE UPDATED MASTER MAP ---------------------------------------------------------------------

      saveRDS(mm, file.path(output_temp, paste0(title,"_MM_08.RDS")))

      # Update the project log with the information that map was updated

      projectLog$last_success <- "MM_08.RDS"

      timeB <- Sys.time() # stop time

      # add performance to log
      projectLog$performance[["add_DTM"]] <- as.numeric(difftime(
         timeB, timeA, units="mins"
      ))

      updateProjectLog(projectLog) # save revised log

      # and delete contents of scratch folder
      cleanUp(scratch_path)


      message(paste0("Finished updating with elevation data. Process took ",
                     round(difftime(timeB, timeA, units = "mins"), digits = 1),
                     " minutes."))



   } else {message("No elevation data input specified.")} # end of running condition


   # Return mm to environment, whether it has been updated or not.
   return({
      invisible({
         mm <<- mm
         projectLog <<- projectLog
      })
   })

} # end of function