5 Slope Correction

When measuring the extent of forest plots, we adopted a sampling methodology according to Jucker et al. (2018). A 30 m radius was measured along the topographic surface when defining the outline of the plot on the ground, i.e the circular plot conformed to the land surface, not a planar surface. Given the availability of high resolution LiDAR across all plot locations, the mean slope was extracted for plot locations which was used to convert the plot area to planimetric area, according to Equation 5.1:

\[ A_{ha} = \frac{A_p \times \cos(S)}{1e^4} \tag{5.1}\] Where \(A_{ha}\) is Area in hectares, \(A_p\) is planemtric area in square meters and \(S\) is slope.

An accurate estimate of forest plot area is essential for calculating the plot-level above ground carbon densities as described in Chapter 6.

Figure 5.1 shows the plot locations with their associated planimetric area values. Also provided is a slope raster, note how plots on steeper slopes have a smaller planimetric area values.

Figure 5.1: Plot locations and their calculated planimetric areas, overlaying the slope raster used to derive these adjusted forest plot areas.

5.1 Code

Slope Correction Code R/SlopeCorrection.R

#' Slope Correction function for plot areas.
#'
#' This file calculates the Area (in hectares) of each Sample Plot (Area_SamplePlot)
#' with consideration of the plot's slope
#'
#'
#' @param dem_path
#' @param points_path
#' @param radius numeric - radius of plot size measured along the ground
#'
#' @details
#' Plots where set up by measuring the distance along the ground
#' Not by horizontal projection
#' Here we slope-correct the area of each plot
#' by multiplying by cos(<U+03B8> ), where <U+03B8> is the average slope of the
#' plot in degrees as calculated from the digital elevation model
#' DEM obtained from the ALS data.
#' See paper on slope correction
#' https://www.nrcs.usda.gov/Internet/FSE_DOCUMENTS/nrcs142p2_041910.pdf
#' Plot on steep slopes have a smaller Horizontal projection
#' Dummy check results with Table 1
#' See also Jucker (2018) - model for Sabah carbon map where the same approach was taken
#' Jucker (2018) Biogeosciences;15(12):3811<U+2013>30.
#' 6% slope difference
#' NB: if results contain NA this means that the buffer doesnt cover even
#' 1 single cell (example at intersection of 4 cells)
#'
#'
#' @return dataframe of locations including Sample plot Area
#' @export
#'
#' @examples
slope_correction <- function(dem_path,
                             points_path,
                             site = "kuamut",
                             radius = 30) {
  # plot area
  area <- pi * radius^2
  # read dtm
  dem_kuamut <- terra::rast(dem_path)
  # calculate slope

  # define slope raster path.
  slp_ras <- set_path(paste0("RoutputProducts/KuamutSlope", terra::res(dem_kuamut)[1], "m.tif"))

  slp <- terra::terrain(dem_kuamut,
    v = "slope", unit = c("radians"),
    neighbors = 8, filename = slp_ras, overwrite = TRUE
  )

  points <- sf::read_sf(points_path)
  # extract slp for plot size
  points_radius <- sf::st_buffer(points, radius)

  kuamuts_sample_plots <- points_radius |>
    mutate(
      slope = exactextractr::exact_extract(slp, points_radius, fun = "mean"),
      area_planimetric = area * cos(slope), # we need planimetric area
      Area = round(area_planimetric / (100 * 100), 3)
    ) |> # add Area of Sample Plot (in Hectares) to data.frame
    dplyr::rename(PlotNo = Pl)

  # Save data file of points including Sample Plot Area
  plot_area_path <- set_path(sprintf("PlotData/%sSamplePlots_TC", site))

  spat_pth <- paste0(plot_area_path, ".fgb")
  csv_pth <- paste0(plot_area_path, ".csv")

  sf::write_sf(kuamuts_sample_plots, spat_pth, delete_dsn = TRUE)
  write.csv(kuamuts_sample_plots, csv_pth, row.names = FALSE)

  return(list(
    slp_correct_spat = spat_pth,
    slp_correct_csv = csv_pth,
    detail_slp_ras = slp_ras
  ))
}