R/capacity_noise_reg.R

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

########################################
### Noise regulation capacity model  ###
### for EcoservR tool                ###
### Sandra Angers-Blondin            ###
### 18 Dec 2020                      ###
########################################

#' Noise Regulation Capacity Model
#'
#' Runs the noise regulation ecosystem service model, generating capacity scores based on the ability of vegetation to reduce noise pollution.

#' @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 res Desired resolution of the raster. Default is 5 m. Range recommended is 5-10m.
#' @param short Radius (m) for focal statistics at short range. Default is 30 m.
#' @param local Radius (m) for focal statistics at local range (maximum distance for effectiveness). Default is 300 m.
#' @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_noise_reg <- function(x = parent.frame()$mm,
                               studyArea = parent.frame()$studyArea,
                               res = 5,
                               short = 30, local = 300,
                               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()
      # NOT using rbindlist here because only keeps the extent of the first tile
   }

   studyArea <- sf::st_zm(studyArea, drop=TRUE)

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

   ### 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') %>%
         merge(hab_lookup[c("Ph1code", "Noise")], by.x = 'HabCode_B', by.y = 'Ph1code', all.x = TRUE)

      message("Loaded hedges from ", projectLog$clean_hedges)
      hedges <- checkcrs(hedges, 27700)

   }


   ### Merge the lookup table -----

   x <- merge(x, hab_lookup, by.x = "HabCode_B", by.y = "Ph1code", all.x = TRUE)

   x <- x[c("HabCode_B", "Noise")] # keep only required columns


   ### 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 -----

   noise_r <- raster::writeRaster(
      fasterize::fasterize(x, r, field = "Noise", background = NA),
      filename = file.path(scratch, "noisereg_score"),
      overwrite = TRUE
   )


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

      noise_r <- raster::mask(noise_r, hedges_noise_r, inverse = T) %>%
         raster::cover(hedges_noise_r)

   }

   ### Focal statistics ----

   message("Calculating noise regulation score at short range")
   noise_short <- focalScore(noise_r, radius = short, type = "sum")

   message("Calculating noise regulation score at local scale")
   noise_long <- focalScore(noise_r, radius = local, type = "sum")



   rm(noise_r)  # remove old raster to clear up some memory

   # Sum the two to create the score layer
   noise_score <- raster::writeRaster(
      sum(noise_short, noise_long, na.rm = TRUE),
      filename = file.path(scratch, "noisereg_score"),
      overwrite = TRUE)

   rm(noise_short, noise_long) # remove old objects


   ### Apply functional threshold -----
   message("Applying functional threshold")
   ## Patches smaller than a certain size are considered not useful and removed

   noise_score <- raster::writeRaster(
      functionalMask(noise_score,  # the score raster we apply the function to
                     local = local,  # the local (wide) search radius set previously
                     res = res,   # the resolution of the raster set previously
                     proportion = 0.10,  # poportion of cells in the search circle which can have a score of 10 without being considered to contribute to the service
                     threshold = threshold)  # the size threshold set previously
      , filename = file.path(scratch, "noisereg_score"),
      overwrite = TRUE)


   ### Clip to study area and save final file
   message("Saving final and standardised scores.")

   final <- raster::writeRaster(
      raster::mask(noise_score, studyArea),
      filename = file.path(save,
                           paste(projectLog$title, runtitle, "noise_regulation_capacity.tif", sep="_")),
      overwrite = TRUE
   )
   rm(noise_score)

   ### 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, "noise_regulation_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_noise"]] <- 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, final, final_scaled, maxval)
      cleanUp(scratch)
      message("Noise regulation 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
      })
   })


}