8  Slope Analysis

This target provides an understanding of the regions of the project area which are broadly less steep. This is key in the spatial formulation of the baseline. It is expected that coupes which fall in those less steep regions would have adopted Mosaic Plantings management and steeper coupes would adopt Natural Forest Management. These “flat” regions therefore contribute to the decision making process for the assignment of management strata to the various logging coupe areas as presented in Chapter 9.

To simplify the description we refer to “flat” and “steep” areas but acknowledge that, given the steep topography at the project, neither region can be considered flat. A k-means self-supervised classification was applied to the slope raster, allowing two centres, resulting in the creation of two slope clusters. Clusters were then vectorised and holes in the spatial vectors were filled.

Figure 8.1 below depicts the “flat” regions in orange overlaying the slope raster used to derive the extent.

Figure 8.1: Kmeans analysis of slopes to identify regions of the site that are broadly flatter and steeper to support qualitative attribution fo strata areas.

8.1 Code

#' Analysis of slope to define flat and steep regions.
#'
#' @param dem_src a raster source or file path readable by {terra}. In this case
#' from `mergeDEM_Sabah.R`.
#' @param pixel_size numeric - resolution of the slope raster in meters.
#' @return file path for flat area polygon.
#' @details
#' Runs k-means cluster analysis to determine steep vs. flatter regions.
#' large flat areas are extracted and holes filled. used to inform selection of
#' coupes that would have experienced Mosaic plantings management under the
#' baseline scenario.
slope_analysis <- function(dem_src, pixel_size = 100) {
  # Load raster.
  demKuamut <- terra::rast(dem_src) # see mergeDEM_Sabah.R

  # use gdal warp to resample to 1ha grid
  demResamp <- ezwarp::ezwarp(
    demKuamut, demKuamut,
    res = pixel_size, engine = "sf"
  )

  # define slope raster path.
  slp_ras <- set_path(
    paste0("RoutputProducts/KuamutSlope", pixel_size, "m.tif")
  )

  # build slope raster
  slp_ha <- terra::terrain(demResamp,
    v = "slope", unit = "degrees", neighbors = 8,
    filename = slp_ras, overwrite = TRUE
  )

  HaBlock <- units::as_units(pixel_size^2, "m^2") |>
    units::set_units(ha) |>
    as.numeric()

  # ----------------------- K means cluster analysis ------------------------- #
  # cluster-analysis to contiguous 'flat areas' and 'steep areas'
  rastVals <- terra::values(slp_ha) # get values of slope raster
  valInd <- which(!is.na(rastVals)) # get index of non NA values.
  set.seed(42) # set seed for k means, 42 is just a nod to Douglas Adams.
  myclust <- kmeans(rastVals[valInd], centers = 2)

  clusterRst <- slp_ha # make copy of slope raster
  rastVals[valInd] <- myclust$cluster # assign values from kmeans to matrix
  terra::values(clusterRst) <- rastVals # assign matrix values to raster

  # convert raster to polygon
  clustPoly <- terra::as.polygons(clusterRst, dissolve = TRUE) |>
    sf::st_as_sf()

  clustPoly <- clustPoly |>
    mutate(AREA = sf::st_area(clustPoly))

  .flatID <- clustPoly$slope[which.max(clustPoly$AREA)]
  .steepID <- unique(clustPoly$slope)[!unique(clustPoly$slope) %in% .flatID]

  clustPoly <- clustPoly |>
    mutate(slope_class = case_when(
      slope == .flatID ~ "flat",
      slope == .steepID ~ "steep",
      TRUE ~ NA_character_
    )) |>
    select(!AREA)

  # ----------------------------- Clean Polygons ----------------------------- #
  Mosaic_flat <- clustPoly |>
    # multi single part
    wk::wk_flatten()
  Mosaic_nohole <- Mosaic_flat |>
    # remove areas < 7k ha
    dplyr::filter(sf::st_area(Mosaic_flat) > units::as_units(7000, "ha")) |>
    # remove holes in polyon geometry
    sfheaders::sf_remove_holes()
  Mosaic <- Mosaic_nohole |>
    mutate(
      areas = sf::st_area(Mosaic_nohole), # get areas for each region
      areas_ha = units::set_units(areas, ha)
    ) |>
    dplyr::filter(slope_class == "flat") # keep only the flat areas.

  # save areas as shp file.
  cluster_dir <- set_path("RoutputProducts/ClusterAnalysisFlatOutline")
  if (!dir.exists(cluster_dir)) dir.create(cluster_dir)

  flat_outline_path <- file.path(
    cluster_dir,
    paste0("TwoFlatAreas", HaBlock, "haRes", ".shp")
  )

  sf::write_sf(Mosaic, flat_outline_path, delete_dsn = TRUE)

  return(list(
    slope_raster = slp_ras,
    flat_areas = flat_outline_path
  ))
}