11  Mosaic Zone Timber Harvest Plan

In this target the baseline timber harvest plan is generated for the planted (Figure 11.1) and unplanted (Figure 11.2) mosaic strata (see Chapter 9). A harvest plan was developed in consultation with forestry exports in Sabah to replicate a realistic approach that might be taken if the Kuamut concession were to be logged. The assumption was made that logging would begin in the Northwest, where access is better. The extent of logging in each year is constrained by legal and practical constraints. In successive years, logging coupes are selected in a south easterly direction on the basis that, a logging campaign would move sequentially through the concession.

For mosaic strata, legislation states that trees with a DBH > 31 cm can be extracted. Therefore the 31 cm Merchantable volume map for 2016 was used to derive the timber volume estimates for each of the years. Ongoing forest growth is considered at a later stage (Chapter 21) therefore only the map for 2016 is required.

Table 11.1: Table of extracted timber volumes from Mosaic strata types for each year. Reported values for Extracted timber volume include a 6% reduction to account for protected tree species, see Chapter 10. Area values are reported following a 1.67% reduction to align with the legal definition of the site area.

The following maps depict the spatial sequencing of the baseline timber harvest plan for the planted and unplanted mosaic strata.

Figure 11.1: Planted mosaic timber harvest plan and corresponding timber volumes

Figure 11.2: Unlanted mosaic timber harvest plan and corresponding timber volumes

11.1 Code

