R/demand_access_nature.R

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

##################################################################
### Accessible nature demand model                             ###
### 01 March 2022                                              ###
### Sandra Angers-Blondin                                      ###
##################################################################

#' Access to Nature Demand Model
#'
#' Runs the access to nature ecosystem service model, generating demand scores based on three indicators: proportion of households meeting the criteria of having accessible greenspaces >2ha within a 300m radius (or specified by user with the "local" argument), population, and health. (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", "health", "GI", and "GIpublic".
#' @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, and threshold for distance to greenspace. Default is 300 m.
#' @param pop_threshold Population density (people per hectare) below which to consider there is no demand. Default is 5.
#' @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 3, with weights for distance to greenspaces, population and health, respectively. Default to c(1, 0.5, 0.5) (so that population and health combined contribute the same amount as the ANGSt indicator).
#' @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_access_nature <- function(x = parent.frame()$mm,
                                 studyArea = parent.frame()$studyArea,
                                 res = 5,
                                 local = 300,
                                 pop_threshold = 5,
                                 indicators = TRUE,
                                 indic_weights = c(1,0.5,0.5),
                                 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 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) != 3){  # if only one weight specified we assign it to each layer
         message("Weights should have 3 values for this model. Defaulting to uniform weights of 1.")
         indic_weights <- c(1,1,1)
      }
   } else stop("Weights must be a numeric vector of length 3 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", "health", "GIpublic", "GI")

   # 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.")
   }



   # Pre-processing ----------------------------------------------------------

   # Keep only fields we need
   x <- x[attributes]

   # join HabClass

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


   # Check that the study area is in the right projection (OS national grid). If not it will be transformed, and in either case we then force it to avoid GDAL system-specific "no colon in init string" error
   studyArea <- suppressWarnings({
      checkcrs(studyArea, 27700) %>%  # make sure study area is 27700, reproject if needed
         sf::st_set_crs(27700)  # and set it manually now we know it's right
   })


   # Create a buffer around the study area (to clip roads with some margin)
   SAbuffer <- sf::st_buffer(studyArea, 500) %>% checkcrs(27700)


   # Create raster with same properties as mastermap to use as a template

   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


   # Indicator 1: distance to greenspaces ------------------------------------------

   if (indic_weights[1] > 0){  # only if user wants this indicator

      message("INDICATOR 1 of 3: Distance from greenspaces")


      message("Extracting greenspaces from basemap")
      ### Extract GI from basemap

      gi <- x %>%
         dplyr::filter(GIpublic %in% c("Public", "Restricted") | (is.na(GIpublic) & GI %in% c("Undetermined Greenspace",
                                                       "Undertermined Greenspace")) # typo fix for legacy maps
                                                                  )

      ## buffer to join up bits separated by paths etc

      gi <- gi %>%
         sf::st_buffer(2.5) %>%
         checkgeometry()

      ## rasterize and then polygonize (faster than dissolving vector data)

      gi_r <- fasterize::fasterize(gi, r)

      gi_poly <- sf::st_as_sf(
         stars::st_as_stars(gi_r),
         merge = TRUE) %>%
         dplyr::mutate(area_ha = as.numeric(sf::st_area(.))/10000)  # calculate area for filtering

      ## extract the GI meeting size requirements
      gi2ha <- gi_poly %>% dplyr::filter(area_ha >= 2)   # woods greater than 2ha

      gi2ha <- checkgeometry(gi2ha)



      ### Extract houses, with the TOID to link up values if needed
      houses <- x %>%
         dplyr::filter(HabCode_B == "J360") %>%
         sf::st_centroid() %>%
         sf::st_intersection(sf::st_geometry(studyArea))  # just those in study area

      message("Calculating distances to houses")

      ### Calculate the indicator: greenspace > 2ha within 300m of house
      ## Reverse process: buffer GI by the ANGST distance threshold and then identify houses within

      gi_local <- sf::st_buffer(gi2ha, local)   # buffer the gi

      # Find out which houses intersect this buffer: assign them values of 0 or 1 accordingly

      index_house <- sf::st_intersects(houses, gi_local)

      # create new pass/fail attribute for houses

      houses$angst <- 0 # initiate col
      houses[lengths(index_house) > 0, ]$angst <- 1  # add a pass score for the houses within buffer


      ## Now rasterize with SUM

      house_r <- terra::rasterize(houses, r, field = "angst", fun = "sum")  # number of houses that pass
      denom_r <- terra::rasterize(houses, r, field = 1, fun=sum)            # count of all houses

      ## Proportion of failed in focal window
      message("Applying ANGST criteria")

      # Set circular window
      fw <- raster::focalWeight(r, local, "circle") #The equivalent Terra function is focalMat
      fw[fw > 0] <- 1   # replacing weights by 1

      #using the Terra version of the focal function
      pass <- terra::focal(x = house_r, w = fw, fun = "sum", na.rm = TRUE)
      denom<- terra::focal(x = denom_r, w = fw, fun = "sum", na.rm = TRUE)

      angst <- 1- pass/denom



      # save raw indicator
      if (indicators){
         angst <- raster::writeRaster(angst,
                                      filename = file.path(indicator_path,
                                                           paste(projectLog$title,runtitle, "access_nature_angst_indic.tif", sep="_")),
                                      overwrite = TRUE)

         message("Distance to greenspace indicator saved.")
      }



   } else {

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




   # Population raster -------------------------------------------------------


   if (indic_weights[2] > 0){  # only if user wants this indicator

      # Population size (sum) within local search distance
      # NOTE - Higher population scores mean there is more societal need

      message("INDICATOR 2 of 3: population density")

      # 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, "access_nature_pop_temp"), overwrite = TRUE)



      message("...calculating population in ", local, " m radius")

      # Calculate the summarised local scores
      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)



      ## Mask to only keep what is close to roads
      #popscore <- raster::mask(popscore, roadmask)

      # Save raw indicator if need be

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

   } else {
      popscore <- r
   }

   # Health indicator --------------------------------------------------------

   if (indic_weights[3] > 0){  # only if user wants this indicator

      # Health deprivation score (mean) within local search distance
      # NOTE - Higher health scores mean an area is more deprived, i.e. there is more societal need for noise regulation
      # the HDD score is centered around 0 so there are negative values which need taking care of

      message("INDICATOR 3 of 3: health score")

      # Create the health raster from the houses data we subsetted in previous indicator
      healthscore <- raster::rasterize(houses, r, field = "health", fun = mean, na.rm = TRUE,
                                       filename = file.path(scratch, "access_nature_health_temp"),
                                       overwrite = TRUE)

      rm(houses)

      # Calculate the summarised local scores
      healthscore <- focalScore(healthscore, radius = local, type = "mean")


      ## OLD FOCAL STATS
      # healthscore <- raster::focal(healthscore, 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.
      #                              },
      #                              pad = TRUE) # pad edges with NAs so raster has full extent)


      ## Mask to only keep what is close to roads

      # healthscore <- raster::mask(healthscore, roadmask)

      # Save raw indicator if need be

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




   } else {

      healthscore <- r
   }

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

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

   angst <- rescale_indicator(angst)
   popscore2 <- rescale_indicator(popscore)
   healthscore <- rescale_indicator(healthscore)

   # COMBINE INDICATORS ------------------------------------------------------

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

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


   ## Apply a population threshold

   # we work out what the score is for a focal sum if the average pop per hectare is what is specified in the function
   # people per hectare times the area in hectares of the focal circle

   popvalue_thres <- pop_threshold * length(fw[fw > 0])*res^2 / 10000

   ## Mask the raster

   popscore[popscore <= popvalue_thres] <- NA

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

   rm(angst, popscore, popscore2, healthscore)

   # 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,
                                                        "access_nature_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,
                                                        "access_nature_demand_rescaled.tif", sep="_")),
                       overwrite = TRUE)



   timeB <- Sys.time() # start time


   # write performance to log
   projectLog$performance[["dem_access"]] <- 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("Access to nature 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
      })
   })

}