ecoserv_tool: Natural Capital Mapping Toolkit (UK)

Documented in add_greenspace

#######################################
### Worflow functions               ###
### for EcoservR                    ###
### Sandra Angers-Blondin           ###
### 05 October 2020                 ###
#######################################

#' Add Greenspace
#'
#' This function adds OS Greenspace and/or OS Open Greenspace to the basemap.

#' @param mm The mm object loaded in the environment, can be at various stages of updating.
#' @param studyAreaBuffer The buffered study area generated during mod01 or reloaded when resuming a session.
#' @param projectLog The RDS project log file generated by the wizard app and containing all file paths to data inputs and model parameters
#' @return Saves a project_title_MM_02.RDS file to project folder
#' @export
#'
add_greenspace <- function(mm = parent.frame()$mm,
                           studyAreaBuffer = parent.frame()$studyAreaBuffer,
                           projectLog = parent.frame()$projectLog){
   # the parent frame bit makes sure the function knows that the argument comes from the users's environment

   timeA <- Sys.time() # start time

   # NOTES ---------------------------------------------------------------------------------------

   # This script first updates the map with OS Greenspace (master version) if available
   # Then it uses OS Open Greenspace on anything that has not been classified by the licensed version.


   ## Extract the file paths and other info from project log ----------------------

   output_temp <- projectLog$output_temp
   title <- projectLog$title
   scratch_path <- file.path(output_temp, "ecoservR_scratch")
   if (!dir.exists(scratch_path)) dir.create(scratch_path)

   # OS GREENSPACE
   greenpath <- projectLog$df[projectLog$df$dataset == "OS_Greenspace", ][["path"]]  # path to the OS MasterMap Greenspace data, if available



   # OS OPEN GREENSPACE
   opengreenpath <- projectLog$df[projectLog$df$dataset == "OS_OpenGreenspace", ][["path"]]  # path to OS Open Greenspace (REQUIRED)

   dsname <- projectLog$df[projectLog$df$dataset == "OS_OpenGreenspace", ][["prettynames"]]

   opengreenlayer <- projectLog$df[projectLog$df$dataset == "OS_OpenGreenspace", ][["layer"]] # layer name
   opengreentype <- guessFiletype(opengreenpath) # data type, automatically determined

   opengreen_cols <- tolower( # make col names lowercase for easier matching
      projectLog$df[projectLog$df$dataset == "OS_OpenGreenspace", ][["cols"]][[1]]
      ) # coerce to named character (remove list nesting)


   # DATA IMPORT ---------------------------------------------------------------------------------

   ## Import the greenspace master map, if present
   # Stored in tiles of 5km and all go in list, regardless of whether or not they are in sub-directories. Only runs if user specified a path.

   if (!is.na(greenpath) & !is.null(greenpath)){ # if greenpath is empty, skipping this step

      green_cols <- tolower(projectLog$df[projectLog$df$dataset == "OS_Greenspace", ][["cols"]][[1]]) # coerce to named character (remove list nesting)

      greentype <- guessFiletype(greenpath) # data type; either gpkg or shp, automatically determined

      green <- loadSpatial(folder = greenpath, filetype = greentype)   # looks for and imports any shapefile within this folder, initially in a list

      # force all names to lowercase for easier matching
      green <- lapply(green, function(x){
         names(x) <- tolower(names(x))
         return(x)
      })
   }


   # GREENSPACE CLEANING -------------------------------------------------------------------------

   # Doing these operations on list elements because faster than using distinct() on the full combined data object.
   # This long section will only run if user has loaded a MasterMap Greenspace dataset.

   if (exists("green")){
      message("Updating MasterMap with OS Greenspace")

      # Clean the green object to keep only what we need

      green <- lapply(green, function(x) dplyr::select(x, tidyselect::all_of((green_cols))) %>%  # keep only needed cols before merge
                         dplyr::mutate(TOID = gsub("[a-zA-Z ]", "", TOID))   # remove the "osgb" characters from the TOID column
      )


      ### Remove duplicated polygons (custom function) -----

      green <- removeDuplicPoly(green, ID = "TOID")


      ### Remove empty features -----
      # If the user downloaded tiles for the whole study area, some of them might be completely empty and cause further steps to crash because of NA coordinates, so they get removed here.

      green <- green[sapply(green, function(x) dim(x)[1]) > 0]

      ### Crop to the study area -----

      green <- checkcrs(green, studyAreaBuffer)  # check CRS and transform if required

      green <- suppressWarnings(  # forcing the crs for compatibility in clipping function
         lapply(green, function(x) sf::st_set_crs(x, 27700)))

      ## Make sure features are clean before cropping

      green <- checkgeometry(green, "POLYGON")


      ### Removing empty objects -----

      # Some list items can become empty after removing duplicates or cropping to the study area (not good coverage in rural zones)

      green <- green[!is.na(green)]      # remove NAs generated when tile was fully out of study area
      green <- green[sapply(green, function(x) dim(x)[1]) > 0]  # drops objects with empty geometries after clipping


      ### Bind all elements of the list together in one object -----

      # Drop geometry
      green <- lapply(green, function(x) sf::st_drop_geometry(x))

      green <- data.table::rbindlist(green)


      ### There could be no features within the study area. We only do the following if there are features to match, and otherwise return a message.

      if (nrow(green) > 0){


         # UPDATE MASTERMAP WITH GREENSPACE ------------------------------------------------------------

         ### Perform merge based on TOID -----

         # for each mm list element, merge (left join) with toid of greenspace, return list of updated objects

         mm <- lapply(mm, function(x) dplyr::left_join(x, green, by = c("TOID")))

         rm(green)

         message("Finished updating with OS Greenspace")
         rm(i)

      } else {
         message("No Greenspace features within your study area.")
         # add a priFunc column anyway (with NA values)
         mm <- lapply(mm, function(x) dplyr::mutate(x, priFunc = NA))
         }

   } else {
      mm <- lapply(mm, function(x) dplyr::mutate(x, priFunc = NA))  # if not using Greenspace, create empty columns so classification rules still apply
      message("No OS Greenspace data provided, skipping to Open Greenspace.")
   }


   if (!is.null(opengreenpath) & !is.na(opengreenpath)){


      # OPEN GREENSPACE CLEANING --------------------------------------------------------------------


      # Import the OPEN greenspace map (required)

      opengreen <- loadSpatial(folder = opengreenpath,
                               layer = opengreenlayer,
                               filetype = opengreentype)
      opengreen <- do.call(rbind, opengreen) %>% sf::st_as_sf()  # remove from list and store into one sf object


      message("Updating MasterMap with OS Open Greenspace")

      # Rename columns
      names(opengreen) <- tolower(names(opengreen)) # force lowercase for easier matching

      ### Crop layer to study area and tidy up ----

      # Check projection: make sure same as studyAreaBuffer so we can clip

      opengreen <- checkcrs(opengreen, studyAreaBuffer)

      opengreen <- opengreen %>%
         dplyr::select(all_of(opengreen_cols)) %>%           # keep only needed cols
         faster_intersect(., studyAreaBuffer)


      opengreen <- checkgeometry(opengreen, "POLYGON")  # check validity and cast to polygon
      # %>%     # crop to study area
      #    sf::st_make_valid() %>%                          # check and repair geometry
      #    sf::st_cast(to = "MULTIPOLYGON") %>%
      #    sf::st_cast(to = "POLYGON", warn = FALSE)        # equivalent of multi-part to single part


      # RASTERIZE OPEN GREENSPACE -------------------------------------------------------------------

      ## Rasterize OPEN greenspace

      # The prepTiles function (custom function) will check if the data can be handled in memory.
      # If so, it will not tile the dataset but just put it in a list. If the dataset is too large, it will get cropped to 10x10km tiles that can then be rasterized in turn.

      opengreen <- prepTiles(mm, opengreen, studyArea = studyAreaBuffer, value = "op_function")

      ## Rasterize and extract into basemap, unless there is no coverage at all in study area (in which case we exit the function early with a message)

      if(is.null(opengreen)){

         projectLog$ignored <- c(projectLog$ignored, dsname)
         updateProjectLog(projectLog)


         return(message("WARNING: OS Open Greenspace not added: No data coverage for your study area."))
      }


      # We then feed this list object of vector tiles into the makeTiles function (another custom function).
      # This will rasterize the vector data into a list of tiles stored temporarily on disk: this prevents running into memory issues for large datasets.

      opgreen_r <- makeTiles(vectlist = opengreen, value = "op_function",
                             scratch = scratch_path,
                             name = "opgreen")  # ref name for the temp files

      # Create a key from the rasters' levels to add the descriptions back into the mastermap after extraction
      key <- as.data.frame(raster::levels(opgreen_r[[1]]))


      rm(opengreen)  # remove vector version



      # ZONAL STATISTICS OPEN GREENSPACE ------------------------------------------------------------

      # Extract raster values into basemap. If working with named tiles, the function matches tiles 1 to 1 and the process is much faster.

      mm <- mapply(function(x, n) extractRaster(x, opgreen_r,
                                                fun = "majority",
                                                tile = n, newcol = "GI_function"),
                   x = mm,
                   n = names(mm), # passing the names of the tiles will allow to select corresponding raster, making function faster. If user is not working with named tiles, will be read as null and the old function will kick in (slower but works)
                   SIMPLIFY = FALSE)  # absolutely necessary

      rm(opgreen_r)  # remove raster tiles

      # Replace the greenspace numeric codes by the text description (also a custom function)

      mm <- lapply(mm, function(x) addAttributes(x, "GI_function", key))

   } else {
      # if not using OS Greenspace, add the GI_function attribute with NAs
      mm <- lapply(mm, function(x) dplyr::mutate(x, GI_function = NA))
   }

   # CREATE COMBINED GREENSPACE CODE -------------------------------------------------------------

   # This section uses information from both greenspace datasets to create a unified greenspace description

   # create a function with all the rules
   classGI <- function(x) {
      dplyr::mutate(x,
                    GI = dplyr::case_when(
                       !is.na(priFunc) ~ as.character(priFunc),    # if there's a code from master green, use it
                       is.na(priFunc) & is.na(GI_function) ~ "Not Greenspace",   # all NAs mean no greenspace
                       is.na(priFunc) & !is.na(GI_function) ~ as.character(GI_function) # if not filled from master and open is present, use that
                    )) %>%

         # when buildings (and paths?) have been classified as greenspace, correct this
         dplyr::mutate(GI = ifelse(
            GI != "Not Greenspace" & Make == "Manmade", "Not Greenspace", GI
         )) %>%

         # when classified non greenspace but Mastermap says natural, create new category
         dplyr::mutate(GI = ifelse(
            GI == "Not Greenspace" & Make == "Natural", "Undetermined Greenspace", GI
         )) %>%

         # remove unnecessary columns
         dplyr::select(-priFunc, - GI_function) %>%

         # lump amenity greenspaces into 1 category
         dplyr::mutate(
            GI = ifelse(grepl("menity",GI), "Amenity", GI)
         )

   }


   mm <- lapply(mm, function(x) classGI(x))


   ## To think about: paths in greenspace area?



   # SAVE UPDATED MASTERMAP ----------------------------------------------------------------------


   saveRDS(mm, file.path(output_temp,  paste0(title,"_MM_02.RDS")))

   rm(key)

   # Update the project log with the information that map was updated

   projectLog$last_success <- "MM_02.RDS"

   timeB <- Sys.time() # stop time

   # add performance to log
   projectLog$performance[["add_greenspace"]] <- as.numeric(difftime(
      timeB, timeA, units="mins"
   ))

   updateProjectLog(projectLog) # save revised log


   # and delete scratch folder
   cleanUp(scratch_path)



   message(paste0("MasterMap updated with Greenspace and saved. Process took ",
                  round(difftime(timeB, timeA, units = "mins"), digits = 1),
                  " minutes."))

   return({
      invisible({
         mm <<- mm
         projectLog <<- projectLog
      })

   })


}