#' Timber Harvest Plan for Mosaic
#'
#' @description This function creates the timber harvest plan for the Mosaic
#' Plantings coupes.
#'
#' @param boundary_src path to the boundary polygon
#' @param MosaicCoupesGo_src path to the Mosaic Coupes spatial polygon
#' @param MosaicPLANT_src path to the Mosaic Planting spatial polygon
#' @param MosaicUNplanted_src path to the Mosaic Unplanted spatial polygon
#' @param dbh31_src path to the 31cm DBH Merchantable Volume raster
#' @param dbh10_src path to the 10cm DBH Merchantable Volume raster
#' @param acd_src path to the Above Ground Carbon (ACD) raster
#' @param setup logging plot setup size in hectares (31ha was average in
#' INFAPRO area see Philipon 2020)
#' @param sp_conserve_factor derived from `R/protected-trees-stats.R`
#' @param LogAfterHowManyYears default is 15.
#'
#' @return
#'
#' @details
#' This function creates the timber harvest plan for the Mosaic Plantings
#' coupes.
timber_harvest_plan_mosaic <- function(
    boundary_src, MosaicCoupesGo_src,
    MosaicPLANT_src, MosaicUNplanted_src,
    dbh31_src, dbh10_src, acd_src,
    setup = 31, sp_conserve_factor,
    LogAfterHowManyYears = 15, Yrs2Plant = 30,
    Rotation = 15) {
  #----------- Read Data ------------------------
  ha <- (100 * 100) # 1 hectare is 100m x 100m
  setup_SqM <- setup * ha # gives size in square meters
  ploty <- sqrt(setup_SqM)
  plotx <- ploty


  # Project Boundary
  Boundary <- sf::read_sf(boundary_src)

  ## Outline of Area for Mosaic concessions
  Mosaic <- sf::read_sf(MosaicCoupesGo_src)

  # strata that will be planted
  MosaicPLANT <- sf::st_make_valid(sf::read_sf(MosaicPLANT_src))

  # strata that will NOT be planted
  MosaicUNplanted <- sf::st_make_valid(sf::read_sf(MosaicUNplanted_src))

  # Areas
  UNpU <- sf::st_union(MosaicUNplanted)
  UNp_ha <- units::set_units(sf::st_area(UNpU), value = ha)


  PlantU <- sf::st_union(MosaicPLANT)
  P_ha <- units::set_units(sf::st_area(PlantU), value = ha)

  MC_U <- sf::st_union(Mosaic)
  MC_ha <- units::set_units(sf::st_area(MC_U), value = ha)

  ## area of each concession
  # area one, five coupes
  Area1coupes <- c("YL1/01", "YL2/00", "YL2/01", "YL2/04", "YL3/02")
  Area1 <- dplyr::filter(Mosaic, BLOCK_N %in% Area1coupes)
  a1_ha <- units::set_units(x = sf::st_area(sf::st_union(Area1)), value = ha)

  # area two, 12 coupes
  Area2coupes <- c(
    "YS9/06", "YL3/03", "YSH7/06", "YSH2/06", "YSH8/06(Area A) Partial",
    "YLH4/05", "YLH4/06", "Steep area YLH406", "YLH3/06(2)Area B", "YLH3/06(1)",
    "YLH2/05", "YLH1/05"
  )
  Area2 <- dplyr::filter(Mosaic, BLOCK_N %in% Area2coupes)
  a2_ha <- units::set_units(sf::st_area(sf::st_union(Area2)), value = ha)


  # ----------- Convert Volume/hectare rasters to Absolute Volume -------------
  #          Run three times for 31cm & 10cm Merch Vol &ACD  (as a check)
  #----------------------------------------------------------------------------#

  fileName_list <- list(
    "31cmDBH" = dbh31_src,
    "10cmDBH" = dbh10_src,
    "ACD_10cm" = acd_src
  )

  # Convert Rasters from volume/ha to absolute volume per cell.
  # function in : 'R/coupes_helpers.R'
  abs_vol_list <- fileName_list |>
    purrr::imap(~ create_abs_volume(.x, .y, Boundary, Mosaic[1]))

  # planted parcels
  plant_polygon_out_dir <- dirname(MosaicPLANT_src)
  plant_polygon_out_shp <- file.path(
    plant_polygon_out_dir, "MosaicPLANT_THP_Grid.shp"
  )

  Plant_MosGrid <- create_parcel_grid(
    MosaicPLANT, "PLANTED",
    plant_polygon_out_shp, plotx, ploty
  )

  # unplanted parcels
  unplant_polygon_out_dir <- dirname(MosaicUNplanted_src)
  unplant_polygon_out_shp <- file.path(
    unplant_polygon_out_dir, "MosaicUNplanted_THP_Grid.shp"
  )
  NoPlant_MosGrid <- create_parcel_grid(
    MosaicUNplanted, "UNPLANTED",
    unplant_polygon_out_shp, plotx, ploty
  )

  # extract values for strata
  # function in : 'R/coupes_helpers.R'
  mosaic_strat_ext <- purrr::map(
    .x = abs_vol_list, ~ extract_strata(.x, Plant_MosGrid)
  )
  mosaic_unplant_ext <- purrr::map(
    .x = abs_vol_list, ~ extract_strata(.x, NoPlant_MosGrid)
  )

  #------------- Define logging rate and select coupes to match.---------------
  LogRate <- (0.15 + 0.11 + 0.25) / 3
  HaPerYr <- round(MC_ha * LogRate)

  # make plots for each logging year
  plant_ggplots <- purrr::imap(
    .x = mosaic_strat_ext,
    ~ plot_setups(.x, .y, "PLANTED", outline_width = 0.6)
  )
  unplant_ggplots <- purrr::imap(
    .x = mosaic_unplant_ext,
    ~ plot_setups(.x, .y, "UNPLANTED", outline_width = 0.2)
  )



  planted_future_parcel <- purrr::imap(
    .x = mosaic_strat_ext,
    ~ summarise_parcels_mos(.x, .y,
      sp_conserve_factor = sp_conserve_factor,
      Yrs2Plant = Yrs2Plant,
      Rotation = Rotation,
      P_ha = P_ha
    )
  )

  unplanted_future_parcel <- purrr::imap(
    .x = mosaic_unplant_ext,
    ~ summarise_parcels_mos(.x, .y,
      Coupe_planted = FALSE,
      sp_conserve_factor = sp_conserve_factor,
      Yrs2Plant = Yrs2Plant,
      Rotation = Rotation,
      P_ha = P_ha
    )
  )

  ParcelsIntoFutureBoth <- planted_future_parcel |>
    purrr::map2(
      .y = unplanted_future_parcel,
      ~ bind_rows(.x, .y)
    )

  # save parcel csvs.
  # function in : 'R/coupes_helpers.R'
  parcel_csv_list <- purrr::imap(
    ParcelsIntoFutureBoth,
    ~ save_parcels(.x, .y, "Mosaic_THP_Parcels")
  )

  return(
    c(
      parcel_csv_list,
      list(planted_mosaic_src = plant_polygon_out_shp),
      list(unplanted_mosaic_src = unplant_polygon_out_shp),
      Species_Conserv_Fac = sp_conserve_factor,
      list(mos_plant_ggplots = plant_ggplots),
      list(mos_unplant_ggplots = unplant_ggplots),
      list(mos_plant_spv = mosaic_strat_ext$`31cmDBH`),
      list(mos_unplant_spv = mosaic_unplant_ext$`31cmDBH`)
    )
  )
}

