R/rep_overlay.R

In ClimaRep: Estimating Climate Representativeness

#' @title Overlay ClimaRep classifications
#'
#' @description Combines multiple single-layer rasters (`tif`), outputs from `ClimaRep::mh_rep()` or `ClimaRep::mh_rep_ch()` for different input polygons, into a multi-layered `SpatRaster`.
#'
#' This function handles inputs from both `ClimaRep::mh_rep()` (which primarily contains **Representative** cells) and `ClimaRep::mh_rep_ch()` (which includes **Retained**, **Lost**, and **Novel** cells). 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 `ClimaRep::mh_rep()` or `ClimaRep::mh_rep_ch()`. These rasters should primarily
#'  contain the categories: `1` (Retained/Representative), `2` (Lost), and `3` (Novel).
#'  Category `0` (Non-representative) will be ignored for the RGB output.
#' @param output_dir `character`. Path to the directory where the output file will be saved.
#'
#' @return Writes the following outputs within the directory specified by `output_dir`:
#' When `ClimaRep::mh_rep()` results are used, the output layers consistently represent counts for **Representative** categories across all input rasters.
#' When `ClimaRep::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.
#' \itemize{
#'  \item A multi-layered `SpatRaster` (`ClimaRep_overlay.tif`) for RGB visualization.
#'  \item Individual `.tif` files for each band (Lost, Retained, Novel) in `Individual_Bands/` subdirectory.
#'  \item `Additional note`: The `ClimaRep::mh_rep()` function analyzes a single period. When its output is used, representative cells are all categorized as Retained (value 1), while other categories (Lost and Novel) have a value of zero.
#' }
#'
#' @details
#' This function streamlines the aggregation of ClimaRep classifications. It is designed to work with outputs from both `ClimaRep::mh_rep()` and `ClimaRep::mh_rep_ch`.
#'
#' For each of the three key categories (Lost, Retained/representative, 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-number of `polygons`) for finer 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::rep_overlay(folder_path = system.file("extdata",
#'                                                                     package = "ClimaRep"),
#'                                          output_dir = file.path(tempdir(), "rep_overlay_output"))
#' terra::plotRGB(ClimaRep_overlay)
#' terra::plot(ClimaRep_overlay)
#' @export
rep_overlay <- function(folder_path,
                        output_dir = file.path(tempdir(), "ClimaRep_overlay")) {
  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.")
  }
  if (!is.character(output_dir) || length(output_dir) != 1) {
    stop("Parameter 'output_dir' must be a single character string.")
  }
  message("Establishing output file structure.")
  dir_individual_bands <- file.path(output_dir, "individual_bands")
  if (!dir.exists(output_dir)) {
    dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  }
  if (!dir.exists(dir_individual_bands)) {
    dir.create(dir_individual_bands, recursive = TRUE, showWarnings = FALSE)
  }
  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) {
    stop("No '.tif files' found in the specified folder: ", folder_path)
  }
  message("Processing ", length(raster_files), " rasters from ", folder_path)
  first_raster <- terra::rast(raster_files[1])
  count_rasters_for_rgb <- list(
    Lost = first_raster * 0,
    Retained = first_raster * 0,
    Novel = first_raster * 0)
  for (i in seq_along(raster_files)) {
    file_path <- raster_files[i]
    current_raster <- terra::rast(file_path)
    for (cat_name in names(rgb_category_map)) {
      category_value <- rgb_category_map[cat_name]
      count_rasters_for_rgb[[cat_name]] <- count_rasters_for_rgb[[cat_name]] +
        terra::ifel(current_raster == category_value, 1, 0)
    }
  }
  for (cat_name in names(rgb_category_map)) {
    band_filename <- paste0(rgb_channel_names[which(names(rgb_category_map) == cat_name)], ".tif")
    terra::writeRaster(
      count_rasters_for_rgb[[cat_name]],
      file.path(dir_individual_bands, band_filename),
      overwrite = TRUE,
      datatype = "INT2U")
  }
  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
  dir_output_file <- file.path(output_dir, "ClimaRep_overlay.tif")
  terra::writeRaster(final_rgb_stack,
                     dir_output_file,
                     overwrite = TRUE,
                     datatype = "INT2U")
  message("All processes were completed")
  message(paste("Output files saved in: ", output_dir))
  return(invisible(final_rgb_stack))
}