R/mod05_add_PHI.R

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

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

#' Add Priority Habitat Inventory
#'
#' This function adds the Priority Habitat Inventory data to the basemap (England).

#' @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_05.RDS file to project folder
#' @export

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

   timeA <- Sys.time() # start time

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

   # Get file path
   phipath <- projectLog$df[projectLog$df$dataset == "phi", ][["path"]]  # path to corine data, if available
   dsname <- projectLog$df[projectLog$df$dataset == "phi", ][["prettynames"]] # dataset name


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

      message("Preparing to update baseline with Priority Habitat data")

      phitype <- guessFiletype(phipath)   # file type, gpkg or shp

      phi_cols <- tolower(  # making all lowercase for easier matching
         projectLog$df[projectLog$df$dataset == "phi", ][["cols"]][[1]]  # attributes
      )


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

# Import PHI data

   phi <- loadSpatial(phipath, filetype = phitype)  # using the loading into list function in case there are many shp

   phi <- do.call(rbind, phi) %>% sf::st_as_sf()  # putting back into one single sf object


# DATA PREP -----------------------------------------------------------------------------------

   # Rename columns

   names(phi) <- tolower(names(phi)) # forcing lowercase attributes

   phi <- dplyr::select(phi, tidyselect::all_of(phi_cols)) %>%  # remove and rename columns
      dplyr::select(bapHab = Main_Habit)

   phi <- checkcrs(phi, studyAreaBuffer)  # check and transform crs

   phi <- faster_intersect(phi, studyAreaBuffer)

   # Tidy up the data
   # now keeping woodlands as they can confirm reclassification

   phi <- phi %>%
         sf::st_make_valid() %>%                     # check and repair geometry
         sf::st_cast(to = "MULTIPOLYGON") %>%        # multi to single part
         sf::st_cast(to = "POLYGON", warn = FALSE)


# Rasterize data ------------------------------------------------------------------------------

      ## These custom functions will check data size and create tiles if required, otherwise will process everything in memory

      phi_v <- prepTiles(mm, phi, studyArea = studyAreaBuffer, value = "bapHab")
      rm(phi)

      if(is.null(phi_v)){

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

         return(message("WARNING: Priority Habitat Inventory data not added: No data coverage for your study area."))
      }

      phi_r <- makeTiles(phi_v, value = "bapHab", name = "PHI")
      rm(phi_v)



# ZONAL STATISTICS TO UPDATE MASTERMAP ----------------------------------------------------------

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


      mm <- mapply(function(x, n) extractRaster(x, phi_r, fun = "majority", tile = n, newcol = "phi"),
                   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(phi_r)  # remove raster tiles

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

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


# SAVE UPDATED MASTER MAP ---------------------------------------------------------------------

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

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

      projectLog$last_success <- "MM_05.RDS"

      timeB <- Sys.time() # stop time

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

      updateProjectLog(projectLog) # save revised log

      # remove objects
      rm(key)

      # and delete contents of scratch folder
      cleanUp(scratch_path)

      message(paste0("Finished updating with Priority Habitat Inventory data. Process took ",
                     round(difftime(timeB, timeA, units = "mins"), digits = 1),
                     " minutes."))

   } else {message("No Priority Habitat data input specified.")} # end of running condition


   # Return mm to environment, whether it has been updated or not.
   return({
      invisible({
         mm <<- mm
         projectLog <<- projectLog
      })
   })

} # end of function