R/capacity_carbon.R

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

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

#' Carbon Storage Capacity Model
#'
#' Runs the carbon storage ecosystem service model, generating capacity scores based on the ability of vegetation (and top 30 cm of soil) to store carbon.

#' @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 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 (in tons C per pixel), and one rescaled 0-100 (where 100 is maximum capacity for the area).
#' @export
#'
capacity_carbon <- function(x = parent.frame()$mm,
                            studyArea = parent.frame()$studyArea,
                            res = 5,
                            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")){
      x <- do.call(rbind, x) %>% sf::st_as_sf()
   }

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

      message("Loaded hedges from ", projectLog$clean_hedges)

   }

   ### Merge the lookup table that contains carbon info
   # hab_lookup is an object built into the package

   message("Adding habitat-specific carbon storage capacity values")

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

   x <- x[c("HabCode_B", "TotCarb")]  # keep only columns we need


   ### Create raster with same properties as mastermap -----

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


   ### Rasterize -----

   # Whenever possible the rasters are saved to disk, which frees up a lot of memory

   carbon_r <- raster::writeRaster(
      fasterize::fasterize(x, r, field = "TotCarb", background = NA),
      filename = file.path(scratch, "carbon_cap"),
      overwrite = TRUE  # because it's a temporary file we don't mind overwriting it if it exists
   )

   if (use_hedges){

      # Create raster of hedge scores
      hedges_r <- raster::writeRaster(
         fasterize::fasterize(hedges, r, field = "TotCarb", background = NA),
         filename = file.path(scratch, "hedges_carbon_stor"), overwrite = TRUE)

      # Overwrite the basemap scores in places with hedges

      carbon_r <- raster::mask(carbon_r, hedges_r, inverse = T) %>% raster::cover(hedges_r)

   }

   ### Compute actual values per cell ----

   # this depends on the resolution of the raster; the value in TotCarb is tons per hectare
   # could eventually adjust by condition multiplier

   ct <- 100^2 / res^2  # this is the constant we use to convert from hectare to cell (to get tons C per cell)

   carbon_r <- raster::writeRaster(carbon_r / ct,
                           filename = file.path(scratch, "carbon_cap"),
                           overwrite = TRUE  # we want to update it
   )

   ### Clip to study area and SAVE to outputs folder
   message("Saving final and standardised scores.")

   final <- raster::writeRaster(
      raster::mask(carbon_r, studyArea),
               filename = file.path(save, paste(projectLog$title, runtitle, "carbon_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, "carbon_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_carbon"]] <- 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, ct, maxval)
      cleanUp(scratch)
      message("Carbon storage 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
      })
   })

}