R/demand_climate_reg.R

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

##################################################################
### Climate regulation demand model                            ###
### 14 October 2021                                            ###
### Sandra Angers-Blondin                                      ###
##################################################################

#' Climate Regulation Demand Model
#'
#' Runs the climate regulation ecosystem service model, generating demand scores based on two indicators: population (weighted by age risk group)  and proportion of manmade surfaces. (Specific indicators can be omitted by setting the appropriate weight to 0.)

#' @param x A basemap, in a list of sf tiles or as one sf object. Must have attributes "HabCode_B", "housePop", "riskgroup".
#' @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 local Radius (m) for focal statistics at local range. Default is 200 m.
#' @param urban Folder path to the Built-up Areas dataset (https://data.gov.uk/dataset/ad0381f2-322b-4699-94ec-32839ec02d33/built-up-areas-december-2011-boundaries-v2).
#' @param urban_size Area threshold (km2) for a builtup area to have a heat island effect (default 10 km2).
#' @param pop_density Population threshold (number of residents within "local" distance) above which we consider there is a demand (default 50).
#' @param indicators Logical; should raw indicators (before transformation into z-scores and rescaling) be saved to the project folder? Default TRUE.
#' @param indic_weights A numeric vector of length 4, with weights for distance to roads, manmade surfaces, population and health, respectively. Default to equal weights of 1 (all indicators contributing equally to final sum).
#' @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 demand scores: one with raw scores (arbitrary units), and one rescaled 0-100 (where 100 is maximum demand for the area).
#' @export
#'
demand_climate_reg <- function(x = parent.frame()$mm,
                            studyArea = parent.frame()$studyArea,
                            res = 5,
                            local = 200,
                            urban = NULL,
                            urban_size = 10,
                            pop_density = 50,
                            indicators = TRUE,
                            indic_weights = c(1, 1),
                            projectLog = parent.frame()$projectLog,
                            runtitle = parent.frame()$runtitle,
                            save = NULL){

   ## Setup ----

   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)
   }


   ## Initial checks ----

   # Is the OS vectormap path valid?
   if(!dir.exists(urban)) stop("Invalid path to builtup areas", call. = FALSE)   # error message if invalid path

   # Are vectormap data in a valid format?
   vectortype <- guessFiletype(urban)  # guess extension of OS vector to support both geopackage and shapefile # CURRENTLY ONLY SUPPORTS SHP

   # Is the indicator argument an expected value (TRUE or FALSE)?
   if(!is.logical(indicators)) {
      warning("indicators argument must be TRUE or FALSE; using TRUE by default")
      indicators <- TRUE
   }

   # Set (and create if needed) the indicator folder
   if (indicators) {
      indicator_path <- file.path(save, "indicators")  # set the indicator path

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

   # Are the weights specified correctly?
   if (is.numeric(indic_weights)){
      if (length(indic_weights) != 2){  # if only one weight specified we assign it to each layer
         message("Weights should have 2 values for this model. Defaulting to uniform weights of 1.")
         indic_weights <- c(1,1)
      }
   } else stop("Weights must be a numeric vector of length 2 for this model. Please refer to documentation.")



   # if mm is stored in list, combine all before proceeding
   if (isTRUE(class(x) == "list")){
      x <- do.call(rbind, x) %>% sf::st_as_sf()
      # NOT using rbindlist here because only keeps the extent of the first tile
   }

   # Set the attributes we need for processing
   attributes <- c("HabCode_B", "housePop", "riskgroup")

   # Are any attributes missing?
   if (any(!attributes %in% names(x))){
      stop("Could not find attribute(s) in your basemap: ",
           paste(attributes[which(!attributes %in% names(x))], collapse = ", "),
           ". Please add or rename the attributes.")
   }

   ## subset x

   x <- x[, attributes]

   # join HabClass

   x <- dplyr::left_join(x, ecoservR::hab_lookup[, c("Ph1code", "HabClass")],
                         by = c("HabCode_B" = "Ph1code"))

   #### Read in data and pre-processing ----

   urban <- loadSpatial(urban, filetype = guessFiletype(urban))
   urban <- do.call(rbind, urban) %>%   # bind in one sf object
      sf::st_geometry() %>%    # keep only the geometry (drop attributes)
      sf::st_as_sf()   # coerce to sf object

   # Transform if projection doesn't match
   urban <- checkcrs(urban, studyArea)

   # Prepare the urban layer
   urban <- sf::st_intersection(urban, studyArea) %>%  # clip to study area
      sf::st_buffer(200) %>%  # Buffer by 200 m
      sf::st_make_valid() %>% sf::st_union() %>% sf::st_as_sf() %>%   # flatten and validate
      dplyr::mutate(area = as.numeric(sf::st_area(.))/1000000) %>%  # calculate shape area, in km2
      dplyr::filter(area > urban_size)  # keep only areas larger than 10km2


   ## Create template raster

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



   ## INDICATOR 1: Population weighted by risk group ----

   if (indic_weights[1] > 0){

      message("INDICATOR 1 of 2: Population at risk")


      message("extracting houses from basemap...")

      # Create the spatial points with socioeconomic data for each household
      # We can't go straight from polygon to raster, because then a big house that spans two raster cells would
      # have its population counted twice. By converting to points first, we make sure that each house is only represented once. If more than house house falls within a raster cell, their population gets summed.

      houses <- dplyr::filter(x, HabCode_B == "J360") %>%   # filter houses from basemap
         checkgeometry("POLYGON") %>%  # need to be single-part before converting to centro
         sf::st_centroid()  # get the centroid of each house


      ## Create the population raster by summing individual house populations in each cell

      popscore <- raster::rasterize(houses,
                                    r,
                                    field = "housePop",
                                    fun = "sum", na.rm = TRUE,
                                    filename = file.path(scratch, "climreg_pop_temp"), overwrite = TRUE)


      riskscore <- raster::rasterize(houses, # risk group proportion: we want the mean here
                                     r,
                                     field = "riskgroup",
                                     fun = mean, na.rm = TRUE,
                                     filename = file.path(scratch, "climreg_risk_temp"), overwrite = TRUE)


      ## FOCAL STATS
      message("calculating population in ", local, " m radius")

      popscore <- focalScore(popscore, radius = local, type = "sum")

      ### OLD FOCAL STATS
      # Create weight matrices based on the search radius for focal stats
      # (automatically considers the res of the raster to calculate distance)
      # window_local <- raster::focalWeight(r, local, "circle")   # search window for the longer dist
      # window_local[window_local > 0] <- 1               # replacing weights by 1 (we want full sum, not mean)
      #
      #
      # message("calculating population in ", local, " m radius")
      #
      # # Calculate the summarised local scores
      # popscore <- raster::focal(popscore,
      #                           w = window_local,
      #                           na.rm = TRUE, pad = TRUE)


      # Reclassify to keep only pop densities above a defined threshold
      popscore <- raster::reclassify(popscore,
                                     cbind(0, pop_density, NA),  # remove anything below pop threshold
                                     include.lowest = TRUE)


      ## Calculate mean riskgroup in local area
      message("calculating health risk in ", local, " m radius")

      riskscore <- focalScore(riskscore, radius = local, type = "mean")

      ### OLD FOCAL STATS
      # riskscore <- raster::focal(riskscore,
      #                            w = window_local,
      #                            fun = function(x){mean(x[window_local != 0], na.rm=TRUE)
      #                               # we select only the values in the circular window (discarding values where window is 0), and perform a mean calculation on this subset rather than the full rectangular extent of the matrix. This ensure the correct denominator while discarding NAs.
      #                            })
      # # custom function defined, as default mean method of focal is flawed with circular window


      ### Multiply them and mask by builtup area
      popscore <- raster::mask(popscore * riskscore,
                               urban)


      rm(riskscore)  # remove as no longer needed


      # ## The score is then masked to show only those areas within 250m of land or houses.
      # ## Mask everything that is further than 250m from land or houses (basically water....)
      # ## currently not applied as too difficult / inefficient
      #
      # housemask <- sf::st_buffer(houses, 250) %>%
      #   checkgeometry() %>%  # must be valid singlepart polygons to rasterize properly
      #   fasterize::fasterize(., r, field = NULL)
      #
      # ## Would be really intensive on a big basemap... can we find a workaround?
      # landmask <- dplyr::filter(x,
      #                           !grepl(paste0(c("Water", "Sea"), collapse = "|"),
      #                                 HabClass,
      #                                 ignore.case = TRUE)) %>%
      #   checkgeometry() %>%  # must be valid singlepart polygons to rasterize properly
      #   fasterize::fasterize(., r, field = NULL)
      #
      #
      # # combine house and land mask and apply



      # Save raw indicator if need be

      if (indicators){
         raster::writeRaster(popscore, file.path(indicator_path,   # we define the path
                                                 paste(
                                                    projectLog$title, runtitle, "clim_reg_pop_indic.tif", sep = "_")),
                             overwrite = TRUE)
      }


   } else {

      popscore <- r # if user wishes to omit this indicator we create empty raster

   }




   ## INDICATOR 2: Manmade surfaces ----

   if (indic_weights[2] > 0){

      message("INDICATOR 2 of 2: manmade surfaces")

      # Proportion (%) of man‐made surface within search distance- higher proportions are correlated with higher air pollution

      # Subset the basemap
      manmade <- dplyr::filter(x, HabClass %in% c("Urban", "Infrastructure"))

      # Make sure we're dealing with polygons
      manmade <- checkgeometry(manmade, "POLYGON")

      # Rasterize
      manmade <- fasterize::fasterize(manmade,
                                      r,
                                      field = NULL)   # manmade features get value of 1


      ### Focal stats on proportion of manmade surfaces
      manmadescore <- focalScore(manmade, radius = local, type = "cover")

      #### OLD FOCAL STATS
      ### Focal stats for proportion of surfaces (because they have value of 1,
      # sum will give the number of cells with manmade features in the radius,
      # and we divide by total ncell for proportion)

      # # Create weight matrices based on the search radius for focal stats
      # # (automatically considers the res of the raster to calculate distance)
      #
      # window_local <- raster::focalWeight(r,
      #                                     local,
      #                                     "circle")   # search window
      # window_local[window_local > 0] <- 1     # replacing weights by 1
      #
      # denominator <- raster::ncell(window_local[window_local>0])  # for custom mean function
      #
      # # We use the default focal stats (sum) function which is way faster than a custom function, THEN we divide the results by the denominator (number of cells in the window) to get the proportion
      # manmade <- raster::focal(manmade, w = window_local, na.rm = TRUE, pad = TRUE)
      # manmadescore <- manmade/denominator

      # Save output: this output has values 0-1
      if (indicators){
         raster::writeRaster(manmadescore,
                             filename = file.path(indicator_path,
                                                  paste(projectLog$title, runtitle, "clim_reg_manmade_indic.tif", sep="_")),
                             overwrite = TRUE)

         message("Manmade surfaces indicator saved.")
      }

      rm(manmade)


   } else {

      manmadescore <- r
   }



   # RESCALE INDICATORS ------------------------------------------------------

   # Each indicator gets transformed into a z-score (mean-centered and deviance-scaled), and then reprojected on a 1-100 scale

   manmadescore <- rescale_indicator(manmadescore)
   popscore <- rescale_indicator(popscore)


   # COMBINE AND SAVE  -------------------------------------------------------

   # All indicators should now be between 0-1 and we can sum them.

   message("Combining indicators...")
   final <- sum(
      indic_weights[1]*popscore,
      indic_weights[2]*manmadescore,
      na.rm = TRUE  # important otherwise only gives raster for where all 4 indicators have values
   )


   rm(manmadescore, popscore)

   # Replace NA with 0
   final[is.na(final)] <- 0

   # Mask to study area

   final <- raster::mask(final, studyArea)

   message("Saving final scores.")

   raster::writeRaster(final,
                       filename = file.path(save, paste(projectLog$title, runtitle,
                                                        "climate_reg_demand.tif",
                                                        sep="_")),
                       overwrite = TRUE)

   ## and a rescaled version

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

   raster::writeRaster(final/maxvalue *100,
                       filename = file.path(save, paste(projectLog$title, runtitle,
                                                        "climate_reg_demand_rescaled.tif", sep="_")),
                       overwrite = TRUE)



   timeB <- Sys.time() # start time


   # write performance to log
   projectLog$performance[["dem_clim"]] <- 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, maxvalue)
      cleanUp(scratch)
      message("Climate regulation demand 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
      })
   })

}