R/demand_air.R

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

##################################################################
### Air purification demand model                              ###
### 12 October 2021                                            ###
### Karl Busdieker and Sandra Angers-Blondin                   ###
##################################################################

#' Air Purification Demand Model
#'
#' Runs the air purification ecosystem service model, generating demand scores based on four indicators: distance to major roads, proportion of manmade surfaces, 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".
#' @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 300 m.
#' @param openroads Path to the folder containing OS Open Roads data or other (linestring) road classification to use as source of pollution.
#' @param roadClassification The attribute name in OS Open Roads containing the classification of road types, default "roadClassification". Only "Motorway" and "A Road" types are selected for analysis. (If using custom road network, please ensure classification matches these terms exactly.)
#' @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_air_purif <- function(x = parent.frame()$mm,
                             studyArea = parent.frame()$studyArea,
                             res = 5,
                             local = 300,
                             openroads = NULL,
                             roadClassification = "roadClassification",
                             indicators = TRUE,
                             indic_weights = c(1,1,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(openroads)) stop("Invalid OS OpenRoads path", call. = FALSE)   # error message if invalid path

   # Are vectormap data in a valid format?
   vectortype <- guessFiletype(openroads)  # 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) != 4){  # if only one weight specified we assign it to each layer
         message("Weights should have 4 values for this model. Defaulting to uniform weights of 1.")
         indic_weights <- c(1,1,1,1)
      }
   } else stop("Weights must be a numeric vector of length 4 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")

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

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

   # NOTE - Distance from roads is assumed to correlate with higher levels of air pollution.
   # Proximity score to nearest busy road (high score = close distance)

   message("INDICATOR 1 of 4: Distance from roads")

   ## Extract roads from basemap
   mmroads <- dplyr::filter(x, HabCode_B == "J511")


   ## Import OS OpenRoads data
   message("Importing OS Open Roads")
   ## Extract only layers we need
   oproads_files <- list.files(path = openroads, pattern = vectortype, recursive = TRUE, full.names = TRUE)   # list of all files in folder


   ## FRAGILE: IMPLEMENT BETTER HANDLING (in case multilinestring etc)

   osroads <- lapply(oproads_files, function(x) {
      # read in the file(s)
      sf::st_read(x,
                  layer = sf::st_layers(oproads_files)$name[sf::st_layers(oproads_files)$geomtype == "Line String"]) %>%
         dplyr::select(classif = !!roadClassification) %>%  # rename classification attribute to set name
         dplyr::filter(classif %in% c("Motorway", "A Road"))  ## Keep only roads of interest
   })


  ## Check the crs

   osroads <- checkcrs(osroads)

   osroads <- do.call(rbind, osroads)


   ## clip to study area

   index <- sf::st_intersects(osroads, SAbuffer) # index for faster clipping

   osroads <- sf::st_intersection(osroads[lengths(index) > 0,], SAbuffer) %>% sf::st_make_valid() # clip


   # buffer the major roads by 40 m
   major_roads40m <- sf::st_buffer(osroads, 40) %>%
      checkcrs(SAbuffer) %>%
      checkgeometry("POLYGON")

   rm(osroads)

   # MasterMap roads are more accurate: we only use OpenRoads classification to identify the major roads in mastermap


   message("Identifying major roads")
   index <- sf::st_intersects(mmroads, major_roads40m)
   mmroads <- mmroads[lengths(index) > 0, ]

   mmroads <- sf::st_intersection(mmroads, major_roads40m) %>% sf::st_make_valid()  # clip because the index also selects minor roads which veer off the major ones

   rm(major_roads40m)



   ## Calculate distance to roads

   message("Computing distance to roads")
   # Make sure roads are polygons
   mmroads <- checkgeometry(mmroads, "POLYGON")


   # calculate distance to roads, writing resulting rasters to scratch folder:
   # (using same function we use in noise model)
   # highest score of 1 on roads, down to 0 at "local" m distance
   roadscore <- distance_from_source(mmroads, r, "air_demand_roads_score", local)

   ## This was actual distances but the score must be inversely proportional.
   ## From ecoservGIS:converted to log10 score and then to inverse score

   if (terra::hasValues(roadscore)){
      # log10 of distance 0 is an error so we add 1 to everything
      roadscore <- 1 - (log10(roadscore + 1) / log10(local))

      # score could be slightly negative (precision), so forcing to zero
      roadscore[roadscore < 0] <- 0

   } else {
      # in some cases there are no major roads in the study area and the outcome of distance_from_source is an empty raster
      terra::values(roadscore) <- 0
   }



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

      message("Road distance indicator saved.")
   }

   rm(mmroads)


   } else {

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


   # ## Produce a roadmask as we only consider demand when it is close to roads
   #
   # if (indic_weights[1] > 0){
   # roadmask <- roadscore  # copy roadscore
   # roadmask[!is.na(roadmask)] <- 1   # transform all values to 1
   #
   # } else {
   #
   #    roadmask <- r
   #    roadmask[] <- 1 # if not considering roads then we don't want a mask
   # }




   # Manmade surfaces indicator ----------------------------------------------

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

   message("INDICATOR 2 of 4: 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 for proportion of surfaces

   manmadescore <- focalScore(manmade, radius = local, type = "cover")


   #### OLD FOCAL STATS
   # because raster has values 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




   ## Mask to only keep what is close to roads: should be there but omitted as buggy?

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


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

      message("Manmade surfaces indicator saved.")
   }

   rm(manmade)

   } else {

      manmadescore <- r
   }


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


   if (indic_weights[3] > 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 3 of 4: 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

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


   # Focal stats (sum)

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

   popscore <- focalScore(population, 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(population, 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, "air_purif_pop_indic.tif", sep = "_")),
                          overwrite = TRUE)
   }

   } else {
      popscore <- r
   }

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

   if (indic_weights[4] > 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 4 of 4: health score")

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

   rm(houses)

   ## Focal statistics
   healthscore <- focalScore(healthmap, radius = local, type = "mean")

   ### OLD FOCAL STATISTICS
   # # Calculate the summarised local scores
   #
   # healthscore <- raster::focal(healthmap, 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, "air_purif_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

   roadscore <- rescale_indicator(roadscore)
   manmadescore <- rescale_indicator(manmadescore)
   popscore <- 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]*roadscore,
      indic_weights[2]*manmadescore,
      indic_weights[3]*popscore,
      indic_weights[4]*healthscore,
      na.rm = TRUE  # important otherwise only gives raster for where all 4 indicators have values
   )


   rm(roadscore, manmadescore, popscore, 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,
                                                        "air_purif_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,
                                                        "air_purif_demand_rescaled.tif", sep="_")),
                       overwrite = TRUE)



   timeB <- Sys.time() # start time


   # write performance to log
   projectLog$performance[["dem_air"]] <- 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("Air purification 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
      })
   })

}