mh_rep: Multivariate Climate Representativeness Analysis

In ClimaRep: Estimating Climate Representativeness

Multivariate Climate Representativeness Analysis

Description

This function calculates Mahalanobis-based Climate Representativeness for input polygon within a defined area.

Representativeness is assessed by comparing the multivariate climate conditions of each cell, of the reference climate space (climate_variables), with the climate conditions within each specific input polygon.

Usage

mh_rep(
  polygon,
  col_name,
  climate_variables,
  th = 0.95,
  dir_output = file.path(tempdir(), "ClimaRep"),
  save_raw = FALSE
)

Arguments

polygon	An sf object containing the defined areas. Must have the same CRS as climate_variables.
col_name	character. Name of the column in the polygon object that contains unique identifiers for each polygon.
climate_variables	A SpatRaster stack of climate variables representing the conditions. Its CRS will be used as the reference system.
th	numeric (0-1). Percentile threshold used to define representativeness. Cells with a Mahalanobis distance below or equal to the th are classified as representative (default: 0.95).
dir_output	character. Path to the directory where output files will be saved. The function will create subdirectories within this path.
save_raw	logical. If TRUE, saves the intermediate continuous Mahalanobis distance rasters calculated for each polygon before binary classification. The final binary classification rasters are always saved (default: FALSE).

Details

This function performs a multivariate analysis using Mahalanobis distance to assess the Climate Representativeness of input polygons for a single time period.

Crucially, this function assumes that all spatial inputs (polygon, climate_variables) are already correctly aligned and share the same Coordinate Reference System (CRS). If inputs do not meet these criteria, the function will stop with an informative error.

Here are the key steps:

1. Checking of spatial inputs: Ensures that polygon and climate_variables have matching CRSs.

2. Calculate the multivariate covariance matrix using climate data from all cells.

3. For each polygon in the polygon object:

   • Crop and mask the climate variables raster (climate_variables) to the boundary of the current polygon.

   • Calculate the multivariate mean using the climate data from the previous step. This defines the climate centroid for the current polygon.

   • Calculate the Mahalanobis distance for each cell relative to the centroid and covariance matrix.

   • Apply the specified threshold (th) to Mahalanobis distances to determine which cells are considered representative. This threshold is a percentile of the Mahalanobis distances within the current polygon.

   • Classify each cell as Representative = 1 (Mahalanobis distance \le th) or Non-Representative = 0 (Mahalanobis distance $>$ th).

4. Saves the binary classification raster (.tif) and generates a corresponding visualization map (.jpeg) for each polygon. These are saved within the specified output directory (dir_output).

It is important to note that Mahalanobis distance is sensitive to collinearity among variables. While the covariance matrix accounts for correlations, it is strongly recommended that the climate_variables are not strongly correlated. Consider performing a collinearity analysis beforehand, perhaps using the vif_filter function from this package.

Value

Writes the following outputs to disk within subdirectories of dir_output:

• Classification (.tif ) rasters: Binary rasters (0 for Non-representative and1 for Representative) for each input polygon are saved in the ⁠Representativeness/⁠ subdirectory.

• Visualization (.jpeg) maps: Image files visualizing the classification results for each polygon are saved in the ⁠Charts/⁠ subdirectory.

• Raw Mahalanobis distance rasters: Optionally saved as .tif files in the ⁠Mh_Raw/⁠ subdirectory if save_raw = TRUE.

Examples

library(terra)
library(sf)
set.seed(2458)
n_cells <- 100 * 100
r_clim_present <- terra::rast(ncols = 100, nrows = 100, nlyrs = 7)
values(r_clim_present) <- c(
   (terra::rowFromCell(r_clim_present, 1:n_cells) * 0.2 + rnorm(n_cells, 0, 3)),
   (terra::rowFromCell(r_clim_present, 1:n_cells) * 0.9 + rnorm(n_cells, 0, 0.2)),
   (terra::colFromCell(r_clim_present, 1:n_cells) * 0.15 + rnorm(n_cells, 0, 2.5)),
   (terra::colFromCell(r_clim_present, 1:n_cells) +
     (terra::rowFromCell(r_clim_present, 1:n_cells)) * 0.1 + rnorm(n_cells, 0, 4)),
   (terra::colFromCell(r_clim_present, 1:n_cells) /
     (terra::rowFromCell(r_clim_present, 1:n_cells)) * 0.1 + rnorm(n_cells, 0, 4)),
   (terra::colFromCell(r_clim_present, 1:n_cells) *
     (terra::rowFromCell(r_clim_present, 1:n_cells) + 0.1 + rnorm(n_cells, 0, 4))),
   (terra::colFromCell(r_clim_present, 1:n_cells) *
     (terra::colFromCell(r_clim_present, 1:n_cells) + 0.1 + rnorm(n_cells, 0, 4)))
)
names(r_clim_present) <- c("varA", "varB", "varC", "varD", "varE", "varF", "varG")
terra::crs(r_clim_present) <- "EPSG:4326"

vif_result <- ClimaRep::vif_filter(r_clim_present, th = 5)
print(vif_result$summary)
r_clim_present_filtered <- vif_result$filtered_raster
hex_grid <- sf::st_sf(
   sf::st_make_grid(
     sf::st_as_sf(
       terra::as.polygons(
         terra::ext(r_clim_present_filtered))),
     square = FALSE)
     )
sf::st_crs(hex_grid) <- "EPSG:4326"
polygons <- hex_grid[sample(nrow(hex_grid), 2), ]
polygons$name <- c("Pol_A", "Pol_B")
terra::plot(r_clim_present_filtered[[1]])
terra::plot(polygons, add = TRUE, color = "transparent", lwd = 3)

ClimaRep::mh_rep(
   polygon = polygons,
   col_name = "name",
   climate_variables = r_clim_present_filtered,
   th = 0.95,
   dir_output = file.path(tempdir(), "ClimaRep"),
   save_raw = TRUE
   )

ClimaRep documentation built on June 28, 2025, 1:07 a.m.