Skip to content

Fix zonal statistics validation and review feedback #838

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
giswqs merged 2 commits into from
Oct 24, 2026
Merged

Fix zonal statistics validation and review feedback #838

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5e775d5
fix: address zonal statistics review feedback
giswqs Oct 24, 2026
bb52f3e
Merge remote-tracking branch 'origin/main' into pr-838
giswqs Oct 24, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
161 changes: 86 additions & 75 deletions crates/analysis/src/rasterops.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -523,29 +523,6 @@ pub fn zonal_means(r: &Raster, bands: &[usize], rings: &[Vec<Coord>]) -> Option<
(n > 0).then(|| (sums.into_iter().map(|s| s / n as f64).collect(), n))
}

/// Transforms one polygon part, keeping its exterior and holes together.
/// A failed coordinate or degenerate ring (< 3 points) rejects the entire polygon part.
pub fn transform_polygon(rings: &[Vec<Coord>], mut transform: impl FnMut(Coord) -> Option<Coord>) -> Option<Vec<Vec<Coord>>> {
if rings.first()?.len() < 3 {
return None;
}
let mut out = Vec::new();
for ring in rings {
if ring.len() < 3 {
return None;
}
let transformed: Option<Vec<Coord>> = ring
.iter()
.map(|&c| {
let p = transform(c)?;
(p.x.is_finite() && p.y.is_finite()).then_some(p)
})
.collect();
out.push(transformed?);
}
Some(out)
}

/// An error encountered while validating or transforming a polygon geometry.
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum PolygonTransformError {
Expand Down Expand Up @@ -643,6 +620,24 @@ pub fn transform_feature_polygons_to_lonlat(
Ok(out)
}

/// Validates a 1-based band number and returns its 0-based index.
pub(crate) fn band_from_f64(value: f64, total_bands: usize) -> Result<usize, String> {
if !value.is_finite() {
return Err("band number must be a finite number".into());
}
if value <= 0.0 {
return Err("band numbers must be positive (1-based)".into());
}
if value.fract() != 0.0 {
return Err("band numbers must be integers, not decimals".into());
}
// Check before casting, including on wasm32, where huge values would saturate.
if value > total_bands as f64 || value >= usize::MAX as f64 {
return Err(format!("band index {value} is invalid: raster has bands 1 to {total_bands}"));
}
Ok(value as usize - 1)
}

/// Parses a band specification string into 0-based band indices.
///
/// Supported inputs:
Expand Down Expand Up @@ -691,20 +686,7 @@ pub fn parse_bands_str(s: &str, total_bands: usize) -> Result<Vec<usize>, String
continue;
}
let f: f64 = subpart.parse().map_err(|_| format!("invalid band number '{subpart}'"))?;
if !f.is_finite() {
return Err("band number must be a finite number".into());
}
if f <= 0.0 {
return Err("band numbers must be positive (1-based)".into());
}
if f.fract() != 0.0 {
return Err("band numbers must be integers, not decimals".into());
}
let b = f as usize;
if b > total_bands {
return Err(format!("band index {b} is invalid: raster has bands 1 to {total_bands}"));
}
let b_idx = b - 1;
let b_idx = band_from_f64(f, total_bands)?;
if !out.contains(&b_idx) {
out.push(b_idx);
}
Expand All @@ -716,53 +698,51 @@ pub fn parse_bands_str(s: &str, total_bands: usize) -> Result<Vec<usize>, String
Ok(out)
}

/// Projects a complete ring to the raster grid without dropping failed vertices.
fn project_ring(r: &Raster, ring: &[Coord]) -> Option<Vec<(f64, f64)>> {
if ring.len() < 3 {
return None;
}
let tr = r.transformer();
ring.iter()
.map(|&c| {
if !c.is_finite() {
return None;
}
let pt = r.transform.to_pixel(tr.from_lonlat(c)?);
(pt.0.is_finite() && pt.1.is_finite()).then_some(pt)
Comment thread
coderabbitai[bot] marked this conversation as resolved.
})
.collect()
}

