R/mod09_classify_habitats.R

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

#######################################
### Worflow functions               ###
### for EcoservR                    ###
### Sandra Angers-Blondin           ###
### 11-01-2021                      ###
#######################################


# Notes -------------------------------------------------------------------

# This module uses the information extracted into the basemap to assign a habitat type to each polygon.
# The module requires at least module 01 to have run (basemap prep), and optionally (and ideally) also considers
# CORINE, Priority Habitats and the Crop Map of England.

## Still to do:
# consolidate original Ecoserv-GIS codes to remove duplication and "u" codes
# add more rules to consider extra datasets (NFI, ...) and sort out discrepancies between datasets (crome vs corine)
# add rules for Scottish datasets once modules written for them
# create a UKHab version




#' Classify Habitats (Phase 1)
#'
#' This function assigns a habitat type (modified from Phase 1 codes) to each polygon in the basemap according to the information available. Requires at least the basemap preparation function to have run, and optionally (and ideally) also considers CORINE, Priority Habitats, the Crop Map of England, hedgerows, and elevation data.

#' @param mm The mm object loaded in the environment, can be at various stages of updating.
#' @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_classified.RDS file to project folder
#' @export

classify_map <- function(mm = parent.frame()$mm,
                    projectLog = parent.frame()$projectLog){

   timeA <- Sys.time() # start time


   ## Extract parameters from project log -----

   params <- projectLog$parameters

   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)


   ## Check that mm is right format

      if (!inherits(mm, "list")){
      mm <- list(mm)
   }


   # Preparation -------------------------------------------------------------

   mm <- checkgeometry(mm, "POLYGON")

   # Make sure we have single-part polygons and calculate area afresh

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


   ### maybe we could work on a non-spatial object (indexed) and merge back into basemap after?

   attributes <- unique(unlist(lapply(mm, function(x) names(x))))


   ## Erase previous classification if necessary

   if ("HabCode_B" %in% attributes){

      mm <- lapply(mm, function(x) dplyr::select(x, -HabCode_B))
   }



# First round of classification -------------------------------------------

# First we need to make sure there is a GI column (sometimes there is no coverage of OS Greenspace / Open Greenspace for small or rural sites) as it is needed in the classif_mastermap function.

   if (!"GI" %in% attributes){
      mm <- lapply(mm, function(x) create_GI(x, params))
   }


# This is the main classification step and based on MasterMap, GI, and parameters only.

   mm <- lapply(mm, function(x) classif_mastermap(x, params))


# Refine with more GI rules

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


# Additional classification -----------------------------------------------

## Depending on datasets, additional classification

# Crome and Corine

   if (any(c("crome", "corine") %in% attributes)){

   mm <- lapply(mm, function(x) classif_agri(x))
   }

# Priority Habitats

if ("phi" %in% attributes){

   mm <- lapply(mm, function(x) classif_phi(x))
}

# National Forest Inventory
   if ("nfi" %in% attributes){

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

   }




# Further decide on agricultural land based on size
   mm <- lapply(mm, function(x) classif_area(x, params))


   # further decisions based on elevation and slope
   if ("elev" %in% attributes){
      mm <- lapply(mm, function(x) classif_elev(x, params))
   }

   ## Get rid of the B4/Bu intermediary code as there is no match for it in models
   ## TO DO add decisions to split, but now all to B4

   mm <- lapply(mm, function(x) dplyr::mutate(x,
                                       HabCode_B = dplyr::case_when(
      HabCode_B == "B4/Bu" ~ "B4",
      TRUE ~ HabCode_B))
   )




# GI accessibility classification -----------------------------------------

# Classify GI accessibility with a set of rules
   mm <- lapply(mm, function(x) classif_access(x))



# Finish and save ---------------------------------------------------------

# Save to disk

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


   ## Create a summary to report

   summary <- unlist(lapply(mm, function(x) x %>% sf::st_drop_geometry() %>%
                               dplyr::select(HabCode_B)
                            ))


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

   projectLog$last_success <- "MM_classified.RDS"

   timeB <- Sys.time() # stop time

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

   updateProjectLog(projectLog) # save revised log

   # and delete contents of scratch folder
   cleanUp(scratch_path)


   message(paste0("Finished classifying map into ", length(unique(summary)), " habitat types. \n
                  Process took ",
                  round(difftime(timeB, timeA, units = "mins"), digits = 1),
                  " minutes."))



   # Return mm to environment
   return({
      invisible({
         mm <<- mm
         projectLog <<- projectLog
      })
   })

} # end of function