ecoserv_tool: Natural Capital Mapping Toolkit (UK)

Documented in create_network

###########################################
### Ecological networks                 ###
### for EcoservR tool                   ###
### Sandra Angers-Blondin               ###
### 12 May 2021                         ###
###########################################

#' Ecological Networks
#'
#' Maps an ecological network to inform connectivity analysis. Six different habitat types can be mapped: woodlands, grasslands, heathlands, lowland wetlands, mires (bogs), and ponds. A generic focal species approach is used and the resulting network is therefore dependent on the habitat requirements (patch size) and dispersal ability of this hypothetical representative species.

#' @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.
#' @param habitat The habitat for which to compute a network. Must be exactly one of the following: "woodland", "grassland", "heath", "pond", "mire", "wetland"
#' @param minsource The minimum patch size considered to be a viable habitat for the generic focal species, in hectares. Default is 0.1 ha (1000 square meters)
#' @param dispersal The dispersal distance (in meters) of generic focal species at three scales: core, meso and macro. Do not set the parameter if you wish to use the sensible defaults (see user guide), or specify as a vector of three numeric values, e.g. c(150, 1000, 1500)
#' @param threshold Functional size threshold (in hectares) under which a mapped network is considered too small to serve its purpose. Default is 1 ha.
#' @param use_hedges Use a separate hedgerow layer? Default FALSE, see clean_hedgerows() for producing a model-ready hedge layer.
#' @param res Desired resolution of the raster. Default is 5 m. Range recommended is 5-10m.
#' @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 "networks_". Do not use this argument unless you need to save the outputs somewhere else.
#' @param export_raster Logical. Save a raster version of the network as well as a polygon layer? Default TRUE.
#' @return A polygon and optionally a raster identifying the core, meso, and macro networks with values of 1, 2, and 3 respectively. Note that the scales build on each other (e.g. meso network is the combination of all 1 and 2 values.)
#' @export
#'
create_network <- function(x = parent.frame()$mm,
                           studyArea = parent.frame()$studyArea,
                           habitat,  # the type of habitat (woodland, mire, wet grassland, heath, grassland, pond)
                           minsource = 0.1, # minimum patch size to be considered source (in ha). Default 0.1 ha (1000 square meters)
                           dispersal = NULL, # dispersal distance (m) for the generic focal species at the core, meso and macro scale. Leave NULL to use the default values, or use a vector of three numeric values.
                           threshold = 1, # minimum patch size (ha) to be considered a functional network
                           use_hedges = FALSE,
                           res = 5, # resolution
                           projectLog = parent.frame()$projectLog,
                           runtitle = parent.frame()$runtitle,
                           export_raster = TRUE,
                           save = NULL
){


   # Create output directory automatically if doesn't already exist
   if (is.null(save)){

      save <- file.path(projectLog$projpath,
                        paste0("networks_", 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()
   }


   # Check that parameters are correctly specified

   if (!habitat %in% c("woodland", "grassland", "heath", "pond", "mire", "wetland")){
      stop("Habitat must be one of: woodland, grassland, heath, wetland, mire, or pond")
   }

   if(!is.null(dispersal) & !(length(dispersal) == 3) & inherits(dispersal, "numeric")){
      stop("Parameter \"dispersal\" must be a vector of 3 numeric values (in meters).
           Refer to user guide or do not set the parameter to use default values.")
   } else if (is.null(dispersal)){

      dispersal <- if (habitat == "woodland"){c(150, 1000, 1500)} else if
      (habitat == "grassland"){c(150, 1000, 1500)} else if
      (habitat == "heath"){c(250, 1000, 1500)} else if
      (habitat == "pond"){c(50, 500, 1500)} else if
      (habitat == "mire"){c(150, 2000, 3000)} else if
      (habitat == "wetland"){c(150, 1000, 1500)}

   }


   # Identify which lookup column is relevant
   costcol <- if (habitat == "woodland"){"CostWood"} else if
   (habitat == "grassland"){"CostGrass"} else if
   (habitat == "heath"){"CostHeath"} else if
   (habitat == "pond"){"CostPond"} else if
   (habitat == "mire"){"CostMire"} else if
   (habitat == "wetland"){"CostLWet"}

   studyArea <- sf::st_zm(studyArea, drop=TRUE) # make sure study area doesn't have Z dim

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

   ### 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", costcol)], by.x = 'HabCode_B', by.y = 'Ph1code', all.x = TRUE)

      message("Loaded hedges from ", projectLog$clean_hedges)
      hedges <- checkcrs(hedges, 27700)

   }


   # Join costs to the basemap keeping only relevant attributes for faster processing
   x <- dplyr::left_join(x[, c("HabCode_B")],
                         ecoservR::hab_lookup[,c("Ph1code", "HabBroad", costcol)], # join the relevant column
                         by = c("HabCode_B" = "Ph1code"))


   # # Check geometry of mm as needs to be clean for rasterizing
   # x <- ecoservR::checkgeometry(x, "POLYGON")

   # Create raster template for all rasters

   r <- raster::raster()
   raster::crs(r) <- sp::CRS(SRS_string = "EPSG:27700")
   raster::extent(r) <- raster::extent(x)
   raster::res(r) <- res


   message("Identifying source habitats for ", habitat)

   # Identify patches of source habitat

   if (habitat == "pond"){
      source <- x
      } else {
   source <- x[x[[costcol]] == 1, ]
      }

   ## Refine further, as cost of 1 is not necessarily the target habitat

   #### TO DO: CONDITIONAL FILTERING BY HABITAT HERE; ESP FOR POND, MIRE AND LWET

   if (habitat == "woodland"){
      source <- dplyr::filter(source, HabBroad %in% c("Woodland, broadleaved", "Woodland, mixed"))
      } else if (habitat == "grassland"){
         source <- dplyr::filter(source, HabBroad %in% c("Grassland, semi-natural", "Maritime cliff and slope"))
         } else if (habitat == "heath"){
            source <- dplyr::filter(source, HabBroad %in% c("Heathland"))
            } else if (habitat == "pond"){

## Calculate shp_area and shp_index

source <- source %>%
   dplyr::mutate(shp_area = as.numeric(sf::st_area(.)),
          shp_length = as.numeric(lwgeom::st_perimeter(.))) %>%
   dplyr::mutate( # calculate shape index
      shp_index = (pi * ((shp_length / (2 * pi)) ^ 2)) / shp_area)

               source <- dplyr::filter(source, HabBroad %in% c("Water, fresh"), shp_area < 10000, shp_index <= 5)

               } else if (habitat == "mire"){
                  source <- dplyr::filter(source, HabBroad %in% c("Grassland, Marshy", "Mire", "Swamp", "Saltmarsh"))

                  } else if (habitat == "wetland"){
                     source <- dplyr::filter(source, HabBroad %in% c("Grassland, marshy"))
                     }



   # Check geometry as needs to be clean for rasterizing
   source <- ecoservR::checkgeometry(source, "POLYGON")

   if (nrow(source) == 0) {
      message(paste0("No ", habitat, " habitat in your study area."))
      return() # exit function without crashing
   }

   # Need to union adjacent patches to calculate area correctly, but unioning a lot of polygons is a serious bottleneck. Instead we rasterize to identify patches, and then polygonize them.

   source_r <- fasterize::fasterize(source, r)

   source_r <- terra::rast(source_r)
   patches <- terra::patches(source_r, direction = 8)
   rm(source_r)

   source <- sf::st_as_sf(
      stars::st_as_stars(patches),
      as_points = FALSE,
      merge = TRUE) %>%
      dplyr::mutate(area_ha = as.numeric(sf::st_area(.))/10000)


   # # Union to keep only large enough patches - BOTTLENECK: could this be done from raster instead?
   # source <- source %>%
   #    rmapshaper::ms_simplify(keep = 0.1) %>%
   #    sf::st_buffer(5) %>% sf::st_buffer(-5) %>% ## buffer out and in to merge close enough patches
   #    sf::st_union() %>%
   #    sf::st_as_sf() %>%
   #    ecoservR::checkgeometry(.)

   # Calculate area and remove small bits

   source <- source %>%
      dplyr::filter(area_ha > minsource)

   if (nrow(source) == 0) {
      message(paste0("No ", habitat, " habitat patch greater than ", minsource, "ha in your study area."))
      return() # exit function without crashing
   }


   # Create source raster
   source_r <- fasterize::fasterize(source, r)

   message("Creating cost raster for ", habitat)
   # Create resistance raster
   cost_r <- fasterize::fasterize(x, r, field = costcol)


   if (use_hedges){

      # Create raster of hedge scores
      cost_hedge <- fasterize::fasterize(hedges, r, field = costcol)

      # Overwrite the basemap scores in places with hedges
      cost_r  <- raster::mask(cost_r, cost_hedge, inverse = T) %>%
         raster::cover(cost_hedge)

   }


   # Convert costs into transition matrix

   #1/mean: reciprocal of cost (resistance) to get permeability (conductance)
   #https://www.rdocumentation.org/packages/gdistance/versions/1.1-3/topics/accCost

   tr <- gdistance::transition(cost_r,
                               transitionFunction=function(x) {1/mean(x)},
                               directions=8)

   # Correct for diagonal movements
   tr <- gdistance::geoCorrection(tr, type="c", multpl=FALSE)


   ## The accCost function requires sets of coordinates as input.
   ## Let's transform each source cell into a Spatial Point

   source_pts <- raster::rasterToPoints(source_r, spatial = TRUE)


   ## Create the accumulated cost distance
   message("Calculating cost distance for ", habitat)
   costdistance <- gdistance::accCost(tr, source_pts)


   ## Reclassify raster to identify core (1), meso (2) and macro (3) networks based on dispersal distances

   costdistance <- raster::reclassify(costdistance,
                                      c(0, dispersal[[1]], 1,  # core network
                                        dispersal[[1]], dispersal[[2]], 2,  # meso network
                                        dispersal[[2]], dispersal[[3]], 3,  # macro network
                                        dispersal[[3]], Inf, NA),
                                      include.lowest = TRUE)

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

   # Polygonize the networks

   message("Creating network polygons...")

   core <- sf::st_as_sf(
      stars::st_as_stars(
         raster::reclassify(costdistance,
                            c(0,1,1,
                              1, Inf, NA))
      ),
      as_points = FALSE,
      merge = TRUE) %>%
      dplyr::mutate(area_ha = as.numeric(sf::st_area(.))/10000,
                    type = "core") %>% # calculate area of networks
      dplyr::filter(area_ha > threshold) # remove those smaller than functional threshold

   meso <- sf::st_as_sf(
      stars::st_as_stars(
         raster::reclassify(costdistance,  # combine core and meso
                            c(0,2,1,
                              2, Inf, NA))
      ),
      as_points = FALSE,
      merge = TRUE) %>%
      dplyr::mutate(area_ha = as.numeric(sf::st_area(.))/10000,
                    type = "meso") %>% # calculate area of networks
      dplyr::filter(area_ha > threshold) # remove those smaller than functional threshold

   macro <- sf::st_as_sf(
      stars::st_as_stars(
         raster::reclassify(costdistance, # combine core, meso and macro
                            c(0,3,1,
                              3, Inf, NA))
      ),
      as_points = FALSE,
      merge = TRUE) %>%
      dplyr::mutate(area_ha = as.numeric(sf::st_area(.))/10000,
                    type = "macro") %>% # calculate area of networks
      dplyr::filter(area_ha > threshold) # remove those smaller than functional threshold


   network_poly <- rbind(core, macro, meso)
   rm(core,macro,meso)

   if (nrow(network_poly) > 0){

      sf::st_write(network_poly,
                   dsn = file.path(save, paste0(projectLog$title, "_",
                                                runtitle, "_",
                                                "network_", habitat, "_poly.gpkg"
                   )),
                   append = FALSE)

      message("Saved ", habitat, " network polygons")


      ## Mask the raster and save

      if (export_raster){

      network_poly <- rmapshaper::ms_dissolve(network_poly) %>% sf::st_as_sf()

      # Save network raster
      message("Saving your ", habitat, " ecological network.")
      raster::writeRaster(
         raster::mask(costdistance, network_poly),
         filename = file.path(save,
                              paste0(projectLog$title, "_", runtitle, "_network_", habitat, ".tif")),
         overwrite = TRUE  # perhaps not desirable but for now prevents error messages
      )

      }


   } else {message("No networks larger than functional threshold (",threshold, " ha) in your study area")}


}