MultiscaleSCP extends prioritizr to
support hierarchical (multiresolution) planning units based on the H3
grid.
What is H3?
H3 is a global hierarchical hexagonal grid system developed by Uber. It
divides the Earth into nested hexagonal cells at discrete resolutions.
Each cell has a unique index that encodes its place in the hierarchy, so
every cell has:
This explicit parent–child structure makes H3 convenient for multiscale conservation planning because it lets you formally link planning units at different spatial resolutions through lineage relationships (ancestor ↔︎ descendant).
MultiscaleSCP leverages these links to:
A typical workflow follows this structure:
Planning Units (sf)
↓
build_h3_maps()
↓
build_crossscale_index()
↓
build_multiscale_connectivity_matrix()
↓
prioritizr problem + add_multiscale_connectivity_penalties()
↓
solve()
↓
post-processing (overlap diagnostics, deduplication, scoring, scenarios, coverage, feature evaluation)This single example illustrates an end-to-end multiscale workflow on a tiny artificial dataset. To keep the vignette fast and solver-independent, we do not run an optimizer here. Instead, we fabricate small 0/1 “solution” columns to demonstrate the post-hoc tools (deduplication, scenarios, coverage, and feature evaluation).
# tiny multiscale H3 system (3 parents at resolution 7; 2 kids at resolution 8 under parent1)
parent1 <- "87182e586ffffff"
kids1 <- unlist(h3jsr::get_children(parent1, 8, simple = TRUE), use.names = FALSE)[1:2]
parent2 <- "87182e584ffffff"
parent3 <- "87182e5b3ffffff"
h3_vec <- c(parent1, parent2, parent3, kids1)
res_vec <- c(7L, 7L, 7L, 8L, 8L) #resolution vector
# Create polygon geometries from H3 cells
geom <- do.call(c, lapply(h3_vec, function(x) h3jsr::cell_to_polygon(x, simple = TRUE)))
df <- data.frame(
h3_address = h3_vec,
res = res_vec,
area_km2 = c(14,14,14,2,2), #not true H3 areas
admin = c("A", "B", "B", "C", "C") #e.g., 2 countries and 1 county within one of them
)
s <- sf::st_sf(df, geometry = geom) |> sf::st_set_crs(4326)
# Feature matrix (resolution-prefixed)
# In a real workflow, these are typically PU-level feature relative coverage of each polygon.
feature_mat <- data.frame(
r7_f1 = c(1.0, 0.2, 0.0, 0.0, 0.0), # only meaningful at r7 rows
r8_f1 = c(0.0, 0.0, 0.0, 0.6, 0.2) # only meaningful at r8 rows
)
# Also attach features as columns (so eval_geom_feature_coverage() can use them directly)
s <- dplyr::bind_cols(s, feature_mat)
# Build hierarchy, cross-scale index, and multiscale connectivity matrix
maps <- build_h3_maps(s, res_levels = c(7L, 8L)) #`maps` stores ancestry mappings.
cs_idx <- build_crossscale_index(maps) #`cs_idx` stores fast row-level indices used in scoring and selection algorithms.
conn <- build_multiscale_connectivity_matrix(maps)
conn #This sparse matrix encodes parent–child connections across resolutions.
#> 5 x 5 sparse Matrix of class "dgCMatrix"
#>
#> [1,] . . . 1 1
#> [2,] . . . . .
#> [3,] . . . . .
#> [4,] 1 . . . .
#> [5,] 1 . . . .
# prioritizr problem (not solving to keep vignette fast)
#p <- problem(s, features = colnames(feature_mat))... #create a problem
#`add_multiscale_connectivity_penalties()` encourages co-selection between connected parent–child cells: increasing the penalty encourages more nested, scale-consistent solutions.
# p <- add_multiscale_connectivity_penalties(
# p, penalty = 1, data = conn)
# s <- solve(p) # any prioritizr-compatible solver may be used
# Fabricate a small "portfolio" of solutions (0/1 columns) to demonstrate post-hoc tools
# In a real workflow these columns would come from solve(p).
# solution_1 co-selects parent1 and kid1 => overlap exists
s$solution_1 <- c(1L, 0L, 0L, 1L, 0L)
s$solution_2 <- c(0L, 0L, 1L, 0L, 1L)
# Overlap diagnostics on raw selections (counts co-selected ancestor–descendant pairs)
compute_overlaps_by_resolution(
s = s,
flag_col = "solution_1",
maps = maps
)
#> 7-8 total
#> 1 1
# Deduplicate overlaps:
# Multiscale solutions may co-select both a coarse-resolution cell and one or more of its finer descendants; `deduplicate_h3_selections()` resolves these overlaps by retaining only one representative planning unit per location according to a chosen rule (i.e., `"finer_first"` or `"coarser_first"`).
# "finer_first": keep children; drop overlapping ancestors (the resulting solution may have less total area than the original solution if not all the children were selected)
s_finer <- deduplicate_h3_selections(
s = 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"
)
# "coarser_first": keep parent; drop overlapping children (total area is the same as the deduplicated solution)
s_coarser <- deduplicate_h3_selections(
s = 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"
)
# Compare what is kept (new column solution_1_deduplicated added, solution_1 column remain untouched)
sf::st_drop_geometry(s_finer)[, c("h3_address","res","solution_1","solution_1_deduplicated")]
#> h3_address res solution_1 solution_1_deduplicated
#> 1 87182e586ffffff 7 1 0
#> 2 87182e584ffffff 7 0 0
#> 3 87182e5b3ffffff 7 0 0
#> 4 88182e5861fffff 8 1 1
#> 5 88182e5863fffff 8 0 0
sf::st_drop_geometry(s_coarser)[, c("h3_address","res","solution_1","solution_1_deduplicated")]
#> h3_address res solution_1 solution_1_deduplicated
#> 1 87182e586ffffff 7 1 1
#> 2 87182e584ffffff 7 0 0
#> 3 87182e5b3ffffff 7 0 0
#> 4 88182e5861fffff 8 1 0
#> 5 88182e5863fffff 8 0 0
# Compute PU scores
# This function calculates a continuous importance score at each planning unit’s native resolution to create an importance ranking that can be used in the "scenario creation" (i.e., by the compute_selection_by_resolution and compute_selection_by_strata functions).
s_scored <- compute_pu_scores(
s = s,
feature_mat = feature_mat,
alpha_freq = 0.5)
sf::st_drop_geometry(s_scored)[, c("h3_address","res","ensemble_score")]
#> h3_address res ensemble_score
#> 1 87182e586ffffff 7 0.6220529
#> 2 87182e584ffffff 7 0.1250000
#> 3 87182e5b3ffffff 7 0.2470529
#> 4 88182e5861fffff 8 0.7470529
#> 5 88182e5863fffff 8 0.4970529
# Cross-scale scores:
# extends the scoring concept across resolutions and for each planning unit, the function computes scores at its "native" resolution and propagates them through the hierarchical structure. As a result, every planning unit with ancestors/descendants receives its native importance score and additional resolution-specific scores propagated through the hierarchy (one column per resolution).
# These scores allow users to compare how different scales prioritize the same location (according to the features defined at that resolution).
s_xscored <- compute_pu_scores_crossscale(
s = s,
feature_mat = feature_mat,
maps = maps,
cs_idx = cs_idx,
agg_mode = "mean",
res_col = "res"
)
sf::st_drop_geometry(s_xscored)[, c("h3_address","res","score_from_r7","score_from_r8")]
#> h3_address res score_from_r7 score_from_r8
#> 1 87182e586ffffff 7 0.75 0.5
#> 2 87182e584ffffff 7 0.25 0.0
#> 3 87182e5b3ffffff 7 0.00 0.0
#> 4 88182e5861fffff 8 0.75 1.0
#> 5 88182e5863fffff 8 0.75 0.0
# Post-hoc scenarios
# These functions rank planning units by importance and select the top units until user-defined area targets are met. Targets may be defined:
# a) By resolution (e.g. 30% of total area at each resolution, each resolution can have its own target)
s_res_scen <- compute_selection_by_resolution(
s = s_scored,
cs_idx = cs_idx,
target = 0.30,
score_col = "ensemble_score",
area_col = "area_km2",
res_col = "res",
out_col = "selected_by_resolution",
target_by_res = c("7" = 0.30, "8" = 0.50)
)
sf::st_drop_geometry(s_res_scen)[, c("h3_address","res","area_km2","ensemble_score","selected_by_resolution")]
#> h3_address res area_km2 ensemble_score selected_by_resolution
#> 1 87182e586ffffff 7 14 0.6220529 0
#> 2 87182e584ffffff 7 14 0.1250000 0
#> 3 87182e5b3ffffff 7 14 0.2470529 1
#> 4 88182e5861fffff 8 2 0.7470529 1
#> 5 88182e5863fffff 8 2 0.4970529 0
# b) By strata (simulate governance constraints, e.g. different targets per administrative unit such as countries, counties, etc)
strata_masks <- list(
A = (sf::st_drop_geometry(s_res_scen)$admin == "A"),
B = (sf::st_drop_geometry(s_res_scen)$admin == "B"),
C = (sf::st_drop_geometry(s_res_scen)$admin == "C"))
s_strat_scen <- compute_selection_by_strata(
s = s_scored,
cs_idx = cs_idx,
strata_masks = strata_masks,
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.14, B = 0.50, C = 0.50)
)
sf::st_drop_geometry(s_strat_scen)[, c("admin","h3_address","res","area_km2","ensemble_score","selected_by_admin")]
#> admin h3_address res area_km2 ensemble_score selected_by_admin
#> 1 A 87182e586ffffff 7 14 0.6220529 0
#> 2 B 87182e584ffffff 7 14 0.1250000 0
#> 3 B 87182e5b3ffffff 7 14 0.2470529 1
#> 4 C 88182e5861fffff 8 2 0.7470529 1
#> 5 C 88182e5863fffff 8 2 0.4970529 0
# Coverage summaries
# These functions account explicitly for hierarchical overlap, which is important when reporting protected areas in multiscale systems where the same location may appear simultaneously as a coarse-resolution cell and multiple finer-resolution descendants (i.e., there is no "deduplication" step). They credit coverage once per location by propagating selections across the parent–child structure before computing area totals. This ensures that reported coverage reflects the actual spatial footprint rather than raw counts of selected planning units.
summarize_coverage_by_resolution(
s = s_res_scen,
flag_col = "selected_by_resolution",
cs_idx = cs_idx
)
#> # A tibble: 2 × 4
#> res area_total_km2 area_selected_km2 coverage_prop
#> <int> <dbl> <dbl> <dbl>
#> 1 7 42 16 0.381
#> 2 8 4 2 0.5
summarize_coverage_by_strata(
s = s_strat_scen,
flag_col = "selected_by_admin",
cs_idx = cs_idx,
strata_col = "admin"
)
#> # A tibble: 3 × 5
#> strata native_res area_total_native_km2 area_selected_km2 coverage_prop
#> <chr> <int> <dbl> <dbl> <dbl>
#> 1 B 7 28 14 0.5
#> 2 C 8 4 2 0.5
#> 3 A 7 14 2 0.143
# Feature representation in a multiscale context
# Features may be defined at only one "native" resolution, but protection can occur via ancestors/descendants. These functions propagate selection through the H3 hierarchy so that representation is credited even when selection occurs in a non-native resolution, and without double-counting overlaps (see "Why multiscale feature evaluation is needed").
eval_geom_feature_coverage(
s = s_strat_scen,
flag_col = "selected_by_admin",
feature_cols = c("r7_f1", "r8_f1"),
res_col = "res",
targets = c(r7_f1 = 0.30, r8_f1 = 0.30)
)
#> # A tibble: 2 × 8
#> feature native_res total held rel_held target abs_shortfall rel_shortfall
#> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 r7_f1 7 1.2 0.143 0.119 0.3 0.217 0.181
#> 2 r8_f1 8 0.8 0.600 0.750 0.3 0 0
# NOTE: If your features are rasters (not PU-level columns), use eval_exact_raster_coverage()
# to compute area-integrated coverage directly from the raster surface.Planning unit object (sf), typically with at least:
| Column | Required | Description |
|---|---|---|
h3_address |
yes | H3 cell ID |
res |
yes | H3 resolution (integer) |
geometry |
yes | sf geometry |
area_km2 |
recommended | Used for coverage summaries |
cost |
optional | Used in optimization |
admin |
optional | Strata for area-constrained scenarios |
locked_in |
optional | Used in scoring |
Feature inputs may be:
MultiscaleSCP assumes:
The admin column represents user-defined spatial strata
(e.g., countries, counties, habitat types, management zones, or
administrative jurisdictions). In multiscale conservation planning,
strata are important because conservation commitments often operate at
political or management levels rather than purely ecological ones. For
example, a global target (e.g., 30% of Europe) may coexist with
national-level commitments (e.g., 30% per country), and local
constraints. The post-hoc functions: -
compute_selection_by_strata() (construct a selection that
meets per-stratum area targets), and -
summarize_coverage_by_strata() (report overlap-aware
coverage per stratum) are designed to mimic and evaluate these
governance constraints.
In multiscale conservation problems, ecological features (as well as costs and other attributes) will be naturally defined at specific spatial resolutions (e.g., fine-scale habitat distributions or coarse-scale species range models). As a result, conservation actions selected at different resolutions must still contribute to the representation of these features. Standard representation calculations in prioritizr (eval_feature_representation_summary()) assume that feature values are evaluated only in the planning units where they are defined. In multiscale systems, however, features are often loaded at a specific resolution while conservation decisions may occur at multiple resolutions simultaneously. For example, a feature defined at resolution r may be represented by selecting either that cell directly, one of its ancestors (coarser cells), or one of its descendants (finer cells). Standard evaluation functions do not account for these cross-scale relationships and therefore would underestimate representation unless the feature were replicated across all resolutions. The evaluation functions in MultiscaleSCP address this issue by propagating selections through the H3 hierarchy so that representation is correctly credited to parent and child cells across resolutions.
This vignette demonstrated a full multiscale workflow on a synthetic example dataset, including: