diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml
new file mode 100644
index 0000000..7367db5
--- /dev/null
+++ b/.github/workflows/core.yml
@@ -0,0 +1,65 @@
+name: Rust core
+
+# The connectivity-core crate on its own: no Python, its minimum supported Rust version, and a
+# crates.io packaging dry run.
+# Core releases are tagged connectivity-core-vX.Y.Z. Those tags never reach the PyPI
+# workflow (publish.yml, vX.Y.Z tags), and this workflow publishes nothing: after a core tag
+# passes, the crate is published manually with `cargo publish -p connectivity-core`.
+on:
+ push:
+ branches:
+ - main
+ tags:
+ - "connectivity-core-v*"
+ pull_request:
+ workflow_dispatch:
+
+permissions:
+ contents: read
+
+jobs:
+ core:
+ name: Core without Python
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v7
+ - name: Drop local CPU tuning
+ run: rm -f .cargo/config.toml
+ - name: Check the core does not depend on Python
+ run: |
+ if cargo tree -p connectivity-core -e normal,build,dev --locked | grep -Ei 'pyo3|numpy'; then
+ echo "connectivity-core must not depend on Python" >&2
+ exit 1
+ fi
+ # Unit, API and doc tests; core tags only run this workflow
+ - name: Test the core
+ run: cargo test -p connectivity-core --locked
+ - name: Run the synthetic example
+ run: cargo run -p connectivity-core --release --example synthetic --locked
+ - name: Check the tag matches the crate version
+ if: startsWith(github.ref, 'refs/tags/connectivity-core-v')
+ env:
+ TAG: ${{ github.ref_name }}
+ run: |
+ version=$(cargo metadata --no-deps --format-version 1 --locked \
+ | jq -r '.packages[] | select(.name == "connectivity-core") | .version')
+ test "$TAG" = "connectivity-core-v$version"
+ - name: Validate crates.io packaging
+ run: cargo publish --dry-run -p connectivity-core --locked
+
+ msrv:
+ name: Core on its minimum Rust version
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v7
+ - name: Drop local CPU tuning
+ run: rm -f .cargo/config.toml
+ # The version comes from rust-version in crates/connectivity-core/Cargo.toml
+ - name: Install the minimum Rust version
+ run: |
+ msrv=$(cargo metadata --no-deps --format-version 1 --locked \
+ | jq -r '.packages[] | select(.name == "connectivity-core") | .rust_version')
+ rustup toolchain install "$msrv" --profile minimal
+ echo "MSRV=$msrv" >> "$GITHUB_ENV"
+ - name: Build the core, its tests and example
+ run: cargo +"$MSRV" check -p connectivity-core --all-targets --locked
diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
index 1f846cc..c6068ef 100644
--- a/.github/workflows/publish.yml
+++ b/.github/workflows/publish.yml
@@ -5,7 +5,8 @@ on:
branches:
- main
tags:
- - "v*"
+ # Python releases (eco-connectivity); connectivity-core uses connectivity-core-v* tags
+ - "v[0-9]*"
pull_request:
workflow_dispatch:
@@ -60,8 +61,8 @@ jobs:
- uses: actions/setup-python@v7
with:
python-version: "3.13"
- - name: Run unit tests
- run: cargo test --locked
+ - name: Run workspace tests
+ run: cargo test --workspace --locked
publish:
name: Publish to PyPI
@@ -83,6 +84,11 @@ jobs:
pattern: cibw-*
path: dist
merge-multiple: true
+ - name: Check the tag matches the package version
+ # Only a vX.Y.Z tag naming the eco-connectivity version may publish
+ env:
+ TAG: ${{ github.ref_name }}
+ run: ls "dist/eco_connectivity-${TAG#v}.tar.gz" dist/eco_connectivity-"${TAG#v}"-*.whl
- uses: actions/attest-build-provenance@v4
with:
subject-path: dist/*
diff --git a/Cargo.lock b/Cargo.lock
index 2c592f9..2957c0e 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -14,6 +14,18 @@ version = "1.4.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26"
+[[package]]
+name = "connectivity-core"
+version = "0.1.0"
+dependencies = [
+ "anyhow",
+ "ndarray",
+ "num_cpus",
+ "pathfinding",
+ "rayon",
+ "rustc-hash",
+]
+
[[package]]
name = "crossbeam-deque"
version = "0.8.6"
@@ -342,16 +354,15 @@ name = "rust_conn"
version = "2.4.0"
dependencies = [
"anyhow",
+ "connectivity-core",
"ndarray",
"num_cpus",
"numpy",
"ordered-float",
- "pathfinding",
"pyo3",
"pyo3-build-config",
"pyo3-built",
"rayon",
- "rustc-hash",
]
[[package]]
diff --git a/Cargo.toml b/Cargo.toml
index 7f786c1..382436a 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -3,6 +3,8 @@ name = "rust_conn"
version = "2.4.0"
edition = "2021"
license = "MIT"
+# Released as the eco-connectivity Python wheels; only connectivity-core goes to crates.io
+publish = false
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
[lib]
@@ -10,23 +12,31 @@ name = "rust_conn"
crate-type = ["cdylib"]
[dependencies]
+connectivity-core = { path = "crates/connectivity-core", version = "0.1.0" }
pyo3 = {version = "0.28", default-features = false, features = ["macros", "abi3-py38"]}
pyo3-built = "0.5"
numpy = "0.28"
-ndarray = { version = "0.16", features = ["rayon"] }
+ndarray = { workspace = true }
ordered-float = "4.2.0"
-rayon = "1.11"
-num_cpus = "1.16"
-anyhow = "1.0.100"
-rustc-hash = "2.1"
-
-[dev-dependencies]
-# Reference Dijkstra for routing tests
-pathfinding = "4.14.0"
+rayon = { workspace = true }
+num_cpus = { workspace = true }
+anyhow = { workspace = true }
[build-dependencies]
pyo3-build-config = { version = "0.28", features = ["resolve-config"] }
+[workspace]
+members = ["crates/connectivity-core"]
+# Plain `cargo test` and `cargo build` at the root also cover the core
+default-members = [".", "crates/connectivity-core"]
+
+# Versions shared by the extension and the core
+[workspace.dependencies]
+anyhow = "1.0.100"
+ndarray = { version = "0.16", features = ["rayon"] }
+num_cpus = "1.16"
+rayon = "1.11"
+
[profile.release]
# Stripped builds fail to load on recent macOS (dyld: "mis-aligned LINKEDIT string pool")
strip = false
diff --git a/NEWS.md b/NEWS.md
index 6d7d0b6..9b3a5fa 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,7 @@
- Batched BERI similarity calculations across all scenarios of each cell, preserving similarity values.
- Fixed each node's outgoing edge order: same-level queen neighbours in the existing order, then fringe links in `get_edge_neighbours` order. Shortest-path costs are unchanged, but integer-cost ties can select different parents and adjusted path sums, and indicator accumulation follows a different destination order. On the bundled example 2-16% of cells change, by up to about 1e-4 (2.4e-4 relative); grids where most edges round to zero integer cost can change more. Repeated runs and different thread counts remain bit-identical on the same build and platform.
- Improved Rust `connectivity()` validation: missing transforms are checked before processing, and skipped powers of two produce a clear panic message instead of an index error. Both previously surfaced when a focal cell first reached the affected level. Python `connectedness()` and `beri()` already validate both.
+- Moved the Rust numerical engine (overviews, window templates, graph building, routing, indicator metrics and pixel coverage) into a separate crate, `connectivity-core` 0.1.0 in `crates/connectivity-core`, with a documented array-based Rust API and no Python, GDAL or file I/O dependencies. The `rust_conn` extension is now a thin PyO3 wrapper around it: Python function signatures, outputs, exceptions and threading are unchanged. The core is versioned and tagged separately (`connectivity-core-vX.Y.Z`); only `vX.Y.Z` tags publish the Python package to PyPI.
# Version 2.3.0
## Changed
diff --git a/README.md b/README.md
index 6fc34ba..c7cf6b2 100644
--- a/README.md
+++ b/README.md
@@ -11,6 +11,7 @@
- [Installation ](#installation-)
- [Install from PyPI ](#install-from-pypi-)
- [Install from Source ](#install-from-source-)
+ - [Rust Core Library ](#rust-core-library-)
- [Connectivity Analysis ](#connectivity-analysis-)
- [Connected Habitat (Connectedness) ](#connected-habitat-connectedness-)
- [Optional resistance surface](#optional-resistance-surface)
@@ -45,6 +46,10 @@ cd connectivity
pip install .
```
+### Rust Core Library
+
+The numerical engine is a separate Rust crate, [`connectivity-core`](https://github.com/csiro/connectivity/tree/main/crates/connectivity-core), which works on in-memory arrays without Python, GDAL or file I/O. Rust programs can use it directly; see its [README](https://github.com/csiro/connectivity/blob/main/crates/connectivity-core/README.md) and crate documentation. `cargo test --workspace` runs every Rust test, and `cargo test -p connectivity-core` needs no Python.
+
## Connectivity Analysis
Resolve the bundled example paths once so they can be reused throughout the
analysis:
diff --git a/crates/connectivity-core/Cargo.toml b/crates/connectivity-core/Cargo.toml
new file mode 100644
index 0000000..cf74901
--- /dev/null
+++ b/crates/connectivity-core/Cargo.toml
@@ -0,0 +1,25 @@
+[package]
+name = "connectivity-core"
+version = "0.1.0"
+edition = "2021"
+# f64::next_up and next_down (distance tables)
+rust-version = "1.86"
+license = "MIT"
+description = "Multi-resolution landscape connectivity, PARC connectedness and BERI on in-memory raster arrays"
+authors = ["Roozbeh Valavi "]
+repository = "https://github.com/csiro/connectivity"
+readme = "README.md"
+keywords = ["connectivity", "landscape-ecology", "biodiversity", "raster", "least-cost-path"]
+categories = ["science::geo", "algorithms"]
+
+[dependencies]
+anyhow = { workspace = true }
+ndarray = { workspace = true }
+rayon = { workspace = true }
+
+[dev-dependencies]
+num_cpus = { workspace = true }
+# Reference Dijkstra for routing tests
+pathfinding = "4.14.0"
+# Hash maps of the frozen reference implementation
+rustc-hash = "2.1"
diff --git a/crates/connectivity-core/LICENSE b/crates/connectivity-core/LICENSE
new file mode 100644
index 0000000..bdf8c54
--- /dev/null
+++ b/crates/connectivity-core/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2026 Commonwealth Scientific and Industrial Research Organisation (CSIRO)
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/crates/connectivity-core/README.md b/crates/connectivity-core/README.md
new file mode 100644
index 0000000..5ff4db5
--- /dev/null
+++ b/crates/connectivity-core/README.md
@@ -0,0 +1,63 @@
+# connectivity-core
+
+Multi-resolution landscape connectivity on in-memory raster arrays: habitat and PARC
+connectedness, the Bioclimatic Ecosystem Resilience Index (BERI), and fractional pixel
+coverage of polygons. This crate is the numerical engine of the
+[`eco-connectivity`](https://pypi.org/project/eco-connectivity/) Python package; see
+[Valavi et al. (2026)](https://doi.org/10.32942/X2S68V) for the method.
+
+It has no Python, GDAL or file I/O dependencies. Callers read rasters, pass `ndarray` arrays
+of one tile, and handle file formats, tiling and post-processing such as Connected Habitat.
+
+```toml
+[dependencies]
+connectivity-core = "0.1"
+```
+
+It requires Rust 1.86 or later.
+
+```rust
+use connectivity_core::ndarray::Array2;
+use connectivity_core::{connectivity, level_transforms, Config, Inputs};
+
+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];
+// rasterio/affine (a, b, c, d, e, f) order; projected grids in metres
+let transforms = level_transforms((100.0, 0.0, 500_000.0, 0.0, -100.0, 6_000_000.0), &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 connectedness = connectivity(&inputs, &config)?;
+```
+
+## Public API
+
+- `connectivity(&Inputs, &Config)`: connectedness, or BERI when `Inputs::transgrids` holds the
+ current climate followed by future scenarios (each `(rows, cols, bands)`). Returns an array
+ with the condition array's shape, NaN where the mask is `true`.
+- `Inputs`: borrowed condition, mask, optional traversal and transgrids, and level transforms.
+- `Config`: levels, lambdas (km), `is_geo`, `max_cost`, window sizes, global tile offsets,
+ window mode and PARC `pa_to_pa` gating. `Config::default()` matches the Python defaults.
+- `pixel_coverage(&[Polygon], &Affine, nrows, ncols)`: fraction of each pixel covered by
+ polygons with holes.
+- `Affine` and `level_transforms`: (row, col) to world transforms and the per-level transforms
+ of a tile, anchored to the full raster's grid.
+- `WindowMode`: circular or square annuli.
+
+The crate documentation covers array shapes, mask semantics, BERI scenario and band order,
+coordinate units, tile offsets, the preparation and post-processing left to callers, and the
+error and panic conditions. Work runs on the current Rayon thread pool; results are
+bit-identical for any number of threads on the same build and platform.
+
+`cargo run --release -p connectivity-core --example synthetic` runs every indicator on
+synthetic arrays.
+
+## Versioning and releases
+
+The crate is versioned independently of the Python package (`eco-connectivity` 2.4.0 uses
+`connectivity-core` 0.1.0). Core releases are tagged `connectivity-core-vX.Y.Z`; Python
+releases keep `vX.Y.Z` tags, and only those publish to PyPI.
+
+## License
+
+MIT
diff --git a/crates/connectivity-core/examples/synthetic.rs b/crates/connectivity-core/examples/synthetic.rs
new file mode 100644
index 0000000..94fbfe1
--- /dev/null
+++ b/crates/connectivity-core/examples/synthetic.rs
@@ -0,0 +1,69 @@
+//! Connectedness, PARC, BERI and pixel coverage on synthetic arrays, with no file I/O.
+//!
+//! `cargo run --release -p connectivity-core --example synthetic`
+use connectivity_core::ndarray::{Array2, Array3, Zip};
+use connectivity_core::{connectivity, level_transforms, pixel_coverage, Config, Inputs};
+use std::error::Error;
+
+fn summary(name: &str, values: &Array2) {
+ let finite: Vec = values.iter().copied().filter(|v| v.is_finite()).collect();
+ let mean = finite.iter().map(|&v| v as f64).sum::() / finite.len().max(1) as f64;
+ let max = finite.iter().copied().fold(f32::NEG_INFINITY, f32::max);
+ println!("{name:24} {} cells evaluated, mean {mean:.4}, max {max:.4}", finite.len());
+}
+
+fn main() -> Result<(), Box> {
+ // A 96 x 128 tile whose top-left pixel is row 40, column 24 of a projected raster with
+ // 100 m pixels; overview cells and transforms are anchored to the full raster's grid
+ let (rows, cols, offsets) = (96, 128, (40, 24));
+ let base = (100.0, 0.0, 400_000.0, 0.0, -100.0, 6_100_000.0);
+ let levels = [1, 2, 4, 8];
+ let transforms = level_transforms(base, &levels, offsets);
+
+ // Habitat patches in [0, 1] around a lake of missing data; masked cells are skipped
+ let condition = Array2::from_shape_fn((rows, cols), |(i, j)| {
+ let (y, x) = (i as f32 - 48.0, j as f32 - 64.0);
+ if x * x + y * y < 120.0 {
+ f32::NAN
+ } else {
+ (0.5 + 0.5 * (i as f32 / 9.0).sin() * (j as f32 / 13.0).cos()).clamp(0.0, 1.0)
+ }
+ });
+ let mask = condition.mapv(f32::is_nan);
+ let inputs = Inputs { condition: &condition, mask: &mask, traversal: None, transgrids: None, transforms: &transforms };
+ let config = Config { levels: &levels, lambdas: &[2.0, 20.0], outer_window: 5, offsets, ..Config::default() };
+
+ let connectedness = connectivity(&inputs, &config)?;
+ summary("connectedness", &connectedness);
+ // Connected Habitat is caller post-processing; Python's default is option 3
+ let mut habitat = connectedness.clone();
+ Zip::from(&mut habitat).and(&condition).for_each(|h, &c| *h = (c * *h).sqrt());
+ summary("connected habitat", &habitat);
+
+ // A road two rows wide that is hard to cross (traversal = 1 - resistance) on otherwise
+ // condition-based traversal
+ let traversal = Array2::from_shape_fn((rows, cols), |(i, j)| if i == 70 || i == 71 { 0.05 } else { condition[[i, j]] });
+ let with_road = connectivity(&Inputs { traversal: Some(&traversal), ..inputs }, &config)?;
+ summary("connectedness, road", &with_road);
+
+ // PARC: protected-area proportions from a reserve polygon with a hole, in world coordinates
+ let tile = transforms[&1];
+ let (x0, y0) = (tile.x_origin, tile.y_origin);
+ let square = |a: f64, b: f64, c: f64, d: f64| vec![(x0 + a, y0 - b), (x0 + c, y0 - b), (x0 + c, y0 - d), (x0 + a, y0 - d)];
+ let reserve = vec![(square(500.0, 300.0, 4_050.0, 3_520.0), vec![square(1_500.0, 1_000.0, 2_500.0, 2_000.0)])];
+ let pa = pixel_coverage(&reserve, &tile, rows, cols);
+ // As numpy.maximum in Python: missing condition stays missing
+ let parc_condition = Zip::from(&condition).and(&pa).map_collect(|&c, &p| if c.is_nan() { c } else { c.max(p) });
+ let parc_mask = Zip::from(&parc_condition).and(&pa).map_collect(|&c, &p| !(p > 0.0 && c.is_finite()));
+ let parc_inputs = Inputs { condition: &parc_condition, mask: &parc_mask, ..inputs };
+ summary("PARC", &connectivity(&parc_inputs, &config)?);
+ summary("PARC, pa_to_pa", &connectivity(&parc_inputs, &Config { pa_to_pa: true, ..config })?);
+
+ // BERI: the current climate first, then two warmer scenarios, with three bands each
+ let climate = |shift: f32| {
+ Array3::from_shape_fn((rows, cols, 3), |(i, j, b)| (i as f32 / 40.0 + j as f32 / 60.0) * (b + 1) as f32 * 0.1 + shift)
+ };
+ let scenarios = [climate(0.0), climate(0.05), climate(0.15)];
+ summary("BERI", &connectivity(&Inputs { transgrids: Some(&scenarios), ..inputs }, &config)?);
+ Ok(())
+}
diff --git a/crates/connectivity-core/src/affine.rs b/crates/connectivity-core/src/affine.rs
new file mode 100644
index 0000000..6f07262
--- /dev/null
+++ b/crates/connectivity-core/src/affine.rs
@@ -0,0 +1,80 @@
+use std::collections::HashMap;
+use std::convert::From;
+
+/// Affine transform from (row, col) pixel indices to world (x, y) coordinates.
+///
+/// The fields follow the rasterio/`affine` tuple order `(a, b, c, d, e, f)`, i.e.
+/// `tuple(transform)[:6]` in Python, not GDAL's `GetGeoTransform` order:
+/// `x = a * col + b * row + c` and `y = d * col + e * row + f`.
+/// Projected coordinates are in metres; geographic coordinates are longitude (x) and
+/// latitude (y) in degrees.
+#[derive(Debug, Clone, Copy)]
+pub struct Affine {
+ pub x_scale: f64, // pixel width
+ pub x_skew: f64, // rotation/skew term (usually 0)
+ pub x_origin: f64, // top-left corner X
+ pub y_skew: f64, // rotation/skew term (usually 0)
+ pub y_scale: f64, // pixel height (negative if north-up)
+ pub y_origin: f64, // top-left corner Y
+}
+
+impl Affine {
+ /// Convert (row, col) to (x, y) at pixel centre by default.
+ pub fn xy(&self, row: i32, col: i32) -> (f64, f64) {
+ let col_f = (col as f64) + 0.5;
+ let row_f = (row as f64) + 0.5;
+
+ let x = self.x_scale * col_f + self.x_skew * row_f + self.x_origin;
+ let y = self.y_skew * col_f + self.y_scale * row_f + self.y_origin;
+ (x, y)
+ }
+
+ /// Convert world (x, y) to continuous (col, row) pixel coordinates.
+ pub fn col_row(&self, x: f64, y: f64) -> (f64, f64) {
+ let det = self.x_scale * self.y_scale - self.x_skew * self.y_skew;
+ let dx = x - self.x_origin;
+ let dy = y - self.y_origin;
+ let col = (self.y_scale * dx - self.x_skew * dy) / det;
+ let row = (-self.y_skew * dx + self.x_scale * dy) / det;
+ (col, row)
+ }
+}
+
+impl From<(f64, f64, f64, f64, f64, f64)> for Affine {
+ /// From an `(a, b, c, d, e, f)` tuple in rasterio/`affine` order.
+ fn from(t: (f64, f64, f64, f64, f64, f64)) -> Self {
+ Self {
+ x_scale: t.0,
+ x_skew: t.1,
+ x_origin: t.2,
+ y_skew: t.3,
+ y_scale: t.4,
+ y_origin: t.5,
+ }
+ }
+}
+
+/// 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.
+///
+/// `base` is the full raster's base-level transform and `offsets` the tile's `(row0, col0)`
+/// in base pixels, so level `L` maps the tile's overview grid, whose cell (0, 0) is global
+/// overview cell `(row0 / L, col0 / L)`, and level 1 maps the tile itself.
+pub fn level_transforms(
+ base: impl Into,
+ levels: &[usize],
+ offsets: (usize, usize),
+) -> HashMap {
+ let Affine { x_scale: a, x_skew: b, x_origin: c, y_skew: d, y_scale: e, y_origin: f } = base.into();
+ levels
+ .iter()
+ .map(|&level| {
+ let l = level as f64;
+ let (sa, sb, sd, se) = (a * l + b * 0.0, a * 0.0 + b * l, d * l + e * 0.0, d * 0.0 + e * l);
+ let (tc, tr) = ((offsets.1 / level) as f64, (offsets.0 / level) as f64);
+ let t = (sa, sb, sa * tc + sb * tr + c, sd, se, sd * tc + se * tr + f);
+ (level as i32, Affine::from(t))
+ })
+ .collect()
+}
diff --git a/src/builder.rs b/crates/connectivity-core/src/builder.rs
similarity index 96%
rename from src/builder.rs
rename to crates/connectivity-core/src/builder.rs
index e124111..791abb6 100644
--- a/src/builder.rs
+++ b/crates/connectivity-core/src/builder.rs
@@ -482,23 +482,7 @@ pub mod tests {
use crate::testkit::Fixture;
pub fn context(f: &Fixture) -> Context<'_> {
- prepare(
- &f.cond,
- f.trav.as_ref(),
- f.trans.as_deref(),
- &f.transforms,
- &f.mask,
- &f.levels,
- &f.lambdas,
- f.is_geo,
- f.max_cost,
- f.window_size,
- f.outer_window,
- f.offsets,
- f.mode,
- f.pa_to_pa,
- )
- .unwrap()
+ prepare(&f.inputs(), &f.config()).unwrap()
}
/// Legacy node id of optimized node u.
@@ -713,14 +697,7 @@ pub mod tests {
let canonical = f.canonical(None);
for threads in [1, 2, num_cpus::get().max(3)] {
let pool = rayon::ThreadPoolBuilder::new().num_threads(threads).build().unwrap();
- let out = pool.install(|| {
- crate::core::conn(
- &Some(f.cond.clone()), f.trav.as_ref(), f.trans.as_deref(), &f.transforms,
- &f.mask, &f.levels, &f.lambdas, f.is_geo, f.max_cost, f.window_size,
- f.outer_window, f.offsets, f.mode, f.pa_to_pa,
- )
- .unwrap()
- });
+ let out = pool.install(|| crate::core::conn(&f.inputs(), &f.config()).unwrap());
assert!(crate::testkit::same_bits(&out, &canonical), "seed {seed}, {threads} threads");
}
}
@@ -860,7 +837,7 @@ pub mod tests {
f.mask.indexed_iter_mut().for_each(|((i, j), m)| *m |= (i * 7 + j * 3) % 5 != 0);
let base_t = f.transforms[&1];
let t = (base_t.x_scale, base_t.x_skew, base_t.x_origin, base_t.y_skew, base_t.y_scale, base_t.y_origin);
- f.transforms = crate::testkit::level_transforms(t, &f.levels, f.offsets);
+ f.transforms = crate::level_transforms(t, &f.levels, f.offsets);
let mut ctx = context(&f);
for budget in [budget, 64 << 10, 0] {
let (tables, templates, retained, peak) = measured(&ctx, &f, budget);
@@ -876,13 +853,7 @@ pub mod tests {
(ctx.tables, ctx.templates) = (tables, templates);
compare_ctx(&f, &ctx);
}
- let stats = crate::core::conn_stats(
- &Some(f.cond.clone()), f.trav.as_ref(), f.trans.as_deref(), &f.transforms, &f.mask,
- &f.levels, &f.lambdas, f.is_geo, f.max_cost, f.window_size, f.outer_window,
- f.offsets, f.mode, f.pa_to_pa,
- )
- .unwrap()
- .1;
+ let stats = crate::core::conn_stats(&f.inputs(), &f.config()).unwrap().1;
assert!(stats.cache_peak <= crate::template::CACHE_BUDGET);
assert!(stats.cache_bytes <= stats.cache_peak && stats.cache_fixed <= stats.cache_bytes);
assert_eq!(stats.scratch.len(), rayon::current_num_threads());
@@ -908,7 +879,7 @@ pub mod tests {
f.levels = (0..=8).map(|k| 1usize << k).collect();
(f.window_size, f.outer_window) = (window_size, outer_window);
f.offsets = (5, 3);
- f.transforms = crate::testkit::level_transforms(
+ f.transforms = crate::level_transforms(
(100.0, 0.0, 400_000.0, 0.0, -100.0, 6_000_000.0),
&f.levels,
f.offsets,
@@ -1024,12 +995,7 @@ pub mod tests {
let ctx = context(&f);
let inject = table_inject(&ctx);
let injected = f.canonical(Some(&inject));
- let optimized = crate::core::conn(
- &Some(f.cond.clone()), f.trav.as_ref(), f.trans.as_deref(), &f.transforms, &f.mask,
- &f.levels, &f.lambdas, f.is_geo, f.max_cost, f.window_size, f.outer_window,
- f.offsets, f.mode, f.pa_to_pa,
- )
- .unwrap();
+ let optimized = crate::core::conn(&f.inputs(), &f.config()).unwrap();
assert!(crate::testkit::same_bits(&optimized, &injected), "fixture {n}");
let reference = f.canonical(None);
for (&a, &b) in injected.iter().zip(&reference) {
diff --git a/src/core.rs b/crates/connectivity-core/src/core.rs
similarity index 76%
rename from src/core.rs
rename to crates/connectivity-core/src/core.rs
index 9868257..d104f20 100644
--- a/src/core.rs
+++ b/crates/connectivity-core/src/core.rs
@@ -11,10 +11,9 @@ use crate::overview;
use crate::tables::{Tables, STATUSES};
use crate::template::{Templates, CACHE_BUDGET};
use crate::utils;
-use crate::window;
+use crate::{Config, Inputs};
use affine::Affine;
use overview::Resampling;
-use window::WindowMode;
/// Call-level resource statistics.
#[derive(Debug, Default, Clone)]
@@ -32,150 +31,88 @@ pub struct Stats {
pub process_s: f64, // focal-cell loop
}
-/// Core connectivity function
+/// Core connectivity function (`crate::connectivity`)
/// # Returns
/// A 2D array of connectivity values for each cell at the native resolution.
-pub fn conn(
- cond_array: &Option>,
- trav_array: Option<&Array2>,
- trans_arrays: Option<&[Array3]>,
- transform_map: &HashMap,
- mask_array: &Array2,
- levels: &Vec,
- lambdas: &[f32],
- is_geo: bool,
- max_cost: f32,
- window_size: i32,
- outer_window: i32,
- offsets: (usize, usize),
- window_mode: WindowMode,
- pa_to_pa: bool,
-) -> Result> {
- conn_stats(
- cond_array,
- trav_array,
- trans_arrays,
- transform_map,
- mask_array,
- levels,
- lambdas,
- is_geo,
- max_cost,
- window_size,
- outer_window,
- offsets,
- window_mode,
- pa_to_pa,
- )
- .map(|(out, _)| out)
+pub fn conn(inputs: &Inputs, config: &Config) -> Result> {
+ conn_stats(inputs, config).map(|(out, _)| out)
}
/// `conn` that also reports scratch memory and setup/processing time.
-pub fn conn_stats(
- cond_array: &Option>,
- trav_array: Option<&Array2>,
- trans_arrays: Option<&[Array3]>,
- transform_map: &HashMap,
- mask_array: &Array2,
- levels: &Vec,
- lambdas: &[f32],
- is_geo: bool,
- max_cost: f32,
- window_size: i32,
- outer_window: i32,
- offsets: (usize, usize),
- window_mode: WindowMode,
- pa_to_pa: bool,
-) -> Result<(Array2, Stats)> {
+pub fn conn_stats(inputs: &Inputs, config: &Config) -> Result<(Array2, Stats)> {
let start = Instant::now();
- // Check condition dictionay was not empty and run the code for level 1 (original resolution)
- if let Some(cond_base) = cond_array {
- let ctx = prepare(
- cond_base,
- trav_array,
- trans_arrays,
- transform_map,
- mask_array,
- levels,
- lambdas,
- is_geo,
- max_cost,
- window_size,
- outer_window,
- offsets,
- window_mode,
- pa_to_pa,
- )?;
+ let (cond_base, mask_array) = (inputs.condition, inputs.mask);
+ let ctx = prepare(inputs, config)?;
- // One scratch per pool thread, sized up front; a row borrows its thread's scratch, so
- // focal cells only allocate if a window outgrows the initial capacity.
- let threads = rayon::current_num_threads();
- let scratch: Vec> =
- (0..threads).map(|_| Mutex::new(Scratch::new(&ctx))).collect();
- let setup_s = start.elapsed().as_secs_f64();
+ // One scratch per pool thread, sized up front; a row borrows its thread's scratch, so
+ // focal cells only allocate if a window outgrows the initial capacity.
+ let threads = rayon::current_num_threads();
+ let scratch: Vec> =
+ (0..threads).map(|_| Mutex::new(Scratch::new(&ctx))).collect();
+ let setup_s = start.elapsed().as_secs_f64();
- // Initialize output with NaNs so masked cells keep the previous output semantics.
- let ncols = cond_base.ncols();
- let mut outarray = Array2::::from_elem(cond_base.dim(), f32::NAN);
- outarray
- .axis_iter_mut(Axis(0))
- .into_par_iter()
- .enumerate()
- .for_each(|(i, mut row_result)| {
- let slot = rayon::current_thread_index().unwrap_or(0) % threads;
- // A panicking cell poisons its scratch, but every cell resets what it uses
- let mut scratch = scratch[slot].lock().unwrap_or_else(PoisonError::into_inner);
- for j in 0..ncols {
- // Skip masked areas
- if mask_array[[i, j]] {
- continue;
- }
- row_result[j] = scratch.eval(&ctx, i, j);
+ // Initialize output with NaNs so masked cells keep the previous output semantics.
+ let ncols = cond_base.ncols();
+ let mut outarray = Array2::::from_elem(cond_base.dim(), f32::NAN);
+ outarray
+ .axis_iter_mut(Axis(0))
+ .into_par_iter()
+ .enumerate()
+ .for_each(|(i, mut row_result)| {
+ let slot = rayon::current_thread_index().unwrap_or(0) % threads;
+ // A panicking cell poisons its scratch, but every cell resets what it uses
+ let mut scratch = scratch[slot].lock().unwrap_or_else(PoisonError::into_inner);
+ for j in 0..ncols {
+ // Skip masked areas
+ if mask_array[[i, j]] {
+ continue;
}
- });
+ row_result[j] = scratch.eval(&ctx, i, j);
+ }
+ });
- let stats = Stats {
- scratch: scratch
- .iter()
- .map(|s| s.lock().unwrap_or_else(PoisonError::into_inner).bytes())
- .collect(),
- cache_bytes: ctx.templates.bytes + ctx.tables.bytes,
- cache_peak: ctx.templates.peak + ctx.tables.bytes,
- cache_fixed: ctx.templates.fixed + ctx.tables.fixed,
- table_bytes: ctx.tables.bytes,
- cached: ctx.templates.cached.to_vec(),
- linked: ctx.templates.linked.to_vec(),
- edges: scratch.iter().fold([0; STATUSES], |mut acc, s| {
- let s = s.lock().unwrap_or_else(PoisonError::into_inner);
- acc.iter_mut().zip(s.counts).for_each(|(a, c)| *a += c);
- acc
- }),
- setup_s,
- process_s: start.elapsed().as_secs_f64() - setup_s,
- };
- Ok((outarray, stats))
- } else {
- anyhow::bail!("No base level (key=1) in condition array");
- }
+ let stats = Stats {
+ scratch: scratch
+ .iter()
+ .map(|s| s.lock().unwrap_or_else(PoisonError::into_inner).bytes())
+ .collect(),
+ cache_bytes: ctx.templates.bytes + ctx.tables.bytes,
+ cache_peak: ctx.templates.peak + ctx.tables.bytes,
+ cache_fixed: ctx.templates.fixed + ctx.tables.fixed,
+ table_bytes: ctx.tables.bytes,
+ cached: ctx.templates.cached.to_vec(),
+ linked: ctx.templates.linked.to_vec(),
+ edges: scratch.iter().fold([0; STATUSES], |mut acc, s| {
+ let s = s.lock().unwrap_or_else(PoisonError::into_inner);
+ acc.iter_mut().zip(s.counts).for_each(|(a, c)| *a += c);
+ acc
+ }),
+ setup_s,
+ process_s: start.elapsed().as_secs_f64() - setup_s,
+ };
+ Ok((outarray, stats))
}
/// Overviews and settings for one call, indexed by ascending level.
-pub fn prepare<'a>(
- cond_base: &Array2,
- trav_array: Option<&Array2>,
- trans_arrays: Option<&[Array3]>,
- transform_map: &HashMap,
- mask_array: &Array2,
- levels: &Vec,
- lambdas: &'a [f32],
- is_geo: bool,
- max_cost: f32,
- window_size: i32,
- outer_window: i32,
- offsets: (usize, usize),
- window_mode: WindowMode,
- pa_to_pa: bool,
-) -> Result> {
+pub fn prepare<'a>(inputs: &Inputs, config: &Config<'a>) -> Result> {
+ let Inputs {
+ condition: cond_base,
+ mask: mask_array,
+ traversal: trav_array,
+ transgrids: trans_arrays,
+ transforms: transform_map,
+ } = *inputs;
+ let Config {
+ levels,
+ lambdas,
+ is_geo,
+ max_cost,
+ window_size,
+ outer_window,
+ offsets,
+ window_mode,
+ pa_to_pa,
+ } = *config;
// If transgrids are provided run BERI, otherwise connectedness.
let run_beri = trans_arrays
.map(|arrays| !arrays.is_empty())
@@ -380,13 +317,14 @@ pub fn prepare<'a>(
#[cfg(test)]
mod tests {
use super::*;
+ use crate::window::WindowMode;
// Repeated runs must give bit-identical output.
#[test]
fn repeats_exactly() {
- let cond = Some(Array2::from_shape_fn((32, 32), |(i, j)| {
+ let cond = Array2::from_shape_fn((32, 32), |(i, j)| {
((i * 7 + j * 13) % 10 + 1) as f32 / 10.0
- }));
+ });
let mask = Array2::from_elem((32, 32), false);
let levels = vec![1, 2, 4];
let transforms: HashMap = levels
@@ -396,13 +334,12 @@ mod tests {
(l as i32, Affine::from((s, 0.0, 0.0, 0.0, -s, 32_000.0)))
})
.collect();
- let run = || {
- conn(
- &cond, None, None, &transforms, &mask, &levels, &[2.0, 20.0], false, 2.0, 3, 5,
- (0, 0), WindowMode::Circular, false,
- )
- .unwrap()
+ let inputs = Inputs { condition: &cond, mask: &mask, traversal: None, transgrids: None, transforms: &transforms };
+ let config = Config {
+ levels: &levels, lambdas: &[2.0, 20.0], is_geo: false, max_cost: 2.0, window_size: 3,
+ outer_window: 5, offsets: (0, 0), window_mode: WindowMode::Circular, pa_to_pa: false,
};
+ let run = || conn(&inputs, &config).unwrap();
let first = run();
assert!(first.iter().any(|&v| v > 0.0));
@@ -568,23 +505,7 @@ mod tests {
(WindowMode::Square, true, true, [0.288273066, 0.334522188, 0.339534730, 0.269480884]),
] {
let f = crate::testkit::Fixture::multilevel(mode, beri, pa_to_pa);
- let output = conn(
- &Some(f.cond.clone()),
- f.trav.as_ref(),
- f.trans.as_deref(),
- &f.transforms,
- &f.mask,
- &f.levels,
- &f.lambdas,
- f.is_geo,
- f.max_cost,
- f.window_size,
- f.outer_window,
- f.offsets,
- f.mode,
- f.pa_to_pa,
- )
- .unwrap();
+ let output = conn(&f.inputs(), &f.config()).unwrap();
for (cell, expected) in cells.iter().zip(expected) {
let actual = output[*cell];
assert!(
diff --git a/src/coverage.rs b/crates/connectivity-core/src/coverage.rs
similarity index 87%
rename from src/coverage.rs
rename to crates/connectivity-core/src/coverage.rs
index d2e7d3a..ad9c5e0 100644
--- a/src/coverage.rs
+++ b/crates/connectivity-core/src/coverage.rs
@@ -2,16 +2,22 @@ use crate::affine::Affine;
use ndarray::parallel::prelude::*;
use ndarray::{s, Array2, Axis};
+/// 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.
+pub type Polygon = (Vec<(f64, f64)>, Vec>);
+
/// Compute fractional pixel coverage for a list of polygons.
///
-/// Each polygon is `(exterior, holes)` with rings as `Vec<(f64, f64)>` of (x, y)
-/// world-coordinate vertices. Polygons are assumed already-unioned on the Python side;
-/// any residual overlap is clipped to 1.0.
+/// Each pixel of the `nrows` x `ncols` grid mapped by `transform` gets the fraction of its
+/// area covered by `polygons`: exterior rings add their area and holes subtract it. Polygons
+/// are assumed already unioned by the caller (the Python wrapper unions them first); any
+/// residual overlap is clipped to 1.0. Rings with fewer than three vertices are ignored.
+/// Rows are processed in parallel on the current Rayon thread pool.
///
/// # Returns
/// A 2D array of shape (nrows, ncols) with values in [0, 1].
-pub fn pixel_coverage_array(
- polygons: &[(Vec<(f64, f64)>, Vec>)],
+pub fn pixel_coverage(
+ polygons: &[Polygon],
transform: &Affine,
nrows: usize,
ncols: usize,
diff --git a/src/distances.rs b/crates/connectivity-core/src/distances.rs
similarity index 100%
rename from src/distances.rs
rename to crates/connectivity-core/src/distances.rs
diff --git a/src/graph.rs b/crates/connectivity-core/src/graph.rs
similarity index 100%
rename from src/graph.rs
rename to crates/connectivity-core/src/graph.rs
diff --git a/src/legacy/builder.rs b/crates/connectivity-core/src/legacy/builder.rs
similarity index 100%
rename from src/legacy/builder.rs
rename to crates/connectivity-core/src/legacy/builder.rs
diff --git a/src/legacy/canonical.rs b/crates/connectivity-core/src/legacy/canonical.rs
similarity index 100%
rename from src/legacy/canonical.rs
rename to crates/connectivity-core/src/legacy/canonical.rs
diff --git a/src/legacy/graph.rs b/crates/connectivity-core/src/legacy/graph.rs
similarity index 100%
rename from src/legacy/graph.rs
rename to crates/connectivity-core/src/legacy/graph.rs
diff --git a/src/legacy/metrics.rs b/crates/connectivity-core/src/legacy/metrics.rs
similarity index 100%
rename from src/legacy/metrics.rs
rename to crates/connectivity-core/src/legacy/metrics.rs
diff --git a/src/legacy/mod.rs b/crates/connectivity-core/src/legacy/mod.rs
similarity index 100%
rename from src/legacy/mod.rs
rename to crates/connectivity-core/src/legacy/mod.rs
diff --git a/src/legacy/routing.rs b/crates/connectivity-core/src/legacy/routing.rs
similarity index 100%
rename from src/legacy/routing.rs
rename to crates/connectivity-core/src/legacy/routing.rs
diff --git a/src/legacy/window.rs b/crates/connectivity-core/src/legacy/window.rs
similarity index 100%
rename from src/legacy/window.rs
rename to crates/connectivity-core/src/legacy/window.rs
diff --git a/crates/connectivity-core/src/lib.rs b/crates/connectivity-core/src/lib.rs
new file mode 100644
index 0000000..48a033c
--- /dev/null
+++ b/crates/connectivity-core/src/lib.rs
@@ -0,0 +1,241 @@
+//! Multi-resolution landscape connectivity on in-memory raster arrays.
+//!
+//! This is the numerical engine of the [`eco-connectivity`](https://pypi.org/project/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](https://docs.rs/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)));
+//! ```
+mod affine;
+mod builder;
+mod core;
+mod coverage;
+mod distances;
+mod graph;
+#[cfg(test)]
+mod legacy;
+mod metrics;
+mod overview;
+mod routing;
+mod tables;
+mod template;
+#[cfg(test)]
+mod testkit;
+mod utils;
+#[cfg(test)]
+mod verify;
+mod window;
+
+pub use affine::{level_transforms, Affine};
+pub use coverage::{pixel_coverage, Polygon};
+/// The `ndarray` version of this crate's array types.
+pub use ndarray;
+pub use window::WindowMode;
+
+use ndarray::{Array2, Array3};
+use std::collections::HashMap;
+
+/// Borrowed arrays of one tile for [`connectivity`]. 2-D arrays share the `(rows, cols)`
+/// shape of `condition`.
+#[derive(Clone, Copy, Debug)]
+pub struct Inputs<'a> {
+ /// Habitat condition, expected in [0, 1] with higher values for better condition. It sets
+ /// each destination's value and habitat area, and the path weights unless `traversal` is
+ /// given. NaN cells are left out of habitat area and of the graph; coarser overview cells
+ /// average the finite cells they contain.
+ pub condition: &'a Array2,
+ /// `true` skips a focal cell, whose output stays NaN. Masked cells with finite condition
+ /// still serve as destinations and path cells. Under [`Config::pa_to_pa`], `false` also
+ /// marks protected cells. Unmasked cells with NaN condition reach nothing and score 0, so
+ /// callers normally mask them.
+ pub mask: &'a Array2,
+ /// Optional traversal values in [0, 1], high for easy passage (e.g. `1 - resistance`),
+ /// used only for path weights. `None` uses condition; cells whose traversal overview is
+ /// NaN fall back to condition.
+ pub traversal: Option<&'a Array2>,
+ /// BERI scenarios: the current climate first, then future scenarios. Each array is
+ /// `(rows, cols, bands)` with condition's rows and columns and the same bands in the same
+ /// order; any memory layout works, e.g. a band-major view. A cell whose bands are not all
+ /// finite has similarity 0. `None` or an empty slice computes connectedness instead.
+ pub transgrids: Option<&'a [Array3]>,
+ /// Transform of each level's grid, keyed by level; see [`level_transforms`].
+ pub transforms: &'a HashMap,
+}
+
+/// Settings of one [`connectivity`] call. `Config::default()` gives the Python defaults,
+/// with level 1 included as the core requires.
+#[derive(Clone, Copy, Debug)]
+pub struct Config<'a> {
+ /// Overview levels in base pixels per overview cell: powers of two forming a continuous
+ /// sequence that includes 1, such as `[1, 2, 4, 8, 16, 32]`, in any order.
+ pub levels: &'a [usize],
+ /// Dispersal kernel bandwidths in kilometres; the output is the mean of their scores. An
+ /// empty slice scores every evaluated cell 0.
+ pub lambdas: &'a [f32],
+ /// Whether transforms are geographic (degrees) rather than projected (metres).
+ pub is_geo: bool,
+ /// Relative cost of crossing a fully degraded cell: a cell with traversal value `t` has
+ /// path weight `(1 - max_cost) * t + max_cost`, so `t = 1` costs 1 and `t = 0` costs
+ /// `max_cost`.
+ pub max_cost: f32,
+ /// Odd window size of every level but the coarsest. The level-`L` annulus spans
+ /// `L * window_size / 2` to `L * window_size` base pixels from the focal pixel centre, and
+ /// level 1 is the disc of radius `window_size`; in square mode these are half-widths.
+ pub window_size: i32,
+ /// Odd window size of the coarsest level `L`, whose annulus reaches `L * outer_window`
+ /// base pixels. Values below `window_size` use `window_size`.
+ pub outer_window: i32,
+ /// Global `(row0, col0)` of the tile's top-left pixel in the full raster, in base pixels;
+ /// `(0, 0)` for a whole raster.
+ pub offsets: (usize, usize),
+ /// Circular or square annuli.
+ pub window_mode: WindowMode,
+ /// PARC target gating: only protected destinations, weighted by each overview cell's
+ /// protected fraction (unmasked share of its valid pixels), count in the connectedness
+ /// numerator. Paths and the denominator are unchanged, and a focal cell that reaches no
+ /// protected destination scores 0. For BERI it only makes focal cells without neighbours
+ /// score 0; Python's `beri()` never sets it.
+ pub pa_to_pa: bool,
+}
+
+impl Default for Config<'_> {
+ fn default() -> Self {
+ Self {
+ levels: &[1, 2, 4, 8, 16, 32],
+ lambdas: &[2.0, 20.0, 200.0],
+ is_geo: false,
+ max_cost: 2.0,
+ window_size: 3,
+ outer_window: 9,
+ offsets: (0, 0),
+ window_mode: WindowMode::Circular,
+ pa_to_pa: false,
+ }
+ }
+}
+
+/// Connectedness, or BERI when [`Inputs::transgrids`] is non-empty, of every unmasked cell.
+///
+/// Returns an array with the shape of [`Inputs::condition`]: NaN where the mask is `true`,
+/// otherwise the mean score over [`Config::lambdas`], normally in [0, 1] and 0 for a cell that
+/// reaches no destination.
+///
+/// # Errors
+/// Returns an error when [`Inputs::traversal`] differs in shape from the condition array
+/// (`"Dimension mismatch for key : vs "`). The
+/// Python binding raises errors as `RuntimeError("Connectivity failed: ")`.
+///
+/// # Panics
+/// Invalid settings panic, with messages the Python binding raises as
+/// `pyo3_runtime.PanicException`; panics on worker threads reach the caller through Rayon.
+/// - A level that is not a power of two (`"Failed average resampling.: Level 3 is not a power
+/// of 2"`), or no level-1 transform (`"Missing base-level transform (key=1)."`).
+/// - When at least one cell is unmasked: no level 1 (`"Base condition level not found!"`),
+/// an even `window_size` (`"Window size must be odd-numbered"`) or an even
+/// `outer_window >= window_size` (`"Outer window size must be odd-numbered"`), a level
+/// without a transform (`"Missing level in Affine set."`), or levels that skip a power of
+/// two once a level past the gap has window cells (`"levels must form a continuous
+/// power-of-two sequence"`).
+/// - An integer path cost beyond `u32::MAX`, e.g. from a very large `max_cost` or unscaled
+/// traversal values (`"Path cost overflow; check max_cost and resistance scaling!"`).
+/// - A mask or transgrid smaller than the condition array can panic on indexing.
+pub fn connectivity(inputs: &Inputs, config: &Config) -> anyhow::Result> {
+ core::conn(inputs, config)
+}
diff --git a/src/metrics.rs b/crates/connectivity-core/src/metrics.rs
similarity index 100%
rename from src/metrics.rs
rename to crates/connectivity-core/src/metrics.rs
diff --git a/src/overview.rs b/crates/connectivity-core/src/overview.rs
similarity index 100%
rename from src/overview.rs
rename to crates/connectivity-core/src/overview.rs
diff --git a/src/routing.rs b/crates/connectivity-core/src/routing.rs
similarity index 100%
rename from src/routing.rs
rename to crates/connectivity-core/src/routing.rs
diff --git a/src/tables.rs b/crates/connectivity-core/src/tables.rs
similarity index 99%
rename from src/tables.rs
rename to crates/connectivity-core/src/tables.rs
index 811d044..7c628ea 100644
--- a/src/tables.rs
+++ b/crates/connectivity-core/src/tables.rs
@@ -324,9 +324,10 @@ impl Tables {
#[cfg(test)]
mod tests {
use super::*;
+ use crate::affine::level_transforms;
use crate::distances::distance_km;
use crate::template::CACHE_BUDGET;
- use crate::testkit::{level_transforms, Rng};
+ use crate::testkit::Rng;
#[test]
fn lsb_and_quantum() {
diff --git a/src/template.rs b/crates/connectivity-core/src/template.rs
similarity index 100%
rename from src/template.rs
rename to crates/connectivity-core/src/template.rs
diff --git a/src/testkit.rs b/crates/connectivity-core/src/testkit.rs
similarity index 96%
rename from src/testkit.rs
rename to crates/connectivity-core/src/testkit.rs
index 10e04c8..0673be1 100644
--- a/src/testkit.rs
+++ b/crates/connectivity-core/src/testkit.rs
@@ -1,7 +1,8 @@
//! Test fixtures, a small deterministic PRNG and output comparisons.
-use crate::affine::Affine;
+use crate::affine::{level_transforms, Affine};
use crate::legacy::{self, canonical, Prepared};
use crate::window::WindowMode;
+use crate::{Config, Inputs};
use ndarray::{Array2, Array3};
use std::collections::HashMap;
@@ -39,27 +40,6 @@ impl Rng {
}
}
-/// Level transforms anchored like `rastio.read_raster`:
-/// `base * Affine.scale(level) * Affine.translation(col0 // level, row0 // level)`,
-/// with the `affine` package's operation order.
-pub fn level_transforms(
- base: (f64, f64, f64, f64, f64, f64),
- levels: &[usize],
- offsets: (usize, usize),
-) -> HashMap {
- let (a, b, c, d, e, f) = base;
- levels
- .iter()
- .map(|&level| {
- let l = level as f64;
- let (sa, sb, sd, se) = (a * l + b * 0.0, a * 0.0 + b * l, d * l + e * 0.0, d * 0.0 + e * l);
- let (tc, tr) = ((offsets.1 / level) as f64, (offsets.0 / level) as f64);
- let t = (sa, sb, sa * tc + sb * tr + c, sd, se, sd * tc + se * tr + f);
- (level as i32, Affine::from(t))
- })
- .collect()
-}
-
#[derive(Clone)]
pub struct Fixture {
pub cond: Array2,
@@ -79,6 +59,32 @@ pub struct Fixture {
}
impl Fixture {
+ /// The fixture's arrays as production inputs.
+ pub fn inputs(&self) -> Inputs<'_> {
+ Inputs {
+ condition: &self.cond,
+ mask: &self.mask,
+ traversal: self.trav.as_ref(),
+ transgrids: self.trans.as_deref(),
+ transforms: &self.transforms,
+ }
+ }
+
+ /// The fixture's settings as a production configuration.
+ pub fn config(&self) -> Config<'_> {
+ Config {
+ levels: &self.levels,
+ lambdas: &self.lambdas,
+ is_geo: self.is_geo,
+ max_cost: self.max_cost,
+ window_size: self.window_size,
+ outer_window: self.outer_window,
+ offsets: self.offsets,
+ window_mode: self.mode,
+ pa_to_pa: self.pa_to_pa,
+ }
+ }
+
pub fn prepared(&self) -> Prepared<'_> {
Prepared::new(
&self.cond,
diff --git a/src/utils.rs b/crates/connectivity-core/src/utils.rs
similarity index 100%
rename from src/utils.rs
rename to crates/connectivity-core/src/utils.rs
diff --git a/src/verify.rs b/crates/connectivity-core/src/verify.rs
similarity index 95%
rename from src/verify.rs
rename to crates/connectivity-core/src/verify.rs
index 3d25d4a..54b029d 100644
--- a/src/verify.rs
+++ b/crates/connectivity-core/src/verify.rs
@@ -1,6 +1,6 @@
//! Verification runs on captured workloads (ignored by default; not timed).
//!
-//! `CONN_INPUTS= CONN_OUT= cargo test --release verify -- --ignored --nocapture`
+//! `CONN_INPUTS= CONN_OUT= cargo test --release -p connectivity-core verify -- --ignored --nocapture`
//! reads every `//rust.txt` (or the names in `CONN_WORKLOADS`) and, per workload:
//! compares Optimized with the Canonical reference bit for bit, runs the structural comparison
//! on every unmasked cell, and writes the Optimized, Canonical and Legacy outputs as raw f32.
@@ -35,23 +35,7 @@ fn workloads() -> (PathBuf, PathBuf, Vec) {
}
fn optimized(f: &Fixture) -> (ndarray::Array2, core::Stats) {
- core::conn_stats(
- &Some(f.cond.clone()),
- f.trav.as_ref(),
- f.trans.as_deref(),
- &f.transforms,
- &f.mask,
- &f.levels,
- &f.lambdas,
- f.is_geo,
- f.max_cost,
- f.window_size,
- f.outer_window,
- f.offsets,
- f.mode,
- f.pa_to_pa,
- )
- .unwrap()
+ core::conn_stats(&f.inputs(), &f.config()).unwrap()
}
/// Structural comparison of every unmasked cell with per-edge distance tallies.
diff --git a/src/window.rs b/crates/connectivity-core/src/window.rs
similarity index 95%
rename from src/window.rs
rename to crates/connectivity-core/src/window.rs
index 1208288..3d829b9 100644
--- a/src/window.rs
+++ b/crates/connectivity-core/src/window.rs
@@ -1,12 +1,18 @@
use ndarray::ArrayView1;
+/// 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.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum WindowMode {
+ /// Source-centred circular annuli (the Python default).
Circular,
+ /// Source-centred square annuli.
Square,
}
impl WindowMode {
+ /// Parse `"circular"` or `"square"`, ignoring case and surrounding whitespace. The error
+ /// is the message the Python binding raises as `ValueError`.
pub fn parse(value: &str) -> Result {
match value.trim().to_ascii_lowercase().as_str() {
"square" => Ok(Self::Square),
diff --git a/crates/connectivity-core/tests/api.rs b/crates/connectivity-core/tests/api.rs
new file mode 100644
index 0000000..bf2324a
--- /dev/null
+++ b/crates/connectivity-core/tests/api.rs
@@ -0,0 +1,210 @@
+//! The public API on the multi-level fixture of the unit tests, through crate paths only.
+use connectivity_core::ndarray::{Array2, Array3};
+use connectivity_core::{
+ connectivity, level_transforms, pixel_coverage, Affine, Config, Inputs, WindowMode,
+};
+use std::collections::HashMap;
+
+struct Tile {
+ condition: Array2,
+ mask: Array2,
+ traversal: Array2,
+ transgrids: Vec>,
+ transforms: HashMap,
+}
+
+const LEVELS: [usize; 3] = [1, 2, 4];
+const OFFSETS: (usize, usize) = (3, 5);
+
+/// `testkit::Fixture::multilevel`, built with the public helpers.
+fn tile() -> Tile {
+ let condition = Array2::from_shape_fn((13, 14), |(i, j)| {
+ if (i * 7 + j * 11) % 29 == 0 {
+ f32::NAN
+ } else {
+ ((i * 7 + j * 13) % 10 + 1) as f32 / 10.0
+ }
+ });
+ let mask = Array2::from_shape_fn((13, 14), |(i, j)| {
+ !condition[[i, j]].is_finite() || (i * 3 + j * 5) % 7 == 0
+ });
+ let traversal = Array2::from_shape_fn((13, 14), |(i, j)| {
+ if (i + j * 3) % 13 == 0 {
+ f32::NAN
+ } else {
+ ((i * 5 + j * 3) % 11) as f32 / 10.0
+ }
+ });
+ let transgrids = (0..2)
+ .map(|scenario| {
+ Array3::from_shape_fn((13, 14, 2), |(i, j, k)| {
+ if (i * 3 + j * 5 + scenario + k) % 31 == 0 {
+ f32::NAN
+ } else {
+ ((i * 7 + j * 3 + scenario * 11 + k * 2) % 23) as f32 / 23.0
+ }
+ })
+ })
+ .collect();
+ let base = (1000.0, 0.0, 500_000.0, 0.0, -1000.0, 6_000_000.0);
+ let transforms = level_transforms(base, &LEVELS, OFFSETS);
+ Tile { condition, mask, traversal, transgrids, transforms }
+}
+
+fn inputs<'a>(t: &'a Tile, beri: bool) -> Inputs<'a> {
+ Inputs {
+ condition: &t.condition,
+ mask: &t.mask,
+ traversal: Some(&t.traversal),
+ transgrids: beri.then_some(t.transgrids.as_slice()),
+ transforms: &t.transforms,
+ }
+}
+
+fn config(window_mode: WindowMode, pa_to_pa: bool) -> Config<'static> {
+ Config {
+ levels: &LEVELS,
+ lambdas: &[2.0, 20.0],
+ is_geo: false,
+ max_cost: 20.0,
+ window_size: 3,
+ outer_window: 5,
+ offsets: OFFSETS,
+ window_mode,
+ pa_to_pa,
+ }
+}
+
+fn same_bits(a: &Array2, b: &Array2) -> bool {
+ a.dim() == b.dim() && a.iter().zip(b).all(|(x, y)| x.to_bits() == y.to_bits())
+}
+
+// The frozen sample values of the unit tests' `multilevel_values`
+#[test]
+fn multilevel_samples() {
+ let t = tile();
+ let cells = [(0, 1), (3, 5), (6, 8), (12, 13)];
+ for (mode, beri, pa_to_pa, expected) in [
+ (WindowMode::Circular, false, false, [0.035036724, 0.060928639, 0.057472181, 0.023413222]),
+ (WindowMode::Circular, false, true, [0.031013189, 0.051940750, 0.050136760, 0.017623197]),
+ (WindowMode::Circular, true, false, [0.287769228, 0.330879152, 0.336351663, 0.270322204]),
+ (WindowMode::Square, false, false, [0.036365397, 0.062699854, 0.057887629, 0.023465034]),
+ (WindowMode::Square, false, true, [0.032288730, 0.053552385, 0.050358951, 0.017685590]),
+ (WindowMode::Square, true, false, [0.288273066, 0.334522188, 0.339534730, 0.269480884]),
+ ] {
+ let out = connectivity(&inputs(&t, beri), &config(mode, pa_to_pa)).unwrap();
+ for (cell, expected) in cells.iter().zip(expected) {
+ let actual = out[*cell];
+ assert!((actual - expected).abs() <= 1.0e-6, "{mode:?} {beri} {pa_to_pa} {cell:?}: {actual}");
+ }
+ for (cell, &masked) in t.mask.indexed_iter() {
+ assert_eq!(out[cell].is_nan(), masked);
+ }
+ }
+}
+
+// Any thread count of the caller's pool gives the same bits
+#[test]
+fn thread_pools() {
+ let t = tile();
+ for beri in [false, true] {
+ let run = |threads| {
+ let pool = rayon::ThreadPoolBuilder::new().num_threads(threads).build().unwrap();
+ pool.install(|| connectivity(&inputs(&t, beri), &config(WindowMode::Circular, true)).unwrap())
+ };
+ let one = run(1);
+ assert!(same_bits(&one, &run(3)) && same_bits(&one, &run(8)));
+ }
+}
+
+#[test]
+fn optional_inputs() {
+ let t = tile();
+ let mut plain = inputs(&t, false);
+ plain.traversal = None;
+ let expected = connectivity(&plain, &config(WindowMode::Circular, false)).unwrap();
+ plain.transgrids = Some(&[]);
+ assert!(same_bits(&connectivity(&plain, &config(WindowMode::Circular, false)).unwrap(), &expected));
+ let empty = Config { lambdas: &[], ..config(WindowMode::Circular, false) };
+ let zeros = connectivity(&plain, &empty).unwrap();
+ assert!(t.mask.iter().zip(&zeros).all(|(&m, &v)| if m { v.is_nan() } else { v == 0.0 }));
+}
+
+#[test]
+fn default_config() {
+ let c = Config::default();
+ assert_eq!(c.levels, &[1, 2, 4, 8, 16, 32]);
+ assert_eq!(c.lambdas, &[2.0, 20.0, 200.0]);
+ assert_eq!((c.max_cost, c.window_size, c.outer_window), (2.0, 3, 9));
+ assert_eq!((c.is_geo, c.offsets, c.window_mode, c.pa_to_pa), (false, (0, 0), WindowMode::Circular, false));
+}
+
+#[test]
+fn window_modes() {
+ assert_eq!(WindowMode::parse(" SqUaRe "), Ok(WindowMode::Square));
+ assert_eq!(WindowMode::parse("circular"), Ok(WindowMode::Circular));
+ assert_eq!(
+ WindowMode::parse("hexagon"),
+ Err("Invalid window_mode 'hexagon'. Expected 'square' or 'circular'.".to_string())
+ );
+}
+
+// Level transforms match rasterio's `base * scale(L) * translation(col0 // L, row0 // L)`
+#[test]
+fn tile_transforms() {
+ let base = Affine::from((250.0, 0.0, 500_000.0, 0.0, -250.0, 6_000_000.0));
+ let transforms = level_transforms(base, &[1, 2, 4], (5, 11));
+ for (level, (x0, y0)) in [(1, (502_750.0, 5_998_750.0)), (2, (502_500.0, 5_999_000.0)), (4, (502_000.0, 5_999_000.0))] {
+ let t = transforms[&level];
+ let s = 250.0 * level as f64;
+ assert_eq!((t.x_scale, t.x_skew, t.x_origin, t.y_skew, t.y_scale, t.y_origin), (s, 0.0, x0, 0.0, -s, y0));
+ assert_eq!(t.xy(0, 0), (x0 + s / 2.0, y0 - s / 2.0));
+ }
+ assert_eq!(base.col_row(500_625.0, 5_999_000.0), (2.5, 4.0));
+}
+
+#[test]
+fn traversal_shape() {
+ let t = tile();
+ let short = t.traversal.slice(connectivity_core::ndarray::s![..12, ..]).to_owned();
+ let mut bad = inputs(&t, false);
+ bad.traversal = Some(&short);
+ let error = connectivity(&bad, &config(WindowMode::Circular, false)).unwrap_err();
+ assert!(error.to_string().starts_with("Dimension mismatch for key"), "{error}");
+}
+
+#[test]
+#[should_panic(expected = "Window size must be odd-numbered")]
+fn even_window_panics() {
+ let t = tile();
+ connectivity(&inputs(&t, false), &Config { window_size: 4, ..config(WindowMode::Circular, false) }).unwrap();
+}
+
+#[test]
+#[should_panic(expected = "Missing base-level transform (key=1).")]
+fn needs_base_transform() {
+ let t = tile();
+ let transforms: HashMap = t.transforms.iter().filter(|(&l, _)| l != 1).map(|(&l, &a)| (l, a)).collect();
+ let all_masked = Array2::from_elem((13, 14), true);
+ let inputs = Inputs { mask: &all_masked, transforms: &transforms, ..inputs(&t, false) };
+ connectivity(&inputs, &config(WindowMode::Circular, false)).unwrap();
+}
+
+#[test]
+fn coverage_fractions() {
+ // 10 m pixels from (1000, 2000); pixel (row 1, col 2) spans x 1020..1030, y 1980..1990
+ let transform = Affine::from((10.0, 0.0, 1000.0, 0.0, -10.0, 2000.0));
+ let square = |x0: f64, y0: f64, x1: f64, y1: f64| vec![(x0, y0), (x1, y0), (x1, y1), (x0, y1)];
+ let polygons = vec![
+ (square(1020.0, 1990.0, 1030.0, 1980.0), vec![]),
+ (square(1040.0, 1990.0, 1045.0, 1980.0), vec![]),
+ (square(1000.0, 1960.0, 1040.0, 1940.0), vec![square(1010.0, 1955.0, 1030.0, 1945.0)]),
+ ];
+ let out = pixel_coverage(&polygons, &transform, 6, 6);
+ assert_eq!((out[[1, 2]], out[[1, 4]]), (1.0, 0.5));
+ for row in 4..6 {
+ assert_eq!(out.row(row).to_vec(), vec![1.0, 0.5, 0.5, 1.0, 0.0, 0.0]);
+ }
+ assert_eq!(out.sum(), 7.5);
+ assert_eq!(pixel_coverage(&[], &transform, 2, 3), Array2::::zeros((2, 3)));
+}
diff --git a/src/affine.rs b/src/affine.rs
deleted file mode 100644
index 94716d2..0000000
--- a/src/affine.rs
+++ /dev/null
@@ -1,48 +0,0 @@
-use std::convert::From;
-
-// Define the Affine struct for transforming ij to xy
-#[derive(Debug, Clone, Copy)]
-pub struct Affine {
- pub x_scale: f64, // pixel width
- pub x_skew: f64, // rotation/skew term (usually 0)
- pub x_origin: f64, // top-left corner X
- pub y_skew: f64, // rotation/skew term (usually 0)
- pub y_scale: f64, // pixel height (negative if north-up)
- pub y_origin: f64, // top-left corner Y
-}
-
-impl Affine {
- /// Convert (row, col) to (x, y) at pixel centre by default.
- pub fn xy(&self, row: i32, col: i32) -> (f64, f64) {
- let col_f = (col as f64) + 0.5;
- let row_f = (row as f64) + 0.5;
-
- let x = self.x_scale * col_f + self.x_skew * row_f + self.x_origin;
- let y = self.y_skew * col_f + self.y_scale * row_f + self.y_origin;
- (x, y)
- }
-
- /// Convert world (x, y) to continuous (col, row) pixel coordinates.
- pub fn col_row(&self, x: f64, y: f64) -> (f64, f64) {
- let det = self.x_scale * self.y_scale - self.x_skew * self.y_skew;
- let dx = x - self.x_origin;
- let dy = y - self.y_origin;
- let col = (self.y_scale * dx - self.x_skew * dy) / det;
- let row = (-self.y_skew * dx + self.x_scale * dy) / det;
- (col, row)
- }
-}
-
-impl From<(f64, f64, f64, f64, f64, f64)> for Affine {
- fn from(t: (f64, f64, f64, f64, f64, f64)) -> Self {
- Self {
- x_scale: t.0,
- x_skew: t.1,
- x_origin: t.2,
- y_skew: t.3,
- y_scale: t.4,
- y_origin: t.5,
- }
- }
-}
-
diff --git a/src/extract.rs b/src/extract.rs
index b3ade69..9772a61 100644
--- a/src/extract.rs
+++ b/src/extract.rs
@@ -4,7 +4,7 @@ use pyo3::Bound;
use pyo3::types::{PyAny, PyAnyMethods, PyList, PyListMethods};
use std::collections::HashMap;
use ndarray::{Array3, Array2};
-use crate::affine::Affine;
+use connectivity_core::Affine;
/// Convert a Python dict of rasterio transforms into a Rust HashMap.
diff --git a/src/lib.rs b/src/lib.rs
index eb26658..06229bf 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,3 +1,4 @@
+use connectivity_core::{Affine, Config, Inputs, WindowMode};
use ndarray::{Array2, Array3};
use numpy::{PyArray2, ToPyArray};
use pyo3::exceptions::{PyRuntimeError, PyValueError};
@@ -6,27 +7,7 @@ use pyo3::types::{PyAny, PyList};
use pyo3::Bound;
use std::collections::HashMap;
// local modules
-mod affine;
-mod builder;
-mod core;
-mod coverage;
-mod distances;
mod extract;
-mod graph;
-#[cfg(test)]
-mod legacy;
-mod metrics;
-mod overview;
-mod routing;
-mod tables;
-mod template;
-#[cfg(test)]
-mod testkit;
-mod utils;
-#[cfg(test)]
-mod verify;
-mod window;
-use affine::Affine;
/// Compute habitat/PARC connectivity and BERI.
///
@@ -49,7 +30,7 @@ use affine::Affine;
/// (e.g. 3 gives an effective 6×6 current-level window).
/// * `outer_window` - Odd coarsest-level window width used to set the long-range
/// search reach. Must be ≥ `window_size`.
-/// * `offsets` -
+/// * `offsets` - (row, col) of the tile's top-left pixel in the full raster, in base pixels.
/// * `window_mode` - "circular" for source-centered circular annuli or
/// "square" for source-centered square annuli.
/// * `traversal` - Optional 2D array (same shape as `condition`) of pre-inverted, [0, 1]
@@ -117,7 +98,7 @@ fn connectivity(
None => None,
};
- let window_mode = window::WindowMode::parse(window_mode).map_err(PyValueError::new_err)?;
+ let window_mode = WindowMode::parse(window_mode).map_err(PyValueError::new_err)?;
// Set the number of cores for parallel processing with Rayon
// Use a local pool to safely control exact number of core
@@ -126,26 +107,34 @@ fn connectivity(
.build()
.map_err(|e| PyRuntimeError::new_err(e.to_string()))?;
+ // The core requires a condition array; a missing one fails as the core call used to
+ let Some(cond_array) = cond_array else {
+ return Err(PyRuntimeError::new_err(
+ "Connectivity failed: No base level (key=1) in condition array",
+ ));
+ };
+ let inputs = Inputs {
+ condition: &cond_array,
+ mask: &mask_array,
+ traversal: trav_array.as_ref(),
+ transgrids: trans_arrays.as_deref(),
+ transforms: &transform_map,
+ };
+ let config = Config {
+ levels: &levels,
+ lambdas: &lambdas,
+ is_geo,
+ max_cost,
+ window_size,
+ outer_window,
+ offsets,
+ window_mode,
+ pa_to_pa,
+ };
+
// Core connectivity function
let outarray: Array2 = pool
- .install(|| {
- core::conn(
- &cond_array,
- trav_array.as_ref(),
- trans_arrays.as_deref(),
- &transform_map,
- &mask_array,
- &levels,
- &lambdas,
- is_geo,
- max_cost,
- window_size,
- outer_window,
- offsets,
- window_mode,
- pa_to_pa,
- )
- })
+ .install(|| connectivity_core::connectivity(&inputs, &config))
.map_err(|er| PyRuntimeError::new_err(format!("Connectivity failed: {er}")))?;
// Convert the array back to Python with gil
@@ -186,7 +175,7 @@ fn pixel_coverage(
.build()
.map_err(|e| PyRuntimeError::new_err(e.to_string()))?;
- let outarray = pool.install(|| coverage::pixel_coverage_array(&rings, &tr, nrows, ncols));
+ let outarray = pool.install(|| connectivity_core::pixel_coverage(&rings, &tr, nrows, ncols));
Python::attach(|py| {
let pyarray = outarray.to_pyarray(py);
diff --git a/tests/test_bindings.py b/tests/test_bindings.py
new file mode 100644
index 0000000..79bd2dc
--- /dev/null
+++ b/tests/test_bindings.py
@@ -0,0 +1,202 @@
+"""Regression tests for the rust_conn bindings on small synthetic arrays.
+
+The expected values are the frozen samples of the Rust ``multilevel_values`` test
+(crates/connectivity-core/src/core.rs), reached here through the Python binding, so they
+also check how arguments are passed to the core.
+"""
+
+import numpy as np
+import pytest
+
+import rust_conn
+
+CELLS = [(0, 1), (3, 5), (6, 8), (12, 13)]
+EXPECTED = {
+ ("circular", False, False): [0.035036724, 0.060928639, 0.057472181, 0.023413222],
+ ("circular", False, True): [0.031013189, 0.051940750, 0.050136760, 0.017623197],
+ ("circular", True, False): [0.287769228, 0.330879152, 0.336351663, 0.270322204],
+ ("circular", True, True): [0.287769228, 0.330879152, 0.336351663, 0.270322204],
+ ("square", False, False): [0.036365397, 0.062699854, 0.057887629, 0.023465034],
+ ("square", False, True): [0.032288730, 0.053552385, 0.050358951, 0.017685590],
+ ("square", True, False): [0.288273066, 0.334522188, 0.339534730, 0.269480884],
+ ("square", True, True): [0.288273066, 0.334522188, 0.339534730, 0.269480884],
+}
+F32 = np.float32
+
+
+def multilevel(mode="circular", beri=False, pa_to_pa=False):
+ """Arguments of the Rust ``Fixture::multilevel`` fixture."""
+ i, j = np.indices((13, 14))
+ cond = ((i * 7 + j * 13) % 10 + 1).astype(F32) / F32(10)
+ cond[(i * 7 + j * 11) % 29 == 0] = np.nan
+ mask = ~np.isfinite(cond) | ((i * 3 + j * 5) % 7 == 0)
+ traversal = ((i * 5 + j * 3) % 11).astype(F32) / F32(10)
+ traversal[(i + j * 3) % 13 == 0] = np.nan
+ grids = None
+ if beri:
+ i3, j3, k3 = np.indices((13, 14, 2))
+ grids = []
+ for s in range(2):
+ grid = ((i3 * 7 + j3 * 3 + s * 11 + k3 * 2) % 23).astype(F32) / F32(23)
+ grid[(i3 * 3 + j3 * 5 + s + k3) % 31 == 0] = np.nan
+ grids.append(grid)
+ levels, offsets = [1, 2, 4], (3, 5)
+ transforms = {}
+ for level in levels:
+ scale = 1000.0 * level
+ row0, col0 = offsets[0] // level * level, offsets[1] // level * level
+ transforms[level] = (scale, 0.0, 500_000.0 + col0 * 1000.0, 0.0, -scale, 6_000_000.0 - row0 * 1000.0)
+ return dict(condition=cond, mask=mask, transgrid_list=grids, transforms=transforms, levels=levels,
+ lambdas=[2.0, 20.0], is_geo=False, max_cost=20.0, window_size=3, outer_window=5,
+ offsets=offsets, traversal=traversal, window_mode=mode, pa_to_pa=pa_to_pa)
+
+
+def run(**changes):
+ args = multilevel()
+ args.update(changes)
+ return rust_conn.connectivity(**args)
+
+
+def bits(a):
+ return np.ascontiguousarray(a).view(np.uint32)
+
+
+@pytest.mark.parametrize("key", sorted(EXPECTED))
+def test_multilevel_values(key):
+ args = multilevel(*key)
+ out = rust_conn.connectivity(**args)
+ assert out.dtype == np.float32 and out.shape == (13, 14)
+ np.testing.assert_allclose([out[c] for c in CELLS], EXPECTED[key], rtol=0, atol=1e-6)
+ np.testing.assert_array_equal(np.isnan(out), args["mask"])
+
+
+def test_band_major():
+ args = multilevel(beri=True)
+ expected = rust_conn.connectivity(**args)
+ # np.moveaxis views, as read_raster returns multi-band rasters
+ args["transgrid_list"] = [np.moveaxis(np.ascontiguousarray(np.moveaxis(g, -1, 0)), 0, -1)
+ for g in args["transgrid_list"]]
+ assert not args["transgrid_list"][0].flags.c_contiguous
+ np.testing.assert_array_equal(bits(rust_conn.connectivity(**args)), bits(expected))
+
+
+@pytest.mark.parametrize("beri", [False, True])
+def test_threads_and_repeats(beri):
+ args = multilevel(beri=beri, pa_to_pa=True)
+ expected = rust_conn.connectivity(**args, n_threads=1)
+ for n_threads in (None, 0, 2, 4, 1):
+ np.testing.assert_array_equal(bits(rust_conn.connectivity(**{**args, "n_threads": n_threads})),
+ bits(expected))
+
+
+def test_optional_inputs():
+ plain = run(traversal=None)
+ np.testing.assert_array_equal(bits(run(traversal=None, transgrid_list=[])), bits(plain))
+ assert not np.array_equal(bits(run()), bits(plain))
+ beri = multilevel(beri=True)
+ grids = beri["transgrid_list"]
+ expected = rust_conn.connectivity(**beri)
+ beri["transgrid_list"] = [None, grids[0], None, grids[1]]
+ np.testing.assert_array_equal(bits(rust_conn.connectivity(**beri)), bits(expected))
+
+
+def test_mode_parsing():
+ np.testing.assert_array_equal(bits(run(window_mode=" SqUaRe ")), bits(run(window_mode="square")))
+ with pytest.raises(ValueError, match=r"^Invalid window_mode 'hexagon'. Expected 'square' or 'circular'.$"):
+ run(window_mode="hexagon")
+
+
+def test_masks_and_nan():
+ args = multilevel()
+ assert np.isnan(run(mask=np.ones((13, 14), dtype=bool))).all()
+ # Unmasked cells without condition reach nothing and score 0
+ out = run(mask=np.zeros((13, 14), dtype=bool))
+ assert (out[np.isnan(args["condition"])] == 0).all()
+ assert np.isfinite(out).all()
+
+
+def test_no_lambdas():
+ out = run(lambdas=[])
+ mask = multilevel()["mask"]
+ assert (out[~mask] == 0).all() and np.isnan(out[mask]).all()
+
+
+@pytest.mark.parametrize("changes, error, message", [
+ (dict(condition=None), RuntimeError, "Connectivity failed: No base level (key=1) in condition array"),
+ (dict(mask=None), RuntimeError, "Reading mask failed: Mask must be provided!"),
+ (dict(condition=np.zeros((13, 14))), RuntimeError, "Reading condition failed: Failed to extract 2D array"),
+ (dict(mask=np.zeros((13, 14), dtype=np.int8)), RuntimeError, "Reading mask failed: Failed to extract mask array"),
+ (dict(traversal=np.zeros((13, 14))), RuntimeError, "Reading traversal failed: Failed to extract 2D array"),
+ (dict(traversal=np.zeros((12, 14), dtype=F32)), RuntimeError,
+ "Connectivity failed: Dimension mismatch for key"),
+ (dict(transforms=[(1, (1.0, 0.0, 0.0, 0.0, -1.0, 0.0))]), RuntimeError,
+ "Faild getting transform information: AttributeError"),
+ (dict(transgrid_list=()), TypeError, "'tuple' object is not an instance of 'list'"),
+])
+def test_input_errors(changes, error, message):
+ with pytest.raises(error) as info:
+ run(**changes)
+ assert type(info.value) is error
+ assert str(info.value).startswith(message)
+
+
+@pytest.mark.parametrize("changes, message", [
+ (dict(window_size=4), "Window size must be odd-numbered"),
+ (dict(window_size=3, outer_window=6), "Outer window size must be odd-numbered"),
+ (dict(levels=[1, 3]), "Failed average resampling.: Level 3 is not a power of 2"),
+ (dict(levels=[2, 4]), "Base condition level not found!"),
+ (dict(levels=[1, 4]), "levels must form a continuous power-of-two sequence"),
+ (dict(transforms={1: (1000.0, 0.0, 0.0, 0.0, -1000.0, 0.0)}), "Missing level in Affine set."),
+ (dict(transforms={2: (2000.0, 0.0, 0.0, 0.0, -2000.0, 0.0)}), "Missing base-level transform (key=1)."),
+ (dict(max_cost=3.0e7), "Path cost overflow; check max_cost and resistance scaling!"),
+ (dict(transgrid_list=[np.zeros((13, 14, 2))]),
+ "called `Result::unwrap()` on an `Err` value: Failed to extract 3D array"),
+])
+def test_setting_panics(changes, message):
+ # pyo3_runtime.PanicException derives from BaseException, not Exception
+ with pytest.raises(BaseException) as info:
+ run(**changes)
+ assert type(info.value).__name__ == "PanicException"
+ assert not isinstance(info.value, Exception)
+ assert str(info.value).startswith(message)
+
+
+def test_small_outer_window():
+ np.testing.assert_array_equal(bits(run(window_size=5, outer_window=3)),
+ bits(run(window_size=5, outer_window=5)))
+
+
+TRANSFORM = (10.0, 0.0, 1000.0, 0.0, -10.0, 2000.0)
+
+
+def ring(*xy):
+ return np.asarray(xy, dtype=np.float64)
+
+
+def test_coverage_fractions():
+ # Pixel (row 1, col 2) spans x 1020..1030 and y 1980..1990
+ full = ring((1020, 1990), (1030, 1990), (1030, 1980), (1020, 1980))
+ half = ring((1040, 1990), (1045, 1990), (1045, 1980), (1040, 1980))
+ outer = ring((1000, 1960), (1040, 1960), (1040, 1940), (1000, 1940), (1000, 1960))
+ hole = ring((1010, 1955), (1030, 1955), (1030, 1945), (1010, 1945))
+ out = rust_conn.pixel_coverage([(full, []), (half, []), (outer, [hole])], TRANSFORM, (6, 6))
+ assert out.dtype == np.float32 and out.shape == (6, 6)
+ assert out[1, 2] == 1 and out[1, 4] == 0.5
+ np.testing.assert_array_equal(out[4:], [[1.0, 0.5, 0.5, 1.0, 0.0, 0.0]] * 2)
+ assert out.sum() == pytest.approx(1 + 0.5 + 8 - 2)
+ for n_threads in (1, 3):
+ np.testing.assert_array_equal(
+ rust_conn.pixel_coverage([(full, []), (half, []), (outer, [hole])], TRANSFORM, (6, 6), n_threads),
+ out)
+ np.testing.assert_array_equal(rust_conn.pixel_coverage([], TRANSFORM, (2, 3)), np.zeros((2, 3), F32))
+
+
+@pytest.mark.parametrize("polygons, message", [
+ ((), "Reading polygons failed: Expected a list of polygons"),
+ ([(np.zeros((4, 3)), [])], "Reading polygons failed: Ring array must have shape (N, 2), got [4, 3]"),
+ ([(np.zeros((4, 2)), ())], "Reading polygons failed: Holes must be a list of arrays"),
+])
+def test_coverage_errors(polygons, message):
+ with pytest.raises(RuntimeError) as info:
+ rust_conn.pixel_coverage(polygons, TRANSFORM, (6, 6))
+ assert str(info.value).startswith(message)
diff --git a/tests/test_workflows.py b/tests/test_workflows.py
new file mode 100644
index 0000000..24d9724
--- /dev/null
+++ b/tests/test_workflows.py
@@ -0,0 +1,173 @@
+"""Regression tests for the Python workflows on the bundled Tasmania example data.
+
+Exact checks (NaN placement, bit-identical repeats, seamless tiles, Connected Habitat
+post-processing) are platform independent; the frozen means come from the 2.4.0 baseline
+and allow for floating-point differences between platforms.
+"""
+
+import numpy as np
+import pytest
+import rasterio
+
+import connectivity as cn
+from connectivity import example_data_path
+from connectivity.rastio import read_raster
+from connectivity.utils import _crop_array
+
+COND = example_data_path("site_condition")
+PA = example_data_path("pa_proportion")
+CURRENT = example_data_path("transgrids/1990")
+FUTURES = [example_data_path(f"transgrids/{n}") for n in ("ACC50_85", "GFD50_85", "IPS50_45")]
+THREADS = 4
+SHAPE = (516, 588)
+RTOL = 1e-5
+
+# Frozen from the 2.4.0 baseline build (default parameters unless noted)
+FROZEN = {
+ "connectedness": 0.3594061735617171,
+ "parc": 0.6910263746707314,
+ "parc_pa_to_pa": 0.5578077798671234,
+ "beri": 0.40609480758632494,
+ "beri_current": 0.46590158319867997,
+ "resistance": 0.3100888245434447,
+ "square_short": 0.4162536757680946, # square windows, levels [2, 4, 8], outer_window 5
+}
+
+
+def bits(a):
+ return np.ascontiguousarray(np.ma.getdata(a)).view(np.uint32)
+
+
+def nanmean(a):
+ return float(np.nanmean(np.asarray(np.ma.getdata(a), dtype=np.float64)))
+
+
+@pytest.fixture(scope="module")
+def condition():
+ return read_raster(COND, levels=[1])[0]
+
+
+@pytest.fixture(scope="module")
+def connectedness():
+ return cn.connectedness(condition_file=COND, option=1, n_threads=THREADS)
+
+
+def test_habitat_options(condition, connectedness):
+ assert connectedness.dtype == np.float32 and connectedness.shape == SHAPE
+ np.testing.assert_array_equal(np.isnan(connectedness), np.isnan(condition))
+ finite = connectedness[np.isfinite(connectedness)]
+ assert finite.min() >= 0 and finite.max() <= 1
+ assert nanmean(connectedness) == pytest.approx(FROZEN["connectedness"], rel=RTOL)
+ # Connected Habitat stays in Python: option 2 = condition x connectedness, 3 = its sqrt
+ option2 = cn.connectedness(condition_file=COND, option=2, n_threads=THREADS)
+ option3 = cn.connectedness(condition_file=COND, n_threads=THREADS)
+ np.testing.assert_array_equal(bits(option2), bits(condition * connectedness))
+ np.testing.assert_array_equal(bits(option3), bits(np.sqrt(np.multiply(condition, connectedness, dtype=np.float32))))
+ with pytest.raises(ValueError, match="option must be one of 1, 2, or 3"):
+ cn.connectedness(condition_file=COND, option=4, levels=[2, 4, 8], outer_window=5, n_threads=THREADS)
+
+
+def test_thread_counts():
+ kwargs = dict(condition_file=COND, levels=[2, 4, 8], outer_window=5, option=1)
+ expected = cn.connectedness(**kwargs, n_threads=1)
+ for n_threads in (2, THREADS, None):
+ np.testing.assert_array_equal(bits(cn.connectedness(**kwargs, n_threads=n_threads)), bits(expected))
+
+
+def test_window_modes(connectedness):
+ square = cn.connectedness(condition_file=COND, option=1, window_mode="square", levels=[2, 4, 8],
+ outer_window=5, n_threads=THREADS)
+ np.testing.assert_array_equal(np.isnan(square), np.isnan(connectedness))
+ assert nanmean(square) == pytest.approx(FROZEN["square_short"], rel=RTOL)
+ with pytest.raises(ValueError, match="window_mode must be one of"):
+ cn.connectedness(condition_file=COND, window_mode="block")
+
+
+def test_parc_and_pa_to_pa(condition):
+ pa = read_raster(PA, levels=[1])[0]
+ parc = cn.connectedness(condition_file=COND, pa_file=PA, n_threads=THREADS)
+ gated = cn.connectedness(condition_file=COND, pa_file=PA, pa_to_pa=True, n_threads=THREADS)
+ # PARC evaluates protected cells only; the output is connectedness (no Connected Habitat)
+ protected = (pa > 0) & ~np.isnan(np.maximum(condition, pa))
+ np.testing.assert_array_equal(np.isfinite(parc), protected)
+ np.testing.assert_array_equal(np.isfinite(gated), protected)
+ assert nanmean(parc) == pytest.approx(FROZEN["parc"], rel=RTOL)
+ assert nanmean(gated) == pytest.approx(FROZEN["parc_pa_to_pa"], rel=RTOL)
+ # Gating only removes unprotected destinations from the numerator
+ assert (gated[protected] <= parc[protected]).all()
+ assert (gated[protected] < parc[protected]).any()
+
+
+def test_beri(connectedness):
+ out = cn.beri(condition_file=COND, current_file=CURRENT, future_files=FUTURES, n_threads=THREADS)
+ current = cn.beri(condition_file=COND, current_file=CURRENT, n_threads=THREADS)
+ for result in (out, current):
+ assert result.dtype == np.float32 and result.shape == SHAPE
+ np.testing.assert_array_equal(np.isnan(result), np.isnan(connectedness))
+ assert (result[np.isfinite(result)] >= 0).all()
+ assert nanmean(out) == pytest.approx(FROZEN["beri"], rel=RTOL)
+ assert nanmean(current) == pytest.approx(FROZEN["beri_current"], rel=RTOL)
+ assert nanmean(out) < nanmean(current)
+
+
+def write_resistance(path):
+ """Resistance rising from west to east, missing in a few rows of valid condition."""
+ with rasterio.open(COND) as src:
+ profile = src.profile | {"nodata": np.nan, "dtype": "float32"}
+ valid = np.isfinite(src.read(1, masked=True).filled(np.nan))
+ cols = np.broadcast_to(np.arange(SHAPE[1], dtype=np.float32), SHAPE)
+ resistance = np.where(valid, np.float32(0.2) + np.float32(0.6) * cols / np.float32(SHAPE[1]), np.nan)
+ resistance[300:303] = np.nan
+ with rasterio.open(path, "w", **profile) as dst:
+ dst.write(resistance.astype(np.float32), 1)
+ return int((valid[300:303]).sum())
+
+
+def test_resistance_surface(tmp_path, connectedness):
+ res_file = str(tmp_path / "resistance.tif")
+ gaps = write_resistance(res_file)
+ with pytest.warns(UserWarning, match=f"^{gaps} cell\\(s\\) have valid condition but missing resistance"):
+ out = cn.connectedness(condition_file=COND, resistance_file=res_file, option=1, n_threads=THREADS)
+ np.testing.assert_array_equal(np.isnan(out), np.isnan(connectedness))
+ assert not np.array_equal(bits(out), bits(connectedness))
+ assert nanmean(out) == pytest.approx(FROZEN["resistance"], rel=RTOL)
+
+
+def test_seamless_tiles():
+ whole = cn.connectedness(condition_file=COND, n_threads=THREADS)
+ with rasterio.open(COND) as src:
+ transform = src.transform
+ tiles = cn.make_tiles(COND, nrows=2, ncols=2, align_to=32)
+ assert list(tiles["row_start"]) == [0, 0, 256, 256] and list(tiles["col_start"]) == [0, 288, 0, 288]
+ for k in range(len(tiles)):
+ tile = tiles.iloc[[k]]
+ out = cn.connectedness(condition_file=COND, polygon_mask=tile, n_threads=THREADS)
+ expected, _ = _crop_array(whole, transform, tile)
+ assert np.ma.isMaskedArray(out) and out.shape == expected.shape
+ np.testing.assert_array_equal(np.ma.getmaskarray(out), np.ma.getmaskarray(expected))
+ # Open-border tiles see the same overview cells and windows as the whole raster
+ inside = ~np.ma.getmaskarray(out)
+ np.testing.assert_array_equal(bits(out)[inside], bits(expected)[inside])
+ closed = cn.connectedness(condition_file=COND, polygon_mask=tile, closed_border=True, n_threads=THREADS)
+ assert np.isfinite(closed).any()
+
+
+def test_pixel_coverage(tmp_path):
+ import geopandas as gpd
+ from shapely.geometry import Point, box
+ from shapely.ops import unary_union
+
+ shapes = [Point(146.5, -42.0).buffer(0.8).difference(Point(146.4, -42.1).buffer(0.3)),
+ box(144.0, -41.0, 145.0, -40.5)]
+ polygons = gpd.GeoDataFrame(geometry=shapes, crs="EPSG:4326")
+ out_file = str(tmp_path / "coverage.tif")
+ coverage = cn.pixel_coverage(polygons, COND, n_threads=THREADS, filename=out_file)
+ assert coverage.dtype == np.float32 and coverage.shape == SHAPE
+ assert coverage.min() == 0 and coverage.max() == 1
+ with rasterio.open(COND) as src:
+ pixel_area = abs(src.transform.a * src.transform.e)
+ assert coverage.sum(dtype=np.float64) * pixel_area == pytest.approx(unary_union(shapes).area, rel=1e-5)
+ for n_threads in (1, None):
+ np.testing.assert_array_equal(cn.pixel_coverage(polygons, COND, n_threads=n_threads), coverage)
+ with rasterio.open(out_file) as written:
+ np.testing.assert_array_equal(written.read(1), coverage)