create_parcel_grid <- function(strata, type_name, out_path, plotx, ploty) {
  MosaicStrat <- strata

  # Make the grid of individual setups to be logged
  grid <- sf::st_make_grid(MosaicStrat, cellsize = c(plotx, ploty))

  gridmask <- sf::st_intersection(grid, MosaicStrat)

  # join attribute tables
  Mosaicgrid <- sf::st_join(
    sf::st_as_sf(gridmask),
    sf::st_as_sf(MosaicStrat),
    join = st_intersects, largest = TRUE
  )

  Mosaicgrid_polys <- sf::st_collection_extract(
    x = Mosaicgrid,
    type = "POLYGON"
  )
  Mosaicgrid <- Mosaicgrid_polys |>
    mutate(id = n()) |>
    mutate(
      LogYear = case_when(
        BLOCK_N %in% c("YL1/01", "YL2/01", "YL2/04") ~ 2016,
        BLOCK_N %in% c("YL3/02") ~ 2017,
        BLOCK_N %in% c("YL2/00", "YS9/06", "YSH7/06") ~ 2018,
        BLOCK_N %in% c("YL3/03", "YSH8/06(Area A) Partial") ~ 2019,
        BLOCK_N %in% c("YLH4/06", "YSH2/06", "YLH4/05", "Steep area YLH406")
        ~ 2020,
        BLOCK_N %in% c("YLH1/05", "YLH2/05", "YLH3/06(1)", "YLH3/06(2)Area B")
        ~ 2021,
        TRUE ~ -9999
      ),
      area_ha = units::set_units(sf::st_area(Mosaicgrid_polys), value = hectare)
    )

  Mosaicgrid <- Mosaicgrid |>
    select(!Area_ha) # dropping old area column.

  if (nrow(Mosaicgrid |> dplyr::filter(LogYear == -9999)) > 0) {
    stop("Check the Block Names!")
  }

  ## Check areas match. Argument join = st_nearest_feature important here
  region_list <- list(gridmask, MosaicStrat, Mosaicgrid)

  areas <- purrr::map_dbl(region_list, ~ get_areas(.x))
  # test that areas are equal
  purrr::walk(areas, ~ testthat::expect_equal(.x, areas[1]))

  sf::write_sf(Mosaicgrid, out_path, delete_dsn = TRUE)

  return(Mosaicgrid)
}