/// The cells of `r` whose centres fall inside any of `polygons` (each an exterior ring and its
/// holes, in longitude/latitude): the column and row window around the polygons, clamped to the
/// raster, and a row-major mask over it. `None` when no cell centre falls inside.
///
/// The rings are projected onto the raster's own grid once and filled by scanline, so a large
/// raster costs a pass over the polygon edges per row rather than a point test per cell.
/// A malformed ring or a failed projection rejects the entire polygon part, including its holes.
pub fn polygon_mask(r: &Raster, polygons: &[Vec<Vec<Coord>>]) -> Option<(Range<u32>, Range<u32>, Vec<bool>)> {
let tr = r.transformer();
polygon_mask_with_failures(r, polygons, &mut Vec::new())
}

/// Builds a polygon mask and appends the zero-based indices of rejected parts to `failed_parts`.
/// Valid parts are retained even if another part cannot be projected onto the raster grid.
pub fn polygon_mask_with_failures(
r: &Raster,
polygons: &[Vec<Vec<Coord>>],
failed_parts: &mut Vec<usize>,
) -> Option<(Range<u32>, Range<u32>, Vec<bool>)> {
let projected: Vec<Vec<Vec<(f64, f64)>>> = polygons
.iter()
.filter_map(|rings| {
let exterior = rings.first()?;
let mut ext_pts = Vec::with_capacity(exterior.len());
for c in exterior {
if !c.is_finite() {
return None;
}
let w = tr.from_lonlat(*c)?;
let pt = r.transform.to_pixel(w);
if !pt.0.is_finite() || !pt.1.is_finite() {
return None;
}
ext_pts.push(pt);
.enumerate()
.filter_map(|(index, rings)| {
// Reject a complete part if any ring fails, so a lost hole cannot enlarge it.
let part = if rings.is_empty() { None } else { rings.iter().map(|ring| project_ring(r, ring)).collect::<Option<Vec<_>>>() };
if part.is_none() {
failed_parts.push(index);
}
if ext_pts.len() < 3 {
return None;
}
let mut valid_rings = vec![ext_pts];
for hole in rings.iter().skip(1) {
let mut hole_pts = Vec::with_capacity(hole.len());
for c in hole {
if !c.is_finite() {
return None;
}
let w = tr.from_lonlat(*c)?;
let pt = r.transform.to_pixel(w);
if !pt.0.is_finite() || !pt.1.is_finite() {
return None;
}
hole_pts.push(pt);
}
if hole_pts.len() < 3 {
return None;
}
valid_rings.push(hole_pts);
}
Some(valid_rings)
part
})
.collect();
let (mut c0, mut r0, mut c1, mut r1) = (f64::INFINITY, f64::INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
Expand Down Expand Up @@ -918,6 +898,16 @@ impl Accumulator {
/// excluded. If no valid cells fall inside the polygon or if the polygon is outside the raster extent,
/// returns `ZoneStats` with `count = 0` and all numeric statistics set to `None`.
pub fn zonal_statistics(r: &Raster, bands: &[usize], polygons: &[Vec<Vec<Coord>>]) -> Result<Vec<ZoneStats>, String> {
zonal_statistics_with_failures(r, bands, polygons, &mut Vec::new())
}

/// Computes zonal statistics, appending raster-side projection failures as zero-based part indices.
pub(crate) fn zonal_statistics_with_failures(
r: &Raster,
bands: &[usize],
polygons: &[Vec<Vec<Coord>>],
failed_parts: &mut Vec<usize>,
) -> Result<Vec<ZoneStats>, String> {
r.require_coordinates()?;
for &b in bands {
if b >= r.band_count() {
Expand All @@ -928,7 +918,7 @@ pub fn zonal_statistics(r: &Raster, bands: &[usize], polygons: &[Vec<Vec<Coord>>
return Ok(Vec::new());
}

let Some((cols, rows, mask)) = polygon_mask(r, polygons) else {
let Some((cols, rows, mask)) = polygon_mask_with_failures(r, polygons, failed_parts) else {
return Ok(vec![ZoneStats::default(); bands.len()]);
};

Expand Down Expand Up @@ -1091,6 +1081,27 @@ mod tests {
assert!(polygon_mask(&r, &[vec![sq(0.0051, 0.0051, 0.0052, 0.0052)]]).is_none());
}

/// A failed hole drops only its containing part and reports that part, even with no mask.
#[test]
fn polygon_mask_reports_rejected_parts_and_retains_valid_parts() {
let r = ramp();
let valid = vec![vec![Coord::new(0.0, 0.0), Coord::new(0.002, 0.0), Coord::new(0.002, 0.002), Coord::new(0.0, 0.002)]];
let mut invalid = valid.clone();
invalid.push(vec![Coord::new(f64::NAN, 0.0); 3]);
let polys = vec![invalid.clone(), valid.clone()];
let mut failures = Vec::new();
let mask = polygon_mask_with_failures(&r, &polys, &mut failures);
assert_eq!(mask, polygon_mask(&r, &[valid]));
assert_eq!(failures, vec![0]);
failures.clear();
let stats = zonal_statistics_with_failures(&r, &[0], &polys, &mut failures).unwrap();
assert_eq!(stats[0].count, 4);
assert_eq!(failures, vec![0]);
failures.clear();
assert!(polygon_mask_with_failures(&r, &[invalid], &mut failures).is_none());
assert_eq!(failures, vec![0]);
}

#[test]
fn terrain_math_and_classes() {
let r = ramp();
Expand Down
108 changes: 39 additions & 69 deletions crates/analysis/src/tools.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -587,11 +587,11 @@ pub fn find_field_case_insensitive(ds: &Dataset, name: &str) -> Result<usize, St
matches.push(idx);
}
}
match matches.len() {
0 => Err(format!("no field named '{name}'")),
1 => Ok(matches[0]),
match matches.as_slice() {
[] => Err(format!("no field named '{name}'")),
[only] => Ok(*only),
_ => {
let names: Vec<String> = matches.iter().map(|&i| ds.fields[i].name.clone()).collect();
let names: Vec<String> = matches.iter().filter_map(|&i| ds.fields.get(i).map(|f| f.name.clone())).collect();
Err(format!(
"field name '{name}' is ambiguous: case-insensitive collision between {}",
names.iter().map(|n| format!("'{n}'")).collect::<Vec<_>>().join(" and ")
Expand Down Expand Up @@ -1546,79 +1546,33 @@ impl StatOp {
}

/// Parses the `bands` parameter for zonal statistics.
/// Null and an empty string default to band 1, matching the workflow form.
pub fn zonal_bands(value: &Json, total_bands: usize) -> Result<Vec<usize>, String> {
let dummy = serde_json::json!({ "bands": value });
parse_bands(&dummy, total_bands)
}

fn parse_bands(p: &Json, total_bands: usize) -> Result<Vec<usize>, String> {
if total_bands == 0 {
return Err("raster has no bands".into());
}
match p.get("bands") {
None | Some(Json::Null) => Ok(vec![0]),
Some(Json::String(s)) => {
if s.trim().is_empty() {
return Err("bands parameter cannot be an empty string".into());
}
crate::rasterops::parse_bands_str(s, total_bands)
}
Some(Json::Number(n)) => {
let f = n.as_f64().ok_or_else(|| "invalid band number".to_string())?;
if !f.is_finite() {
return Err("band number must be a finite number".into());
}
if f <= 0.0 {
return Err("band numbers must be positive (1-based)".into());
}
if f.fract() != 0.0 {
return Err("band numbers must be integers, not decimals".into());
}
let b = f as usize;
if b > total_bands {
return Err(format!("band index {b} is invalid: raster has bands 1 to {total_bands}"));
}
Ok(vec![b - 1])
match value {
Json::Null => Ok(vec![0]),
Json::String(s) => crate::rasterops::parse_bands_str(s, total_bands),
Json::Number(n) => {
Ok(vec![crate::rasterops::band_from_f64(n.as_f64().ok_or_else(|| "invalid band number".to_string())?, total_bands)?])
}
Some(Json::Array(arr)) => {
Json::Array(arr) => {
if arr.is_empty() {
return Err("bands array cannot be empty".into());
}
let mut out = Vec::new();
for item in arr {
match item {
Json::Number(n) => {
let f = n.as_f64().ok_or_else(|| "invalid band number".to_string())?;
if !f.is_finite() {
return Err("band number must be a finite number".into());
}
if f <= 0.0 {
return Err("band numbers must be positive (1-based)".into());
}
if f.fract() != 0.0 {
return Err("band numbers must be integers, not decimals".into());
}
let b = f as usize;
if b > total_bands {
return Err(format!("band index {b} is invalid: raster has bands 1 to {total_bands}"));
}
let b_idx = b - 1;
if !out.contains(&b_idx) {
out.push(b_idx);
}
}
Json::String(s) => {
if s.trim().is_empty() {
return Err("band items cannot be empty strings".into());
}
let sub = crate::rasterops::parse_bands_str(s, total_bands)?;
for b_idx in sub {
if !out.contains(&b_idx) {
out.push(b_idx);
}
}
if !matches!(item, Json::Number(_) | Json::String(_)) {
return Err("band items must be numbers or strings".into());
}
if item.as_str().is_some_and(|s| s.trim().is_empty()) {
return Err("band items cannot be empty strings".into());
}
for band in zonal_bands(item, total_bands)? {
if !out.contains(&band) {
out.push(band);
}
_ => return Err("band items must be numbers or strings".into()),
}
}
Ok(out)
Expand All @@ -1627,6 +1581,10 @@ fn parse_bands(p: &Json, total_bands: usize) -> Result<Vec<usize>, String> {
}
}

fn parse_bands(p: &Json, total_bands: usize) -> Result<Vec<usize>, String> {
zonal_bands(p.get("bands").unwrap_or(&Json::Null), total_bands)
}

fn parse_stats(p: &Json) -> Result<Vec<StatOp>, String> {
match p.get("stats") {
None | Some(Json::Null) => Ok(StatOp::ALL.to_vec()),
Expand Down Expand Up @@ -1750,14 +1708,21 @@ fn zonal_stats(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {

let target_rows = rows(inp, p, poly_param, poly_ds);
let mut skipped_features = 0usize;
let mut skipped_parts = 0usize;
let mut failure_details = Vec::new();
for r_idx in target_rows {
let Some(f) = poly_ds.features.get(r_idx) else { continue };
let (stats_by_band, feature_valid) = match &f.geometry {
Some(Geometry::Polygon(polys)) if !polys.is_empty() => {
match crate::rasterops::transform_feature_polygons_to_lonlat(polys, poly_tr.as_ref()) {
Ok(polys_lonlat) => {
let st = crate::rasterops::zonal_statistics(r, &bands, &polys_lonlat)?;
let mut failed_parts = Vec::new();
let st = crate::rasterops::zonal_statistics_with_failures(r, &bands, &polys_lonlat, &mut failed_parts)?;
skipped_parts += failed_parts.len();
if !failed_parts.is_empty() {
failure_details
.push(format!("feature {}: raster projection failed for polygon part(s) {failed_parts:?}", f.id));
}
(st, true)
}
Err(e) => {
Expand Down Expand Up @@ -1798,6 +1763,11 @@ fn zonal_stats(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
let mut messages = vec![format!("Calculated zonal statistics for {n_features} feature(s)")];
if skipped_features > 0 {
messages.push(format!("{skipped_features} feature(s) had invalid geometry or failed coordinate transformation"));
}
if skipped_parts > 0 {
messages.push(format!("{skipped_parts} polygon part(s) could not be projected onto the raster grid; valid parts were retained"));
}
if !failure_details.is_empty() {
for detail in failure_details.iter().take(5) {
messages.push(detail.clone());
}
Expand Down Expand Up @@ -3140,9 +3110,9 @@ mod tests {
assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "stats": ["median"]})).is_err());
assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "stats": [123]})).is_err());

// 2. Invalid bands inputs
// 2. An empty band string defaults to band 1; malformed inputs fail.
assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "bands": []})).is_err());
assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "bands": ""})).is_err());
assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "bands": ""})).is_ok());
assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "bands": 1.5})).is_err());
assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "bands": 99})).is_err());
assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "bands": -1})).is_err());
Expand Down
Loading
Loading