ClimaRep (version 1.0)

mh_rep: Multivariate climate representativeness analysis

Description

This function calculates Mahalanobis-based climate representativeness for input polygon within a defined area.

The function categorizes cells based on how their climate representativeness, labeling them as Non-representative and Representative.

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,
  study_area,
  th = 0.95,
  dir_output = file.path(tempdir(), "ClimaRep"),
  save_raw = FALSE
)

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.

Arguments

polygon

An sf object containing the defined areas. Its CRS will be used as the reference system.

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.

study_area

An sf object defining the study area. The analysis will be limited to this area.

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.

Key steps:

1. Ensures that climate_variables and study_area have matching CRSs with polygon by automatically transforming them if needed. The CRS of polygon will be used as the reference system.

2. Crops and masks the climate_variables raster to the study_area to limit all subsequent calculations to the area of interest.

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

4. 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).

5. 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.

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")
study_area_polygon <- sf::st_as_sf(terra::as.polygons(terra::ext(r_clim_present_filtered)))
sf::st_crs(study_area_polygon) <- "EPSG:4326"
terra::plot(r_clim_present_filtered[[1]])
terra::plot(polygons, add = TRUE, color = "transparent", lwd = 3)
terra::plot(study_area_polygon, add = TRUE, col = "transparent", lwd = 3, border = "red")

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