Title:	Multiscale Systematic Conservation Planning Across Nested H3 Grids
Version:	0.1.1
Maintainer:	Pablo Merlo <pablo.merlo@universityofgalway.ie>
Description:	Provides tools for multiscale systematic conservation planning using the H3 hierarchical hexagonal grid system (Uber Technologies (2024) https://h3geo.org) and the 'prioritizr' package (Hanson et al. (2025) <doi:10.1111/cobi.14376>). Supports the definition and solution of conservation problems across nested H3 resolutions with resolution-specific features, costs, and management attributes, including cross-scale connectivity penalties derived from parent-child relationships. Also includes utilities to evaluate solutions using multiscale-aware diagnostics and to post-process optimization outputs into alternative area-targeted conservation scenarios.
Depends:	R (≥ 4.2)
Imports:	sf, terra, raster, dplyr, tibble, Matrix, methods, prioritizr, exactextractr, stats, cli, R6, assertthat, h3jsr, rlang
Suggests:	ggplot2, rnaturalearth, rnaturalearthdata, rmapshaper, knitr, rmarkdown, testthat (≥ 3.0.0)
VignetteBuilder:	knitr
License:	MIT + file LICENSE
Encoding:	UTF-8
RoxygenNote:	7.3.3
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2026-03-25 14:35:29 UTC; merlo
Author:	Pablo Merlo [aut, cre], Oisín Callery [aut], Carlos Tighe [ctb], Anthony Grehan [ctb]
Repository:	CRAN
Date/Publication:	2026-03-30 17:00:02 UTC

Add multiscale connectivity penalties to a prioritizr problem

Description

This function mirrors prioritizr::add_connectivity_penalties() but adds a second, independent symmetric penalty for cross-resolution (vertical) connectivity between H3 planning units. It accepts the same input formats (matrix, Matrix::dgCMatrix, data frame, or 4D array) and internally converts them to a sparse connectivity matrix.

Usage

add_multiscale_connectivity_penalties(
  x,
  penalty,
  zones = diag(prioritizr::number_of_zones(x)),
  data,
  normalize = c("none", "sym")
)

## S4 method for signature 'ANY,ANY,ANY,matrix,character'
add_multiscale_connectivity_penalties(
  x,
  penalty,
  zones = diag(prioritizr::number_of_zones(x)),
  data,
  normalize = c("none", "sym")
)

## S4 method for signature 'ANY,ANY,ANY,Matrix,character'
add_multiscale_connectivity_penalties(
  x,
  penalty,
  zones = diag(prioritizr::number_of_zones(x)),
  data,
  normalize = c("none", "sym")
)

## S4 method for signature 'ANY,ANY,ANY,data.frame,character'
add_multiscale_connectivity_penalties(
  x,
  penalty,
  zones = diag(prioritizr::number_of_zones(x)),
  data,
  normalize = c("none", "sym")
)

## S4 method for signature 'ANY,ANY,ANY,dgCMatrix,character'
add_multiscale_connectivity_penalties(
  x,
  penalty,
  zones = diag(prioritizr::number_of_zones(x)),
  data,
  normalize = c("none", "sym")
)

## S4 method for signature 'ANY,ANY,ANY,array,character'
add_multiscale_connectivity_penalties(
  x,
  penalty,
  zones = diag(prioritizr::number_of_zones(x)),
  data,
  normalize = c("none", "sym")
)

Arguments

Value

The modified ConservationProblem object.

Functions

• add_multiscale_connectivity_penalties( x = ANY, penalty = ANY, zones = ANY, data = matrix, normalize = character ): matrix method

• add_multiscale_connectivity_penalties( x = ANY, penalty = ANY, zones = ANY, data = Matrix, normalize = character ): Matrix method

• add_multiscale_connectivity_penalties( x = ANY, penalty = ANY, zones = ANY, data = data.frame, normalize = character ): data.frame (Marxan) method

• add_multiscale_connectivity_penalties( x = ANY, penalty = ANY, zones = ANY, data = dgCMatrix, normalize = character ): dgCMatrix method

• add_multiscale_connectivity_penalties( x = ANY, penalty = ANY, zones = ANY, data = array, normalize = character ): array (4D) method

Build cross-scale index structures for H3-based SCP workflows

Description

Extends the basic hierarchy from build_h3_maps() into full ancestor, descendant and resolution-index mappings used by all multiscale selection and evaluation functions.

Usage

build_crossscale_index(maps)

Arguments

Value

A list with elements:

• res_levels

• rows_by_res

• pos_in_res

• anc_at_res

• desc_at_res

• finer_rows_by_r0cell

Examples

h3_child  <- "8a2a1072b59ffff"
h3_parent <- "872a1072bffffff"

maps <- build_h3_maps(
  s_or_h3 = c(h3_parent, h3_child),
  res_vec = c(7L, 8L)
)
cs_idx <- build_crossscale_index(maps)

names(cs_idx)
cs_idx$res_levels

Build H3 hierarchy maps for multiscale planning units

Description

Creates basic parent–child relationships for a multiresolution H3 planning-unit dataset. This is the core structure used by all cross-scale operations.

Usage

build_h3_maps(s_or_h3, res_vec = NULL, res_levels = NULL)

Arguments

Value

A list with elements h3_vec, res_vec, res_levels, row_idx_by_h3, nearest_parent_row_of, and children_by_row.

Examples

# Minimal example: two resolutions, parent-child relationship
h3_child  <- "8a2a1072b59ffff" # example H3 index
h3_parent <- "872a1072bffffff"

maps <- build_h3_maps(
  s_or_h3  = c(h3_parent, h3_child),
  res_vec  = c(7L, 8L),
  res_levels = c(7L, 8L)
)

str(maps, max.level = 1)
maps$nearest_parent_row_of

Build a multiscale H3 connectivity matrix

Description

Construct a symmetric sparse connectivity matrix linking each finer H3 planning unit to its parent at the next coarser resolution. This is typically used as the data input to add_multiscale_connectivity_penalties().

Usage

build_multiscale_connectivity_matrix(maps, symmetric = TRUE)

Arguments

Value

A Matrix::dgCMatrix connectivity matrix of size ⁠n_pu × n_pu⁠.

Examples

# Minimal 2-resolution parent-child example
h3_child  <- "8a2a1072b59ffff"
h3_parent <- "872a1072bffffff"

maps <- build_h3_maps(
  s_or_h3 = c(h3_parent, h3_child),
  res_vec = c(7L, 8L)
)

conn <- build_multiscale_connectivity_matrix(maps)
conn

Count cross-resolution overlaps in a selection

Description

For a given 0/1 selection column, counts how many planning units are simultaneously selected at a finer resolution and at their ancestor cell at a coarser resolution (using the hierarchy in maps).

For each pair of resolutions r_low < r_high, it reports the number of co-selected ancestor–descendant pairs, plus a total across all pairs.

Usage

compute_overlaps_by_resolution(
  s,
  flag_col = "selected_by_resolution_final",
  maps
)

Arguments

Value

A named integer vector with one entry per resolution pair (e.g. "5-6", "6-7") and a "total" element.

Examples

parent <- "872a1072bffffff"
kids   <- c("882a1072b1fffff", "882a1072b3fffff")

s <- data.frame(
  h3_address = c(parent, kids),
  res = c(7L, 8L, 8L),
  selected_by_resolution_final = c(1L, 1L, 0L)
)
maps <- build_h3_maps(s, res_levels = c(7L, 8L))

compute_overlaps_by_resolution(s, "selected_by_resolution_final", maps)

Compute importance scores for planning units

Description

This function calculates a continuous importance score (ensemble_score) for each planning unit by combining two sources of information: (1) a rarity-weighted aggregation of feature values, and (2) how consistently each planning unit is selected across one or more solutions.

The relative influence of these components is controlled by alpha_freq. When multiple solutions are provided, selection consistency reflects how often each planning unit is selected across the portfolio. When a single solution is provided, it reflects whether the planning unit is selected or not. If no solution columns are present, the score is based entirely on the rarity-weighted feature component.

Optionally, locked_in planning units can be forced to rank above all others.

Usage

compute_pu_scores(
  s,
  feature_mat,
  feature_weights = NULL,
  alpha_freq = 0.25,
  freq_gamma = 1.5,
  lock_mode = c("top", "none"),
  winsor_p = NULL,
  cost_col = NULL,
  cost_beta = 1
)

Arguments

Value

The same sf object s with two new numeric columns: selection_consistency and ensemble_score.

Examples

# Minimal sf with 3 planning units + two solution columns
s <- sf::st_as_sf(
  data.frame(
    x = c(0, 1, 2),
    y = c(0, 0, 0),
    solution_1 = c(1, 0, 1),
    solution_2 = c(1, 1, 0),
    locked_in  = c(FALSE, TRUE, FALSE)
  ),
  coords = c("x", "y"), crs = 4326
)

# Feature matrix must match nrow(s)
feature_mat <- data.frame(
  f1 = c(0.2, 0.0, 0.8),
  f2 = c(0.0, 0.5, 0.1)
)

s2 <- compute_pu_scores(
  s, feature_mat,
  alpha_freq = 0.25,
  lock_mode  = "top"
)

s2[, c("selection_consistency", "ensemble_score", "locked_in")]

Compute cross-scale planning-unit scores by resolution

Description

This function computes a feature-based importance score at each resolution present in a multiscale planning-unit set, and then propagates those resolution-specific scores through the H3 hierarchy so every planning unit receives a score "as seen from" each resolution.

Propagation rules:

• Downward (parent \rightarrow children): undamped broadcast. All descendants of a given parent at resolution r inherit the same score_from_r value.

• Upward (children \rightarrow parent): aggregation of descendant scores at resolution r using mean (default) or max.

Usage

compute_pu_scores_crossscale(
  s,
  feature_mat,
  maps,
  cs_idx,
  feature_weights = NULL,
  winsor_p = NULL,
  cost_col = NULL,
  cost_beta = 1,
  agg_mode = c("mean", "max"),
  res_col = "res",
  res_levels = NULL,
  feature_cols_by_res = NULL,
  prefix = "score_from_r"
)

Arguments

Value

The same sf object s with one new numeric column per resolution, named paste0(prefix, res) (e.g., score_from_r5).

Examples

# Tiny 2-level setup: 2 parents (r7), and 2 children (r8) under the first parent
parent1 <- "872a1072bffffff"
parent2 <- "872a10729ffffff"
kids    <- c("882a1072b1fffff", "882a1072b3fffff")

h3_vec  <- c(parent1, parent2, kids)
res_vec <- c(7L, 7L, 8L, 8L)

s <- sf::st_as_sf(
  data.frame(
    h3_address = h3_vec,
    res        = res_vec,
    cost       = c(1, 1, 1, 1),
    x          = c(0, 1, 2, 3),
    y          = c(0, 0, 0, 0)
  ),
  coords = c("x", "y"), crs = 4326
)

# Feature columns prefixed by resolution (r7_, r8_)
feature_mat <- data.frame(
  r7_f1 = c(1.0, 0.2, 0.0, 0.0),
  r8_f1 = c(0.0, 0.0, 0.6, 0.2)
)

maps   <- build_h3_maps(s)
cs_idx <- build_crossscale_index(maps)

s2 <- compute_pu_scores_crossscale(
  s, feature_mat,
  maps = maps, cs_idx = cs_idx,
  agg_mode = "mean",
  res_col = "res"
)

s2[, c("res", "score_from_r7", "score_from_r8")]

Stratified multiscale selection by resolution

Description

Select planning units to meet area-based targets at each H3 resolution. Selection is based on a single score column (e.g., compute_pu_scores) and a cross-scale index produced by build_crossscale_index.

Resolutions are processed from finest to coarsest. At each resolution, the function accounts for area already covered by selected finer units to avoid double-counting. To keep selections non-overlapping, it prefers units that do not have finer descendants in the input dataset; coarser units with descendants are only used if needed to reach the target.

The result is a resolution-stratified selection that meets per-resolution area targets while minimizing overlap between coarse and fine planning units.

Usage

compute_selection_by_resolution(
  s,
  cs_idx,
  target = 0.3,
  score_col = "ensemble_score",
  area_col = "area_km2",
  res_col = "res",
  out_col = "selected_by_resolution",
  blocked_col = NULL,
  target_by_res = NULL
)

Arguments

Value

The input s with an additional column named by out_col (0/1), and optionally a column named by blocked_col if requested.

Examples

# Build a tiny multiscale set: 2 parents (r7), 2 children (r8) under parent1
parent1 <- "872a1072bffffff"
parent2 <- "872a10729ffffff"
kids1   <- c("882a1072b1fffff", "882a1072b3fffff")

h3_vec  <- c(parent1, parent2, kids1)
res_vec <- c(7L, 7L, 8L, 8L)

s <- data.frame(
  h3_address     = h3_vec,
  res            = res_vec,
  area_km2       = c(14, 14, 2, 2),
  ensemble_score = c(0.2, 0.9, 0.8, 0.1)
)

maps   <- build_h3_maps(s, res_levels = c(7L, 8L))
cs_idx <- build_crossscale_index(maps)

out <- compute_selection_by_resolution(
  s, cs_idx,
  target    = 0.5,
  score_col = "ensemble_score",
  area_col  = "area_km2",
  res_col   = "res",
  out_col   = "selected_by_resolution"
)

out[, c("res", "area_km2", "ensemble_score", "selected_by_resolution")]

Stratified multiscale selection by strata and resolution

Description

Select planning units to meet area-based targets for each strata (e.g. country) at each H3 resolution. Selection is based on a single score column (typically ensemble_score) and a cross-scale index produced by build_crossscale_index.

Resolutions are processed from finest to coarsest. At each strata and resolution, the function accounts for area already covered by selected planning units at other resolutions to avoid double-counting. Within each strata, it prefers native-resolution units that do not have finer descendants in the input dataset, using coarser units with descendants only when needed to reach the target.

The result is a single 0/1 selection that meets per-strata area targets while minimizing cross-scale overlap.

Usage

compute_selection_by_strata(
  s,
  cs_idx,
  strata_masks,
  target = 0.3,
  admin_strata = "admin",
  score_col = "ensemble_score",
  area_col = "area_km2",
  res_col = "res",
  out_col = "selected_by_admin",
  blocked_col = NULL,
  target_by_strata = NULL
)

Arguments

Value

The input s with an additional column out_col (0/1), and optionally blocked_col if requested.

Examples

# Tiny multiscale example with two strata (A and B)
parent1 <- "872a1072bffffff"
kids1   <- c("882a1072b1fffff", "882a1072b3fffff")
parent2 <- "872a10729ffffff"
parent3 <- "872a10774ffffff"

h3_vec  <- c(parent1, parent2, parent3, kids1)
res_vec <- c(7L, 7L, 7L, 8L, 8L)

s <- data.frame(
  h3_address     = h3_vec,
  res            = res_vec,
  admin          = c("A", "B", "B", "A", "A"),
  area_km2       = c(14, 14, 14, 2, 2),
  ensemble_score = c(0.2, 0.9, 0.2, 0.8, 0.1)
)

maps   <- build_h3_maps(s, res_levels = c(7L, 8L))
cs_idx <- build_crossscale_index(maps)

strata_masks <- list(
  A = (s$admin == "A"),
  B = (s$admin == "B")
)

out <- compute_selection_by_strata(
  s, cs_idx,
  strata_masks     = strata_masks,
  target           = 0.3,
  admin_strata     = "admin",
  score_col        = "ensemble_score",
  area_col         = "area_km2",
  res_col          = "res",
  out_col          = "selected_by_admin",
  target_by_strata = c(A = 0.12, B = 0.5)
)

out[, c("admin", "res", "area_km2", "ensemble_score", "selected_by_admin")]

Deduplicate multiscale selections across H3 resolutions

Description

Removes redundant planning-unit selections across nested H3 resolutions. Because finer and coarser H3 cells overlap perfectly (parent–child hierarchy), a selection at one resolution makes selections at other resolutions redundant.

This function enforces a consistent hierarchy using one of two strategies:

• "coarser_first" – keep coarser selected cells and drop all selected descendants (finer cells). Useful if coarse-scale representation should dominate.

• "finer_first" – keep finer selected cells and drop selected ancestors (coarser cells). Useful when fine-scale detail should dominate and coarse cells should not “override” them.

Usage

deduplicate_h3_selections(
  s,
  sel_col,
  h3_vec,
  res_vec,
  res_levels,
  nearest_parent_row_of,
  children_by_row,
  mode = c("coarser_first", "finer_first")
)

Arguments

Details

When mode = "finer_first", removing selected coarser cells can reduce the total selected area if only a subset of their descendant cells were selected (i.e., partial coverage of the parent footprint).

The function operates in a single pass, using the precomputed parent/child lists from build_h3_maps(), and produces a new column paste0(sel_col, "_deduplicated").

Value

The input s with an additional column paste0(sel_col, "_deduplicated") containing a clean 0/1 selection.

Examples

# One parent (res7) with two children (res8)
parent <- "872a1072bffffff"
kids   <- c("882a1072b1fffff", "882a1072b3fffff")

h3_vec  <- c(parent, kids)
res_vec <- c(7L, 8L, 8L)

maps <- build_h3_maps(h3_vec, res_vec = res_vec)

s <- data.frame(solution_1 = c(1L, 1L, 1L))

# Keep the coarser cell (drops children)
out_coarse <- deduplicate_h3_selections(
  s, sel_col = "solution_1",
  h3_vec = maps$h3_vec, res_vec = maps$res_vec, res_levels = maps$res_levels,
  nearest_parent_row_of = maps$nearest_parent_row_of,
  children_by_row = maps$children_by_row,
  mode = "coarser_first"
)

# Keep the finer cells (drops parent)
out_fine <- deduplicate_h3_selections(
  s, sel_col = "solution_1",
  h3_vec = maps$h3_vec, res_vec = maps$res_vec, res_levels = maps$res_levels,
  nearest_parent_row_of = maps$nearest_parent_row_of,
  children_by_row = maps$children_by_row,
  mode = "finer_first"
)

out_coarse
out_fine

Exact feature coverage from rasters and selected PUs

Description

Calculate how well raster-based features are represented by a selected set of planning units.

This function evaluates feature coverage from raster layers using exact area-weighted extraction. For each feature, it calculates the total amount available within its native-resolution planning unit footprint, the amount held by the selected planning units, and the proportion of the feature represented by the solution.

If targets are provided, the function also reports absolute and relative shortfalls from those targets.

Usage

eval_exact_raster_coverage(
  spp_all,
  pu_sf,
  solution_col,
  targets = NULL,
  res_col = "res"
)

Arguments

Value

A data frame with columns:

• feature

• native_res

• total_area Total feature amount within the native-resolution domain (area-integral units).

• held_area Feature amount held by the selected area, restricted to that domain (area-integral units).

• relative_held Relative held amount (held_area / total_area).

• target Relative target proportion (if targets supplied).

• absolute_shortfall Absolute shortfall in feature units (\max(0, target * total\_area - held\_area)).

• relative_shortfall Relative shortfall in proportions (\max(0, target - relative\_held)).

Examples

# Tiny raster with a single layer named with r7_ prefix
r <- terra::rast(nrows = 2, ncols = 2, xmin = 0, xmax = 2, ymin = 0, ymax = 2)
terra::values(r) <- 1
names(r) <- "r7_featA"

# One square PU covering the raster
pu <- sf::st_sf(
  res = 7L,
  selected = 1L,
  geometry = sf::st_sfc(
    sf::st_polygon(list(rbind(c(0, 0), c(2, 0), c(2, 2), c(0, 2), c(0, 0)))),
    crs = 4326
  )
)

eval_exact_raster_coverage(
  spp_all = r,
  pu_sf = pu,
  solution_col = "selected",
  targets = c(r7_featA = 0.3),
  res_col = "res"
)

Geometry-based feature coverage evaluation from sf column features and selected PUs

Description

Calculate how well features are represented by a selected set of planning units.

This function evaluates feature coverage using feature values stored directly in planning unit attributes. For each feature, it calculates the total amount available, the amount held by the selected planning units, and the proportion of the feature represented by the solution.

If targets are provided, the function also reports absolute and relative shortfalls from those targets.

Feature coverage is calculated using fractional overlap between planning units and the union of selected planning units.

Usage

eval_geom_feature_coverage(
  s,
  flag_col = "selected_by_resolution",
  feature_cols,
  res_col = "res",
  targets = NULL
)

Arguments

Value

A tibble with one row per feature and columns:

• feature, native_res

• total Total feature amount across all PUs at the native resolution.

• held Feature amount held by the selected area (via fractional overlap), restricted to the native-resolution footprint.

• rel_held Relative held amount (held / total).

• target Relative target proportion (if targets supplied).

• abs_shortfall Absolute shortfall in feature units (\max(0, target * total - held)).

• rel_shortfall Relative shortfall in proportions (\max(0, target - rel\_held)).

Examples

# Minimal polygons
p1 <- sf::st_polygon(list(rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1), c(0, 0))))
p2 <- sf::st_polygon(list(rbind(c(1, 0), c(2, 0), c(2, 1), c(1, 1), c(1, 0))))