#' convert raster density to absolute values.
#'
#' @param ras_src raster source relative to data folder
#' @param ras_name the desired name of the raster
#' @param aoi the area of interest (spatial polygon or sf)
#' @param overlay coupes areas for plot vis
#' @param scale.factor default is 10000 i.e number of square meters in a ha.
#'
#' @return raster with absolute values.
#' @details
#' This function converts the raster density values to absolute values.
#' The raster density values are per hectare, so we need to convert them to
#' absolute values. The scale factor is the number of square meters in a
#' hectare.
#'
create_abs_volume <- function(
    ras_src, ras_name, aoi, overlay, scale.factor = 1e4) {
  V_BSL <- terra::rast(ras_src)[[1]]
  V_BSL * ((terra::res(V_BSL)[1] * terra::res(V_BSL)[2]) / scale.factor)
}


#' Extract sum
#'
#' Simple function to extract sum of raster within coupes.
#'
#' @param V_BSLpPix a raster object. Either Timber Volume or ACD.
#' @param strata_grid a gridded spatial vector of the strata of interest. sf object.
#'
#' @return sf object with sum of V_BSL per setup.
#' @details
#' extract sum of V_BSL per setup using MosaicGrid ( will be == V_EX when back
#' to per ha) UNITS should be sum per pixel
extract_strata <- function(V_BSLpPix, strata_grid) {
  strata_grid |>
    mutate(VolSum = exactextractr::exact_extract(
      V_BSLpPix, strata_grid, "sum"
    )) # use {exactextractr} for speed.
}



#' plot ACD or Timber Volume
#'
#' @param setups
#' @param dbh
#' @param mosaic_type
#' @param outline_width
#'
#' @return ggplot object.
plot_setups <- function(setups, dbh, mosaic_type, outline_width = 0.03, .ncol = 2) {
  p <- setups |>
    ggplot() +
    aes(fill = VolSum) +
    geom_sf(
      data = sf::st_geometry(dplyr::summarise(setups)), colour = "grey70", fill = "transparent",
      inherit.aes = FALSE, lwd = 0.1
    ) +
    geom_sf(data = setups, colour = "grey30", fill = "transparent", lwd = outline_width) +
    geom_sf(colour = NA) +
    scale_fill_viridis_c() +
    theme_light() +
    facet_wrap(~LogYear, ncol = .ncol) +
    labs(title = dbh) +
    theme(
      axis.title.x = element_blank(),
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks.y = element_blank(),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank()
    )
}


#' save parcel csv
#'
#' @param parcel_df the parcel summary df to save.
#' @param dbh one of: c("ACD_10cm", "40cmDBH", "31cmDBH", "10cmDBH")
#' @param save_pref character prefix for the file path
#'
#' @return file path (character) of the saved csv.
save_parcels <- function(parcel_df, dbh, save_pref) {
  # Save Parcels using 30cm DBH, 10cm DBH (or ACD as a check)
  if (dbh != "ACD_10cm") {
    what <- paste0("V_EX_", dbh)
  } else {
    what <- dbh
  }
  # Save each of the three Parcel files to csv

  parcels_save_path <- set_path(paste0(
    "RoutputProducts/",
    save_pref,
    what, ".csv"
  ))
  write.table(parcel_df, parcels_save_path)

  return(parcels_save_path)
}


#' helper function for initial parcel defs
#'
#' @param setups
#'
#' @return sf object of summarised parcels.
parcel_setup <- function(setups, dbh) {
  Parc1 <- setups |>
    group_by(LogYear) |>
    dplyr::summarise(VolSum = sum(VolSum, na.rm = TRUE))
  Parc <- Parc1 |>
    mutate(Area_ha = units::set_units(sf::st_area(Parc1), ha)) |>
    sf::st_drop_geometry()

  # For Asner Map version ACD, for merchantable, name is V_EX
  if (dbh == "ACD_10cm") {
    Parc <- mutate(Parc, ACD = round(VolSum / as.numeric(Area_ha)))
  } else {
    Parc <- mutate(Parc, V_EX = round(VolSum / as.numeric(Area_ha))) # -> V_EX
  }
  return(Parc)
}

