R/capacity_flood_reg.R

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

###########################################
### Water purification capacity model   ###
### for EcoservR tool                   ###
### Sandra Angers-Blondin               ###
### 09 Dec 2020                         ###
###########################################

#' Flood Risk Mitigation Capacity Model
#'
#' Runs the flood risk mitigation ecosystem service model, generating capacity scores based on the ability of the landscape to slow down the flow of water. Indicators for this model are vegetation roughness and terrain slope (and optionally soil permeability: surface percent runoff data can be obtained from the Hydrology of Soil Types (HOST) dataset.)

#' @param x A basemap, in a list of sf tiles or as one sf object. Must have attribute HabCode_B.
#' @param studyArea The boundaries of the site, as one sf object. The final raster will be masked to this shape. For best results this shape should be smaller than the basemap (which should be buffered by typically 300 m - 1km to avoid edge effects).
#' @param DTM A raster of elevation values for the study area. If left NULL (default), will be accessed from the projectLog data. Alternatively a file path can be specified.
#' @param spr Basemap attribute name (must be quoted) representing Standard Percentage Runoff values (optional). Default is to leave it out, but if specified it will be included in the scoring.
#' @param res Desired resolution of the raster. Default is 5 m. Range recommended is 5-10m.
#' @param threshold Size (m2) below which an isolated patch is not considered able to provide the service. Default is 500 m2.
#' @param use_hedges Use a separate hedgerow layer? Default FALSE, see clean_hedgerows() for producing a model-ready hedge layer.
#' @param projectLog The RDS project log file generated by the wizard app and containing all file paths to data inputs and model parameters
#' @param runtitle A customised title you can give a specific model run, which will be appended to your project title in the outputs. If comparing a basemap to an intervention map, we recommend using "pre" and "post", or a short description of the interventions, e.g. "baseline" vs "tree planting".
#' @param save Path to folder where outputs will be saved. By default a folder will be created using your chosen run title, prefixed by "services_". Do not use this argument unless you need to save the outputs somewhere else.
#' @return Two rasters with capacity scores: one with raw scores (arbitrary units), and one rescaled 0-100 (where 100 is maximum capacity for the area).
#' @export
#'
capacity_flood_reg <- function(x = parent.frame()$mm,
                                 studyArea = parent.frame()$studyArea,
                                 DTM = NULL, spr = NULL,
                                 res = 5,
                                 threshold = 500,
                                 use_hedges = FALSE,
                                 projectLog = parent.frame()$projectLog,
                                 runtitle = parent.frame()$runtitle,
                                 save = NULL
){

   timeA <- Sys.time() # start time

   # Create output directory automatically if doesn't already exist
   if (is.null(save)){

      save <- file.path(projectLog$projpath,
                        paste0("services_", runtitle))

      if (!dir.exists(save)){
         dir.create(save)
      }
   } else {
      # if user specified their own save directory we check that it's ok
      if(!dir.exists(save) | file.access(save, 2) != 0){
         stop("Save directory doesn't exist, or you don't have permission to write to it.")}
   }

   # Create a temp directory for scratch files

   scratch <- file.path(projectLog$projpath,
                        "ecoservR_scratch")

   if(!dir.exists(scratch)){
      dir.create(scratch)
   }

   # if mm is stored in list, combine all before proceeding
   if (isTRUE(class(x) == "list")){
      message("Recombining basemap tiles")
      x <- do.call(rbind, x) %>% sf::st_as_sf()
   }

   x <- checkcrs(x, 27700)
   studyArea <- checkcrs(studyArea, 27700)

   ### Import DTM -----------

   # Conditional import based on whether the user specified a folder or if we are reading from the project log.

   if (is.null(DTM)){
      # if no dtm specified by user, we look for dataset in project log

      dtm <- try(projectLog$df[projectLog$df$dataset == "terrain",][["path"]])

      if (inherits(dtm, "try-error") | is.na(dtm)){ stop("No elevation dataset found in project log, and you did not specify argument DTM")}



   } else {
      dtm <- DTM

      if(!dir.exists(dtm)){stop("No such folder: ", dtm)}
   }

   ## Now no matter what, the object dtm should be a valid directory
   ## Does it contain a raster?

   if (length(list.files(path = dtm, pattern = "tif", recursive = TRUE)) == 0 &&  # any tif?
      length(list.files(path = dtm, pattern = "asc", recursive = TRUE)) == 0      # or any asc?
   ){
      stop("Cannot find elevation rasters at file path: ", dtm)
   }


   # Now we know there is at least one raster in the folder; list them

   dtm <- list.files(dtm, # folder with DTM tiles
                     pattern = paste0(c('.asc$', '.tif$'), collapse ="|"),
                     all.files=TRUE, full.names=TRUE,
                     recursive = TRUE)


   # Read in the dtm tiles
   message("Importing elevation data...")
   dtm <- lapply(dtm, function(x) raster::raster(x))

   # Check whether tiles have a CRS
   if (is.na(raster::crs(dtm[[1]]))){
      message("Warning! No coordinate reference system associated with DTM data. Assuming OSGB 1936 (British National Grid; EPSG 27700). Check with your data provider if results are unexpected.")

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

   }


   # Keep only those intersecting with basemap extent
   ex <- methods::as(raster::extent(x), "SpatialPolygons")

   dtm <- dtm[sapply(dtm, function(x) rgeos::gIntersects(methods::as(raster::extent(x), "SpatialPolygons"), ex))]

   if (length(dtm) == 0) {stop("DTM doesn't intersect your study area.")} else
      if(length(dtm) == 1){
         dtm <- dtm[[1]] # extract from list
      } else {

   # If tiles, they will be merged (takes a while) and merged output saved for future use.
      message("Merging DTM tiles. This may take a while.")

      dtmpath <- file.path(projectLog$projpath, "DTM")  # decide where we'll save merge output

      if (!dir.exists(dtmpath)){ dir.create(dtmpath) }

      # set up the filename as a property
      dtm$filename <- file.path(dtmpath, paste(projectLog$title, "_merged_DTM.tif", sep = ""))
      dtm$overwrite <- TRUE
      dtm <- do.call(raster::merge, dtm)
      message("Merged DTM saved to: ", dtmpath, " .You can use it to speed up future runs of the model.")
   }


   ## Now we should be working with one single dtm for the rest of the model. Phew!




   ### Check and import hedgerows ----
   if (use_hedges){

      if (!file.exists(projectLog$clean_hedges)){stop("use_hedges is TRUE but no file found. Check projectLog$clean_hedges")}

      hedges <- readRDS(projectLog$clean_hedges) %>%
         dplyr::mutate(HabCode_B = 'J21') %>% dplyr::select(HabCode_B) %>%
         merge(hab_lookup[if (is.null(spr)) c("Ph1code",  "HabBroad", "Rough") else c("HabCode_B", "HabBroad", "Rough", spr)],
               by.x = 'HabCode_B', by.y = 'Ph1code', all.x = TRUE)

      message("Loaded hedges from ", projectLog$clean_hedges)
      hedges <- rename_geometry(hedges, attr(x, "sf_column"))
      hedges <- checkcrs(hedges, 27700)
   }

   # delete HabBroad because gets merged again with lookup
   x$HabBroad <- NULL


   ## Check SPR
   if (!is.null(spr) && !spr %in% names(x)) stop("SPR attribute ",spr, " not found in your basemap.")


   ### Merge the lookup table -----
   message("Adding habitat-specific roughness scores")
   x <- merge(x, hab_lookup, by.x = "HabCode_B", by.y = "Ph1code", all.x = TRUE)

   if (is.null(spr)){ #no HOST column

      x <- x[c("HabCode_B", "HabBroad", "Rough")]  # keep only columns we need

   } else {

      x <- x[c("HabCode_B", "HabBroad", "Rough", spr)]

   }

   x <- checkgeometry(x, "POLYGON") # make sure all is polygon before rasterizing

   ### Calculate slopes -----
   message("Calculating slopes indicator")

   slopes <- raster::terrain(dtm, "slope", "degrees",
                     filename = file.path(scratch, "slopes"), overwrite = TRUE)

   rm(dtm)

   ### Reclassify slopes -----
   # Based on Ecoserv-GIS values attributing score based on steepness
   slopes <- raster::reclassify(slopes,
                        c(0, 5, 100,
                          5, 8, 85,
                          8, 15, 70,
                          15, 25, 60,
                          25, 35, 30,
                          35, 45, 20,
                          45, Inf, 10),
                        include.lowest = TRUE,
                        filename = file.path(scratch, "slopesrcl"), overwrite = TRUE)  # note: saving under new name, overwriting may cause unexpected behaviour


   ### Create raster template with same properties as mastermap -----
   r <-  raster::raster()  # create empty raster
   raster::crs(r) <- sp::CRS(SRS_string = "EPSG:27700") # hard-coding datum to preserve CRS
   raster::extent(r) <-  raster::extent(x)  # set same extent as the shapefile
   raster::res(r) <- res  # set resolution


   ### Rasterize -----

   message("Calculating roughness indicator")

   # Raster of roughness
   rough_r <- raster::writeRaster(
      fasterize::fasterize(x, r, field = "Rough", background = NA),
      filename = file.path(scratch, "waterpurif_rough"),
      overwrite = TRUE  # because it's a temporary file we don't mind overwriting it if it exists
   )

   if(use_hedges){
      hedge_rough_r <- raster::writeRaster(
         fasterize::fasterize(hedges, r, field = "Rough", background = NA),
         filename = file.path(scratch, "waterpurif_hedge_rough"), overwrite = TRUE)

      rough_r <- raster::mask(rough_r, hedge_rough_r, inverse = T) %>% raster::cover(hedge_rough_r)
   }


   ## Align slopes raster (from DTM) with raster template from mastermap
   message("Aligning indicators")

   slopes <- raster::resample(slopes, r, "bilinear",
                      filename = file.path(scratch, "slopes"), overwrite = TRUE)


   ### Multiply rasters ----

   # Score is obtained by combining roughness, slope (and maybe someday condition)
   water_score <- raster::writeRaster(
      slopes * rough_r,
      filename = file.path(scratch, "waterpurif_score"),
      overwrite = TRUE  # because it's a temporary file we don't mind overwriting it if it exists
   )

   rm(slopes, rough_r) # free up memory


   ### Rasterize SPR if present and multiply with previous indicators:

   if (!is.null(spr)){

      message("Calculating soil permeability indicator")

      # Soil run off raster
      soil_r <- raster::writeRaster(
         fasterize::fasterize(x, r, field = spr, background = NA),
         filename = file.path(scratch, "waterpurif_soil"),
         overwrite = TRUE  # because it's a temporary file we don't mind overwriting it if it exists
      )

      # flip the values (high runoff means low permeability)

      if (mean(x[[spr]], na.rm = TRUE) > 1){ct <- 100} else {ct = 1} # guess if data range is 0-1 or 0-100

      soil_r <- 100 - soil_r

      # multiply with water score

      water_score <- raster::writeRaster(
         water_score * soil_r,
         filename = file.path(scratch, "waterpurif_score2"),
         overwrite = TRUE  # because it's a temporary file we don't mind overwriting it if it exists
      )

      rm(ct, soil_r)

   }


   # Raster mask - land only ---
   message("Creating mask for bodies of water")

   watermask <- dplyr::filter(x, HabBroad %in% c("Water, fresh", "Water, sea"))

   if (nrow(watermask) > 0){
   watermask <- fasterize::fasterize(watermask, r, background = NA)
   water_score <- raster::mask(water_score, watermask, inverse = TRUE)  # this sets water areas to NA
   }
   rm(watermask)

   ### Create functional threshold mask ----

   # The function and notes on how it works can be found in the 00_new_functions script

   # This function will create a mask raster based on patch statistics, removing all patches smaller than the threshold.

   message("Applying functional threshold")
   water_score <- functionalMask(water_score,  # the score raster we apply the function to
                                 res = res,  # resolution of the raster
                                 proportion = 0.10,  # % of cells in search radius that are allowed a low score but not considered providing the service
                                 threshold = threshold)  # the size threshold set previously


   ### Clip to study area and save final raster ----
   message("Saving final and standardised scores.")
   final <- raster::writeRaster(
      raster::mask(water_score, studyArea),
      filename = file.path(save, paste(projectLog$title, runtitle, "flood_mitigation_capacity.tif", sep="_")),
      overwrite = TRUE  # perhaps not desirable but for now prevents error messages
   )

   ### Also create a standardised version

   maxval <- max(raster::values(final), na.rm = TRUE)

   final_scaled <- raster::writeRaster(
      final/maxval*100,  # rescale from 0-100
      filename = file.path(save, paste(projectLog$title, runtitle, "flood_mitigation_capacity_rescaled.tif", sep="_")),
      overwrite = TRUE  # perhaps not desirable but for now prevents error messages
   )

   timeB <- Sys.time() # stop time

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


   updateProjectLog(projectLog) # save revised log


   # Delete all the stuff we don't need anymore
   on.exit({
      rm(r, maxval)
      cleanUp(scratch)
      message("Flood risk mitigation capacity model finished. Process took ", round(difftime(timeB, timeA, units = "mins"), digits = 1), " minutes. Please check output folder for your maps.")
   })

   return({
      ## returns the objects in the global environment
      invisible({
         projectLog <<- projectLog
      })
   })

}