s <- sf::st_sf(
  res = c(7L, 7L),
  area_km2 = c(1, 1),
  selected_by_resolution = c(1L, 0L),
  r7_featA = c(10, 5),
  geometry = sf::st_sfc(p1, p2, crs = 4326)
)

eval_geom_feature_coverage(
  s,
  flag_col = "selected_by_resolution",
  feature_cols = "r7_featA",
  res_col = "res"
)

Cross-scale coverage summary by resolution

Description

Summarize how much area is selected at each resolution in a multiscale H3 system. This function adjusts for overlap between coarse and fine planning units so that selected area is not double-counted across resolutions.

The output reports, for each resolution, the total available area, the effective selected area (after cross-scale adjustment), and the proportion selected.

Usage

summarize_coverage_by_resolution(
  s,
  flag_col = "selected_by_resolution_final",
  cs_idx
)

Arguments

Value

A tibble with columns:

• res – resolution,

• area_total_km2 – total area at that resolution,

• area_selected_km2 – cross-scale selected area, and

• coverage_prop – proportion selected.

Examples

parent <- "872a1072bffffff"
kids   <- c("882a1072b1fffff", "882a1072b3fffff")

s <- data.frame(
  h3_address = c(parent, kids),
  res = c(7L, 8L, 8L),
  area_km2 = c(14, 2, 2),
  selected_by_resolution_final = c(0L, 1L, 0L)
)
maps   <- build_h3_maps(s, res_levels = c(7L, 8L))
cs_idx <- build_crossscale_index(maps)

