R/fun_export.R

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

#######################################
### Worflow functions               ###
### for EcoservR                    ###
### Sandra Angers-Blondin           ###
### 08 April 2021                   ###
#######################################

#' Finalise Basemap
#'
#' This function tidies up the final basemap and exports it to the "final" folder. It also outputs summary statistics. (NB: these should be checked, and recalculated if necessary, as MasterMap data may contain some overlapping polygons which will affect the area calculations.)

#' @param mm The mm object loaded in the environment. It should have all the desired data added and classified.
#' @param projectLog The RDS project log file generated by the wizard app and containing all file paths to data inputs and model parameters
#' @param export Should a clipped, recombined version of the map be saved as a GIS file? Default TRUE
#' @param format The desired output format, only used if export = TRUE. Currently supported are "gpkg" (geopackage, the default), "ESRI Shapefile" and "MapInfo File".
#' @param tiles Should the basemap be exported as 10x10km tiles? Defaults to FALSE but may be required for large study areas which may take a long time to render in a GIS software. Only used when export = TRUE.
#' @return Saves a project_title_final.RDS and optionally a project_title_final GIS file to project folder.
#' @export
#'
finalise_basemap <- function(mm = parent.frame()$mm,
                             projectLog = parent.frame()$projectLog,
                             export = TRUE,
                             format = "gpkg",
                             tiles = FALSE){
   # the parent frame bit makes sure the function knows that the argument comes from the users's environment

   timeA <- Sys.time() # start time

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

   studypath <- projectLog$df[projectLog$df$dataset == "studyArea", ][["path"]]

   final_folder <- file.path(projectLog$projpath, "final")
   title <- projectLog$title


   ### Create output folder if doesn't exist
   if (!dir.exists(final_folder)){
      dir.create(final_folder)
   }



# Tidy up geometries and save RDS -----------------------------------------

## This will be used in ES models so should still be the buffered version

mm <- checkgeometry(mm, "POLYGON")
mm <- checkcrs(mm, 27700)

saveRDS(mm, file.path(final_folder, paste0(title,"_basemap_final.RDS")))

message("RDS version saved as ", file.path(final_folder, paste0(title,"_basemap_final.RDS")),
        ". This is the version to use in ES models.")


if (export){

# Clip final map to study area --------------------------------------------

## For the final presentation version and for calculating statistics,
## the map needs to be clipped to the actual study area


   ### Reimport the study area outline (specifying OSGB as crs)

   studyArea <- try(loadSpatial(studypath,
                                layer = NULL,
                                filetype = guessFiletype(studypath)) %>%
                       do.call(rbind, .))

   if (inherits(studyArea, "try-error")) stop("Could not import study area. Check file path and format.") else message("Clipping to study area")

   studyArea <- suppressWarnings({
      checkcrs(studyArea, 27700) %>%   # check that CRS is Brit National Grid, and transform if not
         sf::st_set_crs(., 27700) %>%  # set the crs manually to get rid of GDAL errors with init string format
         sf::st_zm(drop = TRUE) %>%
         sf::st_make_valid() %>% sf::st_geometry() %>% sf::st_as_sf()  # retain only geometry
   })


   ## identify core vs edge tiles
   SAgrid <- sf::st_set_crs(ecoservR:::grid,27700)[
      lengths(sf::st_intersects(sf::st_set_crs(ecoservR:::grid,27700), studyArea))>0,]
   is_core <- unlist(sf::st_contains(studyArea,SAgrid))

   ## Clip where necessary
   # (if a whole 10k tile is included within study area, no need)

   for (i in 1:length(mm)){


      # we only need to perform clipping on edge tiles
      if (names(mm)[[i]] %in% names(SAgrid)[is_core]){next}

      ## if we have an edge tile, we apply the faster_intersect custom function
      ## which deletes polygons outside the boundary, protects those inside
      ## and clips those on the edge

      mm[[i]] <- faster_intersect(mm[[i]], studyArea)
      message("Clipped tile ", names(mm)[[i]])

   }

   message("Repairing geometries")

   mm <- checkgeometry(mm, target = "POLYGON")


   message("Calculating area")

   mm <- lapply(mm, function(x) dplyr::mutate(x, area = as.numeric(sf::st_area(x))))


   ## Add additional info from lookup table

   info <- ecoservR::hab_lookup %>% dplyr::select(Ph1code, HabNmPLUS, HabBroad, HabClass)

   mm <- lapply(mm, function(x) dplyr::left_join(x, info, by = c("HabCode_B" = "Ph1code")))



   ## Recombine map in one big object

   ## if too many tiles, issue warning and suggest using tiled outputs instead

   if (tiles == FALSE){

      ## Check how many tiles
      ntiles <- length(mm)

      if (ntiles >= 20){message("Warning: Large study area. You might want to use tiles = TRUE to output map tiles instead.")}

      mm <- do.call(rbind, mm)

      if (format == "gpkg"){
         sf::st_write(mm,
                  dsn = paste0(final_folder, "/", title, "_basemap_final.gpkg"),
                  append = FALSE)
      } else {
         sf::st_write(mm,
                  dsn = final_folder,
                  layer = paste0(title, "_basemap_final"),
                  driver = format,
                  append = FALSE)
      }

   } else {

      # save series of named tiles
      if (format == "gpkg"){
         mapply(function(x, n){
            st_write(x,
                     dsn = file.path(final_folder, paste0(title, "_", n,"_basemap_final.gpkg")),
                     append = FALSE)
         },
         x = mm,
         n = names(mm))

      } else {
         mapply(function(x, n){
            st_write(x,
                     dsn = final_folder,
                     layer = paste0(title, "_", n,"_basemap_final"),
                     driver = format,
                     append = FALSE)
         },
         x = mm,
         n = names(mm))

      }
   }


   message("Saved final map. Calculating summary statistics...")

   # Statistics ----------------------------------------------------------------------------

if (inherits(mm, "list")){
   summary <- data.table::rbindlist(lapply(mm, function(x) sf::st_drop_geometry(x)))
} else {
   summary <- sf::st_drop_geometry(mm)
}

   tot_area_ha <- sum(unlist(as.numeric(sf::st_area(studyArea))))/10000

   summary <- summary %>% dplyr::group_by(HabBroad) %>%
      dplyr::summarise(area_ha = sum(area)/10000) %>%
      dplyr::ungroup() %>%
      dplyr::mutate(cover_pct = area_ha/tot_area_ha*100)

#summary$area_ha <- round(summary$area_ha, digits = 2)
#summary$cover_pct <- round(summary$cover_pct, digits = 2)

## Sort summary in decreasing order

   summary <- summary[order(summary$area_ha, decreasing = TRUE),]

# Save a HTML summary table
   kableExtra::save_kable(
      kableExtra::kbl(summary, format = "html",
                      digits = 2,
                      col.names = c("Habitat type", "Area (ha)", "Cover (%)"),
                      caption = paste0(toupper(title),
                                       ": Natural capital asset register (Automatically generated by EcoservR on ", Sys.Date(),")")
                      ) %>%
         kableExtra::kable_styling(),
      file=file.path(final_folder, paste0(title, "_natcap_stats.html")))

} # end of export

   timeB <- Sys.time() # stop time

   message(paste0("Export completed. Process took ",
                  round(difftime(timeB, timeA, units = "mins"), digits = 1), " minutes."))


   on.exit(invisible(gc())) # garbage collection - return some memory to computer

}