#' function to create mosaic parcels
#'
#' @param setups
#' @param dbh
#' @param Coupe_planted
#'
#' @return sf object of summarised parcels.
summarise_parcels_mos <- function(setups, dbh, Coupe_planted = TRUE,
                                  sp_conserve_factor, Yrs2Plant, Rotation,
                                  P_ha) {
  ### Here we add up the volume of all setups into the yearly parcels.
  Parcels <- parcel_setup(setups, dbh)

  # In the Mosaics we will assume a baseline of extracting only 70% of all
  # merchantable timber this is a more conservative way to account for species
  # that shouldn't be logged.

  if (dbh == "31cmDBH") {
    #  # Don't extract all merchantable volume based on sp_conserve_factor
    Parcels <- mutate(Parcels, V_EX = round(V_EX * sp_conserve_factor, 1))
  }
  strata_def <- function(plant_type) {
    paste("Mosaic", plant_type, sep = "-")
  }
  Parcels <- Parcels |>
    mutate(
      Strata = case_when(
        isTRUE(Coupe_planted) ~ strata_def("planted"),
        TRUE ~ strata_def("unplanted")
      ),
      Plant = case_when(
        isTRUE(Coupe_planted) ~ "planted",
        TRUE ~ "unplanted"
      ),
      Rotation = "First"
    )

  #------------ Future rounds of logging ------------------------------------
  # For transparency and simplicity We keep predictions for Future logging
  # periods as straightforward as possible as these will be revisited
  #-------------------------------------------------------------------------#
  ## Modelling overall growth rate based on planting rate and rotation
  plantingArea <- as.numeric(P_ha)
  AReaPlantPA <- round(plantingArea / Yrs2Plant)

  ParcelsNextRotation <- data.frame(
    LogYear = seq(from = (2017 + Rotation), length.out = Yrs2Plant),
    VolSum = rep(NA, length.out = Yrs2Plant),
    Area_ha = rep(AReaPlantPA, length.out = Yrs2Plant),
    V_EX = rep(95 * sp_conserve_factor, length.out = Yrs2Plant),
    Strata = paste("Mosaic-planted"),
    Plant = paste("planted"),
    Rotation = paste("Second")
  ) |>
    mutate(VolSum = Area_ha * V_EX)

  # colname for ACD carbon map check
  if (dbh == "ACD_10cm") {
    colnames(ParcelsNextRotation)[4] <- "ACD"
  }

  # If the area is planted, there will be a future rotation
  # (according to above), if not planted, then not during this baseline period
  if (isTRUE(Coupe_planted)) {
    return(rbind(Parcels, ParcelsNextRotation))
  } else {
    return(Parcels)
  }
}


#' Function to create Natural Forest Management Parcels.
#'
#' @param setups
#' @param dbh
#'
#' @return sf object of summarised parcels.
summarise_parcels_nfm <- function(setups, dbh, sp_conserve_factor, R2ReEntry) {
  Parcels <- parcel_setup(setups, dbh)

  # In # the NFM we will assume a baseline of extracting only 70% of all
  # merchantable timber. This is a more conservative way to account for
  # species that shouldn't be logged and would fit with 'selective logging'.

  if (dbh == "40cmDBH") {
    # Don't extract all merchantable volume based on sp_conserve_factor
    Parcels <- mutate(Parcels, V_EX = round(V_EX * sp_conserve_factor, 1))
  }

  Parcels$Strata <- "NFM"
  Parcels$Plant <- "unplanted"
  Parcels$Rotation <- "First"

  # Future rounds of logging:
  # For transparency and simplicity We keep predictions for Future logging
  # periods as straightforward as possible as these will be revisited
  ParcelsNextRentry <- Parcels |>
    mutate(
      LogYear = LogYear + R2ReEntry,
      Rotation = "Second"
    )

  # given we have a conservative baseline extraction rate, we estimate the same
  # volume could and would be extracted after a further 15 years
  rbind(Parcels, ParcelsNextRentry)
}