R/mh_overlay.R

In ClimaRep: Estimating Climate Representativeness

#' @title Overlay Mahalanobis-based Climate Representativeness Classifications
#'
#' @description Combines multiple single-layer rasters (`tif`), outputs from `mh_rep` or `mh_rep_ch` for different input polygons, into a multi-layered `SpatRaster`.
#'
#' This function handles inputs from both `mh_rep` (which primarily contains **Represented** areas) and `mh_rep_ch` (which includes **Retained**, **Lost**, and **Novel** areas). The output layers consistently represent counts of each input.
#'
#' @param folder_path `character`. The path to the directory containing the classification rasters (`.tif`) generated by `mh_rep` or `mh_rep_ch`. These rasters should primarily
#'  contain the categories: `1` (Retained/Represented), `2` (Lost), and `3` (Novel).
#'  Category `0` (Non-represented) will be ignored for the RGB output.
#'
#' @return Writes the multi-layered (`ClimaRep_overlay.tif`) outputs to disk in a new `overlay` subfolder within the `folder_path`.
#' When `mh_rep_ch` results are used, the output layers consistently represent counts for **Lost** (Red), **Retained** (Green), and **Novel** (Blue) categories across all input rasters. Designed for direct RGB plotting.
#' When `mh_rep` results are used, the output layers consistently represent counts for **Represented** categories across all input rasters.
#'
#' @details
#' This function streamlines the aggregation of Climate Representativeness classifications. It is designed to work with outputs from both `mh_rep` and `mh_rep_ch`.
#'
#' For each of the three key categories (Lost, Retained/Represented, Novel), the function:
#' \enumerate{
#'  \item Identifies and reads all `.tif` files within the `folder_path`.
#'  \item For each input raster, it creates a binary layer: `1` if the cell's value matches the target category (e.g., `2` for 'Lost'), and `0` otherwise.
#'  \item Sums these binary layers to generate a cumulative count for that specific category at each grid cell.
#' }
#'
#' The three resulting count layers (Lost, Retained, Novel) are then consistently stacked in the following order:
#' \itemize{
#'  \item First layer (Red): Cumulative count of **Lost**.
#'  \item Second layer (Green): Cumulative count of **Retained**.
#'  \item Third layer (Blue): Cumulative count of **Novel**.
#' }
#' This fixed order ensures that the output `SpatRaster` is immediately ready for direct RGB visualization using `terra::plotRGB()`, where the color mixtures will intuitively reflect
#' the spatial agreement of these change types.
#'
#' The output `SpatRaster` contains raw counts. While `terra::plotRGB()` often handles stretching for visualization, users might normalize these counts manually (e.g., to 0-255) for finer control over visual contrast.
#'
#' A new subfolder named `overlay/` will be created within the `folder_path`. The resulting three-layered RGB will be saved as `ClimaRep_overlay.tif` inside this new `overlay/` subfolder.
#'
#' @importFrom terra rast ifel app writeRaster nlyr values
#'
#' @examples
#' ClimaRep_overlay <- ClimaRep::mh_overlay(folder_path = system.file("extdata", package = "ClimaRep"))
#' terra::plotRGB(ClimaRep_overlay)
#' terra::plot(ClimaRep_overlay)
#' @export
mh_overlay <- function(folder_path) {
  if (!is.character(folder_path) ||
      length(folder_path) != 1 || !dir.exists(folder_path)) {
    stop("Parameter 'folder_path' must be a character string and a valid directory.")
  }
  rgb_category_map <- c("Lost" = 2,
                        "Retained" = 1,
                        "Novel" = 3)
  rgb_channel_names <- c("Lost_Count_R", "Retained_Count_G", "Novel_Count_B")
  raster_files <- list.files(folder_path, pattern = "\\.tif$", full.names = TRUE)
  if (length(raster_files) == 0) {
    warning("No '.tif files' found in the specified folder: ",
            folder_path,
            ". Returning NULL.")
    return(invisible(NULL))
  }
  message(
    "Processing ",
    length(raster_files),
    " rasters from ",
    folder_path
  )
  first_raster <- terra::rast(raster_files[1])
  count_rasters_for_rgb <- vector("list", length = length(rgb_category_map))
  names(count_rasters_for_rgb) <- names(rgb_category_map)
  all_binary_layers <- list()
  for (cat_name in names(rgb_category_map)) {
    category_value <- rgb_category_map[cat_name]
    message("Calculating counts for category: ",
            cat_name,
            " (value = ",
            category_value,
            ") ")
    current_category_binary_layers <- list()
    for (i in seq_along(raster_files)) {
      file_path <- raster_files[i]
      current_raster <- terra::rast(file_path)
      binary_layer <- terra::ifel(current_raster == category_value, 1, 0)
      current_category_binary_layers[[length(current_category_binary_layers) + 1]] <- binary_layer
    }
    if (length(current_category_binary_layers) > 0) {
      stacked_binary_layers <- terra::rast(current_category_binary_layers)
      current_category_sum_raster <- terra::app(stacked_binary_layers, fun =
                                                  "sum")
    } else {
      current_category_sum_raster <- first_raster * 0
    }
    count_rasters_for_rgb[[cat_name]] <- current_category_sum_raster
  }
  final_rgb_stack <- c(count_rasters_for_rgb[["Lost"]],
                       count_rasters_for_rgb[["Retained"]],
                       count_rasters_for_rgb[["Novel"]])
  names(final_rgb_stack) <- rgb_channel_names
  overlay_dir <- file.path(folder_path, "overlay")
  if (!dir.exists(overlay_dir)) {
    dir.create(overlay_dir,
               recursive = TRUE,
               showWarnings = FALSE)
  }
  dir_output <- file.path(overlay_dir, "ClimaRep_overlay.tif")
  terra::writeRaster(final_rgb_stack,
                     dir_output,
                     overwrite = TRUE,
                     datatype = "INT2U")
  message("All processes were completed")
  message(paste("Output files in: ", dir_output))
  return(invisible(final_rgb_stack))
}