summarize_coverage_by_resolution(s, "selected_by_resolution_final", cs_idx)

Cross-scale coverage summary by strata

Description

Summarize how much area is selected within each stratum (e.g. country or region) in a multiscale H3 system. This function reports a cross-scale coverage metric that avoids double-counting when selected planning units overlap across resolutions.

For each stratum, coverage is reported at a single "native" resolution chosen from the planning units labeled with that stratum.

Usage

summarize_coverage_by_strata(
  s,
  flag_col = "selected_by_resolution_final",
  cs_idx,
  strata_col = "strata"
)

Arguments

Value

A tibble with one row per stratum and the following columns:

• strata

• native_res

• area_total_native_km2

• area_selected_km2

• coverage_prop

Examples

parent <- "872a1072bffffff"
kids   <- c("882a1072b1fffff", "882a1072b3fffff")

s <- data.frame(
  h3_address = c(parent, kids),
  res = c(7L, 8L, 8L),
  strata = c("A", "B", "B"),
  area_km2 = c(14, 2, 2),
  selected_by_resolution_final = c(0L, 1L, 0L)
)
maps   <- build_h3_maps(s, res_levels = c(7L, 8L))
cs_idx <- build_crossscale_index(maps)

summarize_coverage_by_strata(
  s,
  "selected_by_resolution_final",
  cs_idx,
  strata_col = "strata"
)