### Also need a sample_basemap function that exports either the full map or a sample of it (using radius in meters)


#' Preview Basemap
#'
#' This function exports the basemap to a spatial format suitable for previewing in a GIS software. It is useful to frequently check on the outputs during the basemapping process to make sure that the coverage is as expected and that the data extracted look reasonable. NB: This is NOT the way to export the final map: use finalise_basemap for that.

#' @param mm The mm object loaded in the environment. It should have all the desired data added and classified.
#' @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
#' @param sample A length (in meters) for the side of a sampling square centered on the geographic middle of the map. Useful to export a small chunk of the map for a rapid check rather than exporting the full object. We recommend 5000 m.
#' @return Saves a "project_title_map_preview.gpkg" to project folder.
#' @export
#'
preview_basemap <- function(mm = parent.frame()$mm,
                            studyAreaBuffer = parent.frame()$studyAreaBuffer,
                            projectLog = parent.frame()$projectLog,
                            sample = NULL){
   # the parent frame bit makes sure the function knows that the argument comes from the users's environment


   if (!is.null(sample)){
      if (!inherits(sample, "numeric")){stop("sample must be a numeric value (in m)")}

      # get centroid of study area
      centro <- suppressWarnings(sf::st_centroid(studyAreaBuffer))

      # create (circle) buffer around it
      centro <- sf::st_buffer(centro, sample/2)

      # create square from bounding box
      square <- sf::st_bbox(centro) %>% sf::st_as_sfc() %>% sf::st_make_valid()


      # identify which mm tiles intersect
      SAgrid <- ecoservR::grid_study(studyAreaBuffer)  # create gridded study area
      index <- names(SAgrid) # extract names

      SAgrid <- do.call(rbind, SAgrid)  # unlist the grid squares

      index <- index[lengths(sf::st_intersects(SAgrid, square)) >0]  # extract tile names when tile intersects the sampling area


      x <- mm[index]
      x <- do.call(rbind, x) %>% sf::st_make_valid()
      x <- suppressWarnings(sf::st_intersection(x, square))


   } else {

      x <- do.call(rbind, mm)
   }

   ## Export ----

   x <- checkgeometry(x, "POLYGON")

  sf::st_write(x,
           dsn = file.path(projectLog$projpath, paste0(projectLog$title,"_basemap_preview.gpkg")),
           append = FALSE)

   message("Preview exported to ", file.path(projectLog$projpath, paste0(projectLog$title,"_basemap_preview.gpkg")))


   on.exit(invisible(gc())) # garbage collection - return some memory to computer

}