R/capacity_pollination.R

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

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


#' Pollination Capacity Model
#'
#' Runs the pollination ecosystem service model, generating capacity scores based on the suitability of semi-natural habitats to support insect pollinators and travel distances of pollinators.

#' @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 dist Distance within which pollinators are found in edge habitats (woodland) from a core habitat. Default 20 m.
#' @param elev A character string giving the name of the basemap attribute containing a polygon's mean elevation. The output raster will be masked to reflect the fact that pollinators are not abundant at high elevations. If the information has not been extracted in the basemap, a DTM can be supplied separately.
#' @param DTM A raster of elevation values for the study area (specified as file path to a folder). If left NULL (default), will use the basemap "elev" information instead (quicker).
#' @param elev_threshold The elevation (in m) threshold above which the model considers that pollinators are absent. Default 250 m.
#' @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 (0-1: likelihood that a pollinator will visit a pixel), and one rescaled 0-100 (where 100 is maximum capacity for the area).
#' @export
#'
capacity_pollination <- function(x = parent.frame()$mm,
                                 studyArea = parent.frame()$studyArea,
                                 res = 5,
                               dist = 20,
                               elev = NULL,
                               DTM = NULL,
                               elev_threshold = 250,
                               use_hedges = FALSE,
                               projectLog = parent.frame()$projectLog,
                               runtitle = parent.frame()$runtitle,
                               save = NULL
){

   timeA <- Sys.time() # start time

   ## issue a warning if no elevation has been specified:

   if (is.null(DTM) && is.null(elev)){
      message("Warning! You have not specified an elevation attribute or a DTM. The model will run but the output will not consider the effect of elevation on pollinator distribution.")
   }


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

   ### If user chooses to use DTM, check and import

   if (!is.null(DTM)){

      dtm <- DTM

      if(grepl(pattern = paste0(c(".tif", ".asc", ".gr"), collapse = "|"), DTM, ignore.case = TRUE)){

         dtm <- dirname(DTM)  # if user gave full file name, we cut back to folder
      }

      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",ignore.case = TRUE)) == 0 &&   # any tif?
       length(list.files(path = dtm, pattern = "asc",ignore.case = 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))

   # 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 raster 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 <- do.call(raster::merge, dtm)
      message("Merged DTM saved to: ", dtmpath, " .You can use it to speed up future runs of the model.")
   }

   } # end of dtm import and processing

   if (!is.null(elev) && !elev %in% names(x)) stop(paste0("Elevation attribute \"", elev, "\" not found in your basemap. Please specify attribute name with the \"elev\" argument, or set elev = NULL to ignore elevation."))



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

   ### 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[c("Ph1code",'HabBroad', 'HabClass')], 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"))
   }

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

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

   x$HabBroad <- NULL
   x$HabClass <- NULL

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

   if (!is.null(elev)){
      x <- x[c("HabCode_B", "HabBroad", "HabClass", elev)]  # keep only required columns
   } else {
      x <- x[c("HabCode_B", "HabBroad", "HabClass")]  # 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


   ### Create core habitat layer -----
   message("Creating layer of core habitats")
   # list of habitats that support pollinators
   corehabs <- c("J21", "J22", "J23", "J12v", "Linear",
                 "H65", "H66", "H84", "H85",
                 "H24", "H26", "H2u",
                 "GW", "GR", # green wall / green roof interventions
                 "Gardens / Parks / Brownfield",
                 "Garden",
                 "Grassland, semi-natural",
                 "Heathland",
                 "Scrub")

   core <- dplyr::filter(x, HabCode_B %in% corehabs | HabBroad %in% corehabs)

   # Create buffer around core with distance specified
   corebuffer <- sf::st_buffer(core, dist) %>%
      sf::st_make_valid() %>% sf::st_geometry() %>% sf::st_as_sf() %>%
      sf::st_cast("MULTIPOLYGON") %>% sf::st_cast("POLYGON")

   if (use_hedges){
      # add hedgerows to the corebuffer object
      hedges <- sf::st_buffer(hedges, dist) %>%
         sf::st_make_valid() %>% sf::st_geometry() %>% sf::st_as_sf() %>%
         sf::st_cast("MULTIPOLYGON") %>% sf::st_cast("POLYGON")

      hedges <- rename_geometry(hedges, attr(corebuffer, "sf_column"))
      hedges <- checkcrs(hedges, 27700)
      corebuffer <- rbind(corebuffer, hedges)
   }


   ### Create edge habitat layer -----
   message("Creating layer of edge habitats")

   edgehabs <- c("Woodland, coniferous",
                 "Woodland, broadleaved",
                 "Woodland, mixed")    # selection of woodlands as edge habitats

   edge <- dplyr::filter(x, HabBroad %in% edgehabs)

   # use buffer around core habitats to clip edge habitats used by pollinators
   edgeclip <- sf::st_intersection(edge, corebuffer) %>% sf::st_make_valid()
   edgeclip <- checkgeometry(edgeclip, "POLYGON")


   ### Merge core and edge habitats -----
   cols <- intersect(colnames(edgeclip), colnames(core))  # get index of shared columns (so we don't keep the duplicated columns from the intersection)

   # bind all together and tidy up
   fullhab <- rbind(core[,cols], edgeclip[,cols])  %>%
      checkgeometry(., "POLYGON")
      # st_make_valid() %>%
      # st_cast(to = "MULTIPOLYGON") %>%  # equivalent of multi-part to single part
      # st_cast(to = "POLYGON", warn = FALSE)           # equivalent of multi-part to single part


   rm(core, corebuffer, corehabs, edge, edgeclip, edgehabs, cols)



   ### Rasterize -----

   # creates a raster where habitat has value of 1 and the rest 8888 (8888 is essential non-habitat code flag for the calculateDistances function)

   pollin_r <- raster::writeRaster(
      fasterize::fasterize(fullhab, r, field = NULL, background = 8888),
      filename = file.path(scratch, "pollination_habs"),
      overwrite = TRUE  # because it's a temporary file we don't mind overwriting it if it exists
   )

   if (use_hedges){

      # add hedge habitats to the pollin raster
      hedges_r <- raster::writeRaster(
         fasterize::fasterize(hedges, r, background = NA),
         filename = file.path(scratch, "hedges_r"), overwrite = TRUE)

      # Overwrite the basemap
      pollin_r <- raster::mask(pollin_r, hedges_r, inverse = T) %>%
         raster::cover(hedges_r)
   }

   ### Mask by elevation -----

   ## Mask by elevation (if elev attribute supplied)
   if (!is.null(elev)){

      message("Masking by elevation")

      #elevmask <- fasterize::fasterize(x[x[[elev]] < elev_threshold,], r)


      elevmask <- fasterize::fasterize(x, r, field = elev)  # create raster of elevation

      ## gap fill elevation for NA cells only
      elevmask <- raster::focal(elevmask, w=matrix(1,nrow=5,ncol=5), fun = mean, na.rm = TRUE, NAonly = TRUE)

      elevmask <- raster::reclassify(elevmask,
                                     cbind(elev_threshold, Inf, NA), right=FALSE)

      pollin_r <- raster::mask(pollin_r,
                               elevmask)

   } else if (!is.null(DTM)){

      message("Masking by elevation")

      # Reclassify dtm
      dtm[dtm < elev_threshold] <- 1
      dtm[dtm >= elev_threshold] <- NA

      # Make sure dtm aligns with indicator raster
      if (!identical(raster::extent(dtm), raster::extent(r)) |
          !identical(raster::crs(dtm), raster::crs(r))|
          !identical(raster::res(dtm), raster::res(r))){
         dtm <- raster::projectRaster(dtm, r)
      }

      pollin_r <- raster::mask(pollin_r, dtm)
      pollin_r[is.na(pollin_r)] <- 8888 # for calculatedistance function
   }


   ### Calculate distances ----

   # Buffer around the area and simplify (used as mask to reduce computation time of distances)

   mask <- sf::st_geometry(fullhab) %>%  # drop all columns except geometry for faster and more efficient handling
      sf::st_buffer(668) %>%  # draws buffer around the selected habitats, default 668m is the distance away from a habitat where not likely to find pollinators
      sf::st_union() %>%   # union to dissolve overlaps
      sf::st_sf()  # return to sf object

   # Mask the raster: we're effectively setting the max search radius for the next steps
   pollin_r <- raster::writeRaster(
      raster::mask(pollin_r, mask),
      filename = file.path(scratch, "pollination_habs"),
      overwrite = TRUE  # because it's a temporary file we don't mind overwriting it if it exists
   )

   rm(mask, fullhab)

   message("Calculating distances...")

   ## Now we have a much reduced object to convert to spatstats object to calculate distances
   # All the NAs (cells outside search area) can be safely ignored, we're only interested in the 1 (source) and 8888 (to calculate distances from source)

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

   # Because Euclidean distances can still be slightly greater than the size of the buffer, we set all distances above threshold to NA
   pollin_r <- raster::reclassify(pollin_r,
                          rcl = c(668,Inf,NA),
                          filename = file.path(scratch, "pollination_score2"),
                          overwrite = TRUE)

   message("Distance raster created. Calculating likelihood of pollinator visit...")

   ### Calculate likelihoods -----
   pollin_r <- raster::writeRaster(
      exp(pollin_r*(-0.00104)),
      filename = file.path(scratch, "pollination_score3"),
      overwrite = TRUE  # because it's a temporary file we don't mind overwriting it if it exists
   )

   pollin_r <- raster::reclassify(pollin_r, cbind(NA, NA, 0), right=FALSE)  # replace NA with 0


   ### Clip to study area and save to outputs folder -----

   message("Saving final and standardised scores.")

   final <- raster::writeRaster(
      raster::mask(pollin_r, studyArea),
      filename = file.path(save, paste(projectLog$title, runtitle, "pollination_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, "pollination_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_pollin"]] <- 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, pollin_r, maxval)
      cleanUp(scratch)
      message("Pollination 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
      })
   })


}