Skip to main content

Crate connectivity_core

Crate connectivity_core 

Source
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) and mask = !(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 option 1 (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 connectivity call. 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 of condition.

Enums§

WindowMode
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::transgrids is 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 the affine package’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.