Expand description
Multi-resolution landscape connectivity on in-memory raster arrays.
This is the numerical engine of the eco-connectivity
Python package. For every unmasked focal cell of a raster tile, connectivity computes
- habitat or PARC connectedness, when no transgrids are given, or
- the Bioclimatic Ecosystem Resilience Index (BERI), when transgrids are given,
and pixel_coverage computes the fraction of each pixel covered by polygons, e.g. to
build protected-area proportion rasters. The crate has no Python, GDAL or file I/O
dependencies: callers read rasters and pass plain ndarray arrays.
§Method
Each level L of Config::levels aggregates L x L base pixels into overview cells.
Around a focal pixel, every level contributes the overview cells of an annulus (see
Config::window_size), so detail is kept near the cell and distant areas are coarse.
These cells form a graph with queen-neighbour edges within a level and links from the rim
of each level to the next coarser one. Two integer Dijkstra searches from the focal cell
give, for every reached cell v, the condition-adjusted path length a_v and the intact
path length g_v in kilometres: an edge into v of length d costs w_v * d and d,
with w_v = (1 - max_cost) * t_v + max_cost for traversal value t_v, rounded to integer
units of 10 m for routing. With q_v the condition, A_v the habitat area weight and
p_v the protected fraction (1 unless Config::pa_to_pa), connectedness for one lambda
is
sum_v exp(-a_v / lambda) q_v p_v A_v / sum_v exp(-g_v / lambda) A_v.
BERI uses Gaussian kernels exp(-(x / lambda)^2 / 283.465) instead, without p_v, and
weights each scenario’s numerator by that scenario’s climate similarity and the
denominator by the current climate’s. A similarity is exp(-L1 distance) between the
focal cell’s current-climate bands and the cell’s scenario bands, and scenario
numerators are aggregated as 0.5 * (mean + min) over all scenarios, current climate
included. The output is the mean over Config::lambdas.
§Inputs
All 2-D arrays are (rows, cols) grids of one tile with the same shape, indexed
[[row, col]] from the tile’s top-left pixel; NaN marks missing values. Transgrids add a
trailing band axis. See Inputs for each array.
Coordinates and units. Affine maps (row, col) pixel centres to world coordinates.
Projected grids (is_geo = false) must be in metres; geographic grids (is_geo = true)
are longitude/latitude in degrees, with great-circle distances on a 6371 km sphere and
habitat area weighted by cos(latitude). Distances and lambdas are in kilometres.
Transform tuples follow the rasterio/affine order (a, b, c, d, e, f), i.e.
x = a * col + b * row + c and y = d * col + e * row + f, not GDAL’s geotransform order.
Tiles and offsets. Config::offsets is the global (row0, col0) of the tile’s
top-left pixel in the full raster. Overview cells are anchored to the full raster’s grid:
level-L cell (R, C) covers global base rows R * L .. (R + 1) * L and columns
C * L .. (C + 1) * L, and a tile’s level-L grid covers every such cell that intersects
the tile, averaging only the tile’s pixels in cells cut by its edge. The level-L
transform must map that grid (level_transforms builds it from the full raster’s base
transform), so tiles evaluated separately with enough context around their cores join
without seams.
§Caller responsibilities
The core aggregates overviews from base-resolution arrays: averages of condition, traversal and transgrids, habitat area weights (valid-cell counts, or cos-latitude areas) and protected fractions. Callers, such as the Python package, remain responsible for:
- reading rasters, detecting geographic CRSs, converting nodata to NaN and scaling condition, protected-area (PA) and resistance values to [0, 1];
- planning tiles, their context padding and read windows, and passing matching offsets and level transforms;
- building the mask from polygons, nodata and, for PARC, non-protected cells;
- PARC preparation:
condition = max(condition, pa)andmask = !(pa > 0 && condition is finite); - resistance surfaces:
traversal = clip(1 - resistance, 0, 1), falling back to condition where resistance is missing and NaN where condition is missing; - BERI inputs: current-climate transgrid first, then future scenarios, bands last;
- normalising levels (add 1, round to powers of two, require a continuous sequence) and
outer_window >= window_size, which the Python wrappers check before calling; - post-processing: Connected Habitat for habitat (not PARC) connectedness with Python’s
option1 (connectedness), 2 (condition x connectedness) or 3 (the default, sqrt(condition x connectedness)); cropping tiles back to their polygons; writing rasters.
§Threads and memory
Work runs on the current Rayon thread pool: the global pool, or
the pool of an enclosing ThreadPool::install, which is how the Python binding applies
n_threads. Outputs are bit-identical for any thread count on the same build and
platform. Besides the overviews and the output, a call keeps a shared cache of window
templates and distance tables of at most 128 MiB and one scratch buffer per pool thread;
focal cells reuse those buffers.
§Example
use connectivity_core::ndarray::Array2;
use connectivity_core::{connectivity, level_transforms, Config, Inputs};
// A 24 x 24 tile at the top-left of a projected raster with 100 m pixels
let condition = Array2::from_shape_fn((24, 24), |(i, j)| ((i * 7 + j * 3) % 10) as f32 / 10.0);
let mask = condition.mapv(f32::is_nan);
let levels = [1, 2, 4];
let base = (100.0, 0.0, 500_000.0, 0.0, -100.0, 6_000_000.0);
let transforms = level_transforms(base, &levels, (0, 0));
let inputs = Inputs {
condition: &condition,
mask: &mask,
traversal: None,
transgrids: None,
transforms: &transforms,
};
let config = Config { levels: &levels, outer_window: 5, ..Config::default() };
let out = connectivity(&inputs, &config).unwrap();
assert_eq!(out.dim(), (24, 24));
assert!(out.iter().all(|v| (0.0..=1.0).contains(v)));Re-exports§
pub use ndarray;
Structs§
- Affine
- Affine transform from (row, col) pixel indices to world (x, y) coordinates.
- Config
- Settings of one
connectivitycall.Config::default()gives the Python defaults, with level 1 included as the core requires. - Inputs
- Borrowed arrays of one tile for
connectivity. 2-D arrays share the(rows, cols)shape ofcondition.
Enums§
- Window
Mode - Shape of the annulus that each level contributes around a focal cell. Cells cut by the annulus boundary count with their covered area fraction in both modes.
Functions§
- connectivity
- Connectedness, or BERI when
Inputs::transgridsis non-empty, of every unmasked cell. - level_
transforms - Level transforms of a tile, keyed by level, anchored like
rastio.read_raster:base * Affine.scale(level) * Affine.translation(col0 // level, row0 // level), with theaffinepackage’s operation order. - pixel_
coverage - Compute fractional pixel coverage for a list of polygons.
Type Aliases§
- Polygon
- One polygon as
(exterior, holes): rings of (x, y) world-coordinate vertices, in the coordinates of the grid’s transform. Rings may be open or closed and of either orientation.