diff --git a/crates/analysis/src/lib.rs b/crates/analysis/src/lib.rs
index 37ef54706..eda905cb1 100644
--- a/crates/analysis/src/lib.rs
+++ b/crates/analysis/src/lib.rs
@@ -24,7 +24,7 @@ pub mod topology;
 pub mod trend;
 pub mod warp;
 
-pub use tools::{TOOLS, holds_coordinates, simplify_geometry};
+pub use tools::{TOOLS, find_field_case_insensitive, holds_coordinates, simplify_geometry, zonal_bands};
 
 use rustgis_data::Dataset;
 use serde_json::Value;
diff --git a/crates/analysis/src/rasterops.rs b/crates/analysis/src/rasterops.rs
index ef44aae51..786d4b9cb 100644
--- a/crates/analysis/src/rasterops.rs
+++ b/crates/analysis/src/rasterops.rs
@@ -4,7 +4,7 @@
 
 use std::ops::Range;
 
-use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster};
+use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster, Transformer};
 use rustgis_geom::{Coord, Extent};
 
 /// The no-data value written by these operations.
@@ -523,6 +523,199 @@ 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 {
+    EmptyPolygon,
+    InvalidExterior(String),
+    InvalidHole { index: usize, reason: String },
+}
+
+impl std::fmt::Display for PolygonTransformError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            Self::EmptyPolygon => write!(f, "polygon has no rings"),
+            Self::InvalidExterior(r) => write!(f, "exterior ring is invalid: {r}"),
+            Self::InvalidHole { index, reason } => write!(f, "hole {index} is invalid: {reason}"),
+        }
+    }
+}
+
+impl std::error::Error for PolygonTransformError {}
+
+/// Transforms a single polygon's rings (exterior ring + optional holes) to WGS 84 lon/lat coordinates.
+///
+/// - The exterior ring must have at least 3 points and all points must project successfully.
+/// - Every interior ring (hole) must have at least 3 points and all points must project successfully.
+/// - If any hole fails or is invalid, the entire polygon is rejected with an error rather than
+///   silently discarding the hole or altering the polygon geometry.
+pub fn transform_polygon_to_lonlat(rings: &[Vec<Coord>], tr: Option<&Transformer>) -> Result<Vec<Vec<Coord>>, PolygonTransformError> {
+    let exterior = rings.first().ok_or(PolygonTransformError::EmptyPolygon)?;
+    let mut ext_lonlat = Vec::with_capacity(exterior.len());
+    for &c in exterior {
+        if !c.is_finite() {
+            return Err(PolygonTransformError::InvalidExterior("non-finite coordinate".into()));
+        }
+        let ll = match tr {
+            Some(t) => t.to_lonlat(c).ok_or_else(|| PolygonTransformError::InvalidExterior("coordinate transformation failed".into()))?,
+            None => c,
+        };
+        if !ll.is_finite() {
+            return Err(PolygonTransformError::InvalidExterior("transformed coordinate is non-finite".into()));
+        }
+        ext_lonlat.push(ll);
+    }
+    if ext_lonlat.len() < 3 {
+        return Err(PolygonTransformError::InvalidExterior("fewer than 3 vertices".into()));
+    }
+
+    let mut out_rings = vec![ext_lonlat];
+    for (idx, hole) in rings.iter().skip(1).enumerate() {
+        let hole_num = idx + 1;
+        let mut hole_lonlat = Vec::with_capacity(hole.len());
+        for &c in hole {
+            if !c.is_finite() {
+                return Err(PolygonTransformError::InvalidHole { index: hole_num, reason: "non-finite coordinate".into() });
+            }
+            let ll = match tr {
+                Some(t) => t.to_lonlat(c).ok_or_else(|| PolygonTransformError::InvalidHole {
+                    index: hole_num,
+                    reason: "coordinate transformation failed".into(),
+                })?,
+                None => c,
+            };
+            if !ll.is_finite() {
+                return Err(PolygonTransformError::InvalidHole { index: hole_num, reason: "transformed coordinate is non-finite".into() });
+            }
+            hole_lonlat.push(ll);
+        }
+        if hole_lonlat.len() < 3 {
+            return Err(PolygonTransformError::InvalidHole { index: hole_num, reason: "fewer than 3 vertices".into() });
+        }
+        out_rings.push(hole_lonlat);
+    }
+
+    Ok(out_rings)
+}
+
+/// Transforms and validates all polygon parts of a (multi-)polygon feature to WGS 84 lon/lat coordinates.
+pub fn transform_feature_polygons_to_lonlat(
+    polys: &[Vec<Vec<Coord>>],
+    tr: Option<&Transformer>,
+) -> Result<Vec<Vec<Vec<Coord>>>, PolygonTransformError> {
+    if polys.is_empty() {
+        return Err(PolygonTransformError::EmptyPolygon);
+    }
+    let mut out = Vec::with_capacity(polys.len());
+    for (part_idx, rings) in polys.iter().enumerate() {
+        let valid_poly = transform_polygon_to_lonlat(rings, tr).map_err(|e| match e {
+            PolygonTransformError::InvalidHole { index, reason } => {
+                PolygonTransformError::InvalidHole { index, reason: format!("part {part_idx}: {reason}") }
+            }
+            PolygonTransformError::InvalidExterior(reason) => PolygonTransformError::InvalidExterior(format!("part {part_idx}: {reason}")),
+            other => other,
+        })?;
+        out.push(valid_poly);
+    }
+    Ok(out)
+}
+
+/// Parses a band specification string into 0-based band indices.
+///
+/// Supported inputs:
+/// - `""` (empty string): defaults to band 1 (`vec![0]`)
+/// - `"all"`: all bands `0..total_bands`
+/// - comma, semicolon, or space-separated list of 1-based band numbers: `"1, 2, 3"`, `"[1, 2]"`, etc.
+///
+/// Validates that:
+/// - Every entry is a valid, finite positive integer.
+/// - Decimals (e.g. `"1.5"`), negative numbers (e.g. `"-1"`), zero (`"0"`), non-numeric (e.g. `"abc"`),
+///   or non-finite values (e.g. `"NaN"`, `"Infinity"`) are strictly rejected.
+/// - Every band number is within `1..=total_bands`.
+/// - Empty parts (such as `"1,,2"` or `","`) are rejected.
+pub fn parse_bands_str(s: &str, total_bands: usize) -> Result<Vec<usize>, String> {
+    if total_bands == 0 {
+        return Err("raster has no bands".into());
+    }
+    let txt = s.trim();
+    if txt.is_empty() {
+        return Ok(vec![0]);
+    }
+    if txt.eq_ignore_ascii_case("all") {
+        return Ok((0..total_bands).collect());
+    }
+    let clean = txt.trim_matches(|c| c == '[' || c == ']');
+    if clean.trim().is_empty() {
+        return Err("band selection cannot be empty brackets".into());
+    }
+    let mut out = Vec::new();
+    for part in clean.split([',', ';']) {
+        let part = part.trim().trim_matches('"').trim();
+        if part.is_empty() {
+            return Err("band selection contains an empty band entry".into());
+        }
+        for subpart in part.split_whitespace() {
+            let subpart = subpart.trim();
+            if subpart.is_empty() {
+                continue;
+            }
+            if subpart.eq_ignore_ascii_case("all") {
+                for b in 0..total_bands {
+                    if !out.contains(&b) {
+                        out.push(b);
+                    }
+                }
+                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;
+            if !out.contains(&b_idx) {
+                out.push(b_idx);
+            }
+        }
+    }
+    if out.is_empty() {
+        return Err("no valid band numbers found in input".into());
+    }
+    Ok(out)
+}
+
 /// 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.
@@ -533,14 +726,44 @@ pub fn polygon_mask(r: &Raster, polygons: &[Vec<Vec<Coord>>]) -> Option<(Range<u
     let tr = r.transformer();
     let projected: Vec<Vec<Vec<(f64, f64)>>> = polygons
         .iter()
-        .map(|rings| {
-            rings
-                .iter()
-                .map(|ring| ring.iter().filter_map(|c| tr.from_lonlat(*c)).map(|w| r.transform.to_pixel(w)).collect::<Vec<_>>())
-                .filter(|ring| ring.len() >= 3)
-                .collect::<Vec<_>>()
+        .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);
+            }
+            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)
         })
-        .filter(|rings| !rings.is_empty())
         .collect();
     let (mut c0, mut r0, mut c1, mut r1) = (f64::INFINITY, f64::INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY);
     for &(x, y) in projected.iter().flatten().flatten() {
@@ -601,6 +824,141 @@ pub fn polygon_mask(r: &Raster, polygons: &[Vec<Vec<Coord>>]) -> Option<(Range<u
     Some((cols, rows, trimmed))
 }
 
+/// Summary statistics of raster cells within a zone.
+#[derive(Clone, Debug, Default, PartialEq)]
+pub struct ZoneStats {
+    /// Number of valid raster cells inside the zone.
+    pub count: usize,
+    /// Arithmetic sum of valid raster cells inside the zone.
+    pub sum: Option<f64>,
+    /// Arithmetic mean of valid raster cells inside the zone.
+    pub mean: Option<f64>,
+    /// Minimum value among valid raster cells inside the zone.
+    pub min: Option<f64>,
+    /// Maximum value among valid raster cells inside the zone.
+    pub max: Option<f64>,
+    /// Population standard deviation of valid raster cells inside the zone: sqrt(sum((x - mean)^2) / n).
+    pub stddev: Option<f64>,
+}
+
+/// Numerically stable online accumulator using Welford's algorithm.
+#[derive(Clone, Debug)]
+pub struct Accumulator {
+    count: usize,
+    sum: f64,
+    min: f64,
+    max: f64,
+    mean: f64,
+    m2: f64,
+}
+
+impl Default for Accumulator {
+    /// Constructs an empty accumulator with `min = INFINITY` and `max = NEG_INFINITY`.
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl Accumulator {
+    /// Constructs a new empty accumulator with `min = INFINITY` and `max = NEG_INFINITY`.
+    pub fn new() -> Self {
+        Self { count: 0, sum: 0.0, min: f64::INFINITY, max: f64::NEG_INFINITY, mean: 0.0, m2: 0.0 }
+    }
+
+    /// Feeds a valid finite measurement into the accumulator, updating count, sum, extrema,
+    /// and running mean / variance using Welford's algorithm. Non-finite values are ignored.
+    pub fn add(&mut self, val: f64) {
+        if !val.is_finite() {
+            return;
+        }
+        self.count += 1;
+        self.sum += val;
+        if val < self.min {
+            self.min = val;
+        }
+        if val > self.max {
+            self.max = val;
+        }
+        let delta = val - self.mean;
+        self.mean += delta / (self.count as f64);
+        let delta2 = val - self.mean;
+        self.m2 += delta * delta2;
+    }
+
+    /// Number of valid finite samples accumulated so far.
+    pub fn count(&self) -> usize {
+        self.count
+    }
+
+    /// Computes and returns the accumulated [`ZoneStats`].
+    ///
+    /// If no valid values were added (`count == 0`), returns `ZoneStats::default()`
+    /// where `count` is 0 and all statistical metrics are `None`.
+    pub fn stats(&self) -> ZoneStats {
+        if self.count == 0 {
+            ZoneStats::default()
+        } else {
+            let variance = (self.m2 / (self.count as f64)).max(0.0);
+            ZoneStats {
+                count: self.count,
+                sum: Some(self.sum),
+                mean: Some(self.mean),
+                min: Some(self.min),
+                max: Some(self.max),
+                stddev: Some(variance.sqrt()),
+            }
+        }
+    }
+}
+
+/// Computes summary statistics (count, sum, mean, min, max, stddev) for each requested band
+/// over the raster cells whose centres fall inside `polygons` (exterior ring and holes in WGS 84 lon/lat).
+///
+/// Returns one `ZoneStats` per requested band. Cells with NoData, NaN, or non-finite values are
+/// 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> {
+    r.require_coordinates()?;
+    for &b in bands {
+        if b >= r.band_count() {
+            return Err(format!("band index {} is out of range: {} has {} bands", b + 1, r.name, r.band_count()));
+        }
+    }
+    if bands.is_empty() {
+        return Ok(Vec::new());
+    }
+
+    let Some((cols, rows, mask)) = polygon_mask(r, polygons) else {
+        return Ok(vec![ZoneStats::default(); bands.len()]);
+    };
+
+    let w = (cols.end - cols.start) as usize;
+    let mut accs: Vec<Accumulator> = (0..bands.len()).map(|_| Accumulator::new()).collect();
+
+    for (ri, row) in rows.enumerate() {
+        let line_offset = ri * w;
+        let line = match mask.get(line_offset..line_offset + w) {
+            Some(l) => l,
+            None => continue,
+        };
+        for (ci, &inside) in line.iter().enumerate() {
+            if !inside {
+                continue;
+            }
+            let col = cols.start + ci as u32;
+            for (acc_idx, &band) in bands.iter().enumerate() {
+                if let Some(v) = r.value(band, col as i64, row as i64)
+                    && let Some(acc) = accs.get_mut(acc_idx)
+                {
+                    acc.add(f64::from(v));
+                }
+            }
+        }
+    }
+
+    Ok(accs.into_iter().map(|acc| acc.stats()).collect())
+}
+
 /// A tiny arithmetic expression over band values.
 enum Expr {
     Num(f64),
@@ -785,4 +1143,329 @@ mod tests {
         assert_eq!(p.len(), 20);
         assert!(p.last().and_then(|x| x.2).unwrap_or(0.0) > p[0].2.unwrap_or(0.0) + 150.0);
     }
+
+    #[test]
+    fn accumulator_exact_math_and_numerical_stability() {
+        // 1. Exact math on [10.0, 20.0, 30.0, 40.0]
+        let mut acc = Accumulator::new();
+        for &x in &[10.0, 20.0, 30.0, 40.0] {
+            acc.add(x);
+        }
+        let st = acc.stats();
+        assert_eq!(st.count, 4);
+        assert_eq!(st.sum, Some(100.0));
+        assert_eq!(st.mean, Some(25.0));
+        assert_eq!(st.min, Some(10.0));
+        assert_eq!(st.max, Some(40.0));
+        let expected_stddev = 125.0f64.sqrt(); // sqrt(500 / 4)
+        assert!((st.stddev.unwrap() - expected_stddev).abs() < 1e-12);
+
+        // 2. Empty accumulator
+        let empty_st = Accumulator::new().stats();
+        assert_eq!(empty_st.count, 0);
+        assert_eq!(empty_st.sum, None);
+        assert_eq!(empty_st.mean, None);
+        assert_eq!(empty_st.min, None);
+        assert_eq!(empty_st.max, None);
+        assert_eq!(empty_st.stddev, None);
+
+        // 3. Single value
+        let mut single_acc = Accumulator::new();
+        single_acc.add(42.5);
+        let single_st = single_acc.stats();
+        assert_eq!(single_st.count, 1);
+        assert_eq!(single_st.sum, Some(42.5));
+        assert_eq!(single_st.mean, Some(42.5));
+        assert_eq!(single_st.min, Some(42.5));
+        assert_eq!(single_st.max, Some(42.5));
+        assert_eq!(single_st.stddev, Some(0.0));
+
+        // 4. Constant values
+        let mut const_acc = Accumulator::new();
+        for _ in 0..100 {
+            const_acc.add(7.0);
+        }
+        let const_st = const_acc.stats();
+        assert_eq!(const_st.count, 100);
+        assert_eq!(const_st.mean, Some(7.0));
+        assert_eq!(const_st.stddev, Some(0.0));
+
+        // 5. Large magnitude values (Welford's algorithm avoids catastrophic cancellation)
+        let base = 1e12;
+        let mut large_acc = Accumulator::new();
+        for &x in &[base + 1.0, base + 2.0, base + 3.0, base + 4.0] {
+            large_acc.add(x);
+        }
+        let large_st = large_acc.stats();
+        assert_eq!(large_st.count, 4);
+        assert_eq!(large_st.mean, Some(base + 2.5));
+        let std_large = large_st.stddev.unwrap();
+        // Variance of [1, 2, 3, 4] is 1.25, stddev is sqrt(1.25) = 1.118033988749895
+        assert!((std_large - 1.25f64.sqrt()).abs() < 1e-10, "large magnitude precision: {std_large}");
+
+        // 6. Non-finite values (NaN and inf) are rejected
+        let mut defensive_acc = Accumulator::new();
+        defensive_acc.add(10.0);
+        defensive_acc.add(f64::NAN);
+        defensive_acc.add(f64::INFINITY);
+        defensive_acc.add(f64::NEG_INFINITY);
+        defensive_acc.add(20.0);
+        let def_st = defensive_acc.stats();
+        assert_eq!(def_st.count, 2);
+        assert_eq!(def_st.sum, Some(30.0));
+        assert_eq!(def_st.mean, Some(15.0));
+    }
+
+    #[test]
+    fn accumulator_default_and_new_consistency_and_extrema() {
+        // default() and new() produce identical initial empty state
+        let d = Accumulator::default();
+        let n = Accumulator::new();
+        assert_eq!(d.count(), 0);
+        assert_eq!(n.count(), 0);
+        assert_eq!(d.stats(), ZoneStats::default());
+        assert_eq!(n.stats(), ZoneStats::default());
+
+        // Negative-only values: min must be -30, max must be -10 (not 0.0)
+        let mut neg_acc = Accumulator::default();
+        neg_acc.add(-30.0);
+        neg_acc.add(-20.0);
+        neg_acc.add(-10.0);
+        let neg_st = neg_acc.stats();
+        assert_eq!(neg_st.count, 3);
+        assert_eq!(neg_st.min, Some(-30.0));
+        assert_eq!(neg_st.max, Some(-10.0));
+        assert_eq!(neg_st.sum, Some(-60.0));
+        assert_eq!(neg_st.mean, Some(-20.0));
+
+        // Positive-only values: min must be 10 (not 0.0), max must be 30
+        let mut pos_acc = Accumulator::default();
+        pos_acc.add(10.0);
+        pos_acc.add(20.0);
+        pos_acc.add(30.0);
+        let pos_st = pos_acc.stats();
+        assert_eq!(pos_st.count, 3);
+        assert_eq!(pos_st.min, Some(10.0));
+        assert_eq!(pos_st.max, Some(30.0));
+
+        // Mixed positive and negative values
+        let mut mixed_acc = Accumulator::default();
+        mixed_acc.add(-5.0);
+        mixed_acc.add(5.0);
+        let mixed_st = mixed_acc.stats();
+        assert_eq!(mixed_st.count, 2);
+        assert_eq!(mixed_st.min, Some(-5.0));
+        assert_eq!(mixed_st.max, Some(5.0));
+        assert_eq!(mixed_st.sum, Some(0.0));
+        assert_eq!(mixed_st.mean, Some(0.0));
+    }
+
+    #[test]
+    fn zonal_statistics_known_cells_nodata_and_holes() {
+        let r = ramp(); // 20x20 raster, band 0: (i % 20) * 10, band 1: i / 20.
+        let sq = |x0: f64, y0: f64, x1: f64, y1: f64| {
+            vec![Coord::new(x0, y0), Coord::new(x1, y0), Coord::new(x1, y1), Coord::new(x0, y1), Coord::new(x0, y0)]
+        };
+
+        // 1. Single known cell: col 5, row 5.
+        // origin_x = 0.0, origin_y = 0.02, pixel_width = 0.001, pixel_height = -0.001
+        // col 5 x in [0.005, 0.006], centre at 0.0055.
+        // row 5 y in [0.014, 0.015], centre at 0.0145.
+        let cell_poly = vec![vec![sq(0.005, 0.014, 0.006, 0.015)]];
+        let stats = zonal_statistics(&r, &[0, 1], &cell_poly).unwrap();
+        assert_eq!(stats.len(), 2);
+        assert_eq!(stats[0].count, 1);
+        assert_eq!(stats[0].sum, Some(50.0)); // col 5 * 10 = 50
+        assert_eq!(stats[0].mean, Some(50.0));
+        assert_eq!(stats[0].min, Some(50.0));
+        assert_eq!(stats[0].max, Some(50.0));
+        assert_eq!(stats[0].stddev, Some(0.0));
+
+        assert_eq!(stats[1].count, 1);
+        assert_eq!(stats[1].sum, Some(5.0)); // row 5 = 5
+        assert_eq!(stats[1].mean, Some(5.0));
+        assert_eq!(stats[1].stddev, Some(0.0));
+
+        // 2. Polygon outside raster: empty zone result
+        let outside_poly = vec![vec![sq(1.0, 1.0, 2.0, 2.0)]];
+        let outside_stats = zonal_statistics(&r, &[0], &outside_poly).unwrap();
+        assert_eq!(outside_stats.len(), 1);
+        assert_eq!(outside_stats[0], ZoneStats::default());
+
+        // 3. Polygon with holes
+        let poly_with_hole = vec![vec![sq(0.005, 0.005, 0.015, 0.015), sq(0.009, 0.009, 0.011, 0.011)]];
+        let hole_stats = zonal_statistics(&r, &[0], &poly_with_hole).unwrap();
+        // 100 cells - 4 hole cells = 96 cells
+        assert_eq!(hole_stats[0].count, 96);
+        assert!(hole_stats[0].sum.is_some());
+        assert!(hole_stats[0].mean.is_some());
+
+        // 4. MultiPolygon
+        let multipoly = vec![
+            vec![sq(0.001, 0.018, 0.003, 0.019)], // 2 cells
+            vec![sq(0.005, 0.018, 0.007, 0.019)], // 2 cells
+        ];
+        let multi_stats = zonal_statistics(&r, &[0], &multipoly).unwrap();
+        assert_eq!(multi_stats[0].count, 4);
+
+        // 5. NoData handling: create a raster with NoData values
+        let (w, h) = (4u32, 4u32);
+        // cells:
+        // [ 10,   20, -9999,   40 ]
+        // [ 50,  NaN,   70,   80 ]
+        // [ -9999, -9999, -9999, -9999 ]
+        // [ 130, 140, 150, 160 ]
+        let data: Vec<f32> =
+            vec![10.0, 20.0, -9999.0, 40.0, 50.0, f32::NAN, 70.0, 80.0, -9999.0, -9999.0, -9999.0, -9999.0, 130.0, 140.0, 150.0, 160.0];
+        let gt = GeoTransform { origin_x: 0.0, origin_y: 4.0, pixel_width: 1.0, pixel_height: -1.0 };
+        let r_nodata = Raster::new("nd", w, h, vec![data], gt, Crs::LonLat, Some(-9999.0), ColorHint::Gray, "Float32").unwrap();
+
+        // Zone covering entire row 0: cells are 10, 20, -9999 (NoData), 40
+        // Valid values: 10, 20, 40 (count = 3, sum = 70, mean = 70/3, min = 10, max = 40)
+        let row0_poly = vec![vec![sq(0.0, 3.0, 4.0, 4.0)]];
+        let r0_stats = zonal_statistics(&r_nodata, &[0], &row0_poly).unwrap();
+        assert_eq!(r0_stats[0].count, 3);
+        assert_eq!(r0_stats[0].sum, Some(70.0));
+        assert_eq!(r0_stats[0].min, Some(10.0));
+        assert_eq!(r0_stats[0].max, Some(40.0));
+        assert!((r0_stats[0].mean.unwrap() - 70.0 / 3.0).abs() < 1e-10);
+
+        // Zone covering row 2 where ALL cells are NoData (-9999)
+        let row2_poly = vec![vec![sq(0.0, 1.0, 4.0, 2.0)]];
+        let r2_stats = zonal_statistics(&r_nodata, &[0], &row2_poly).unwrap();
+        assert_eq!(r2_stats[0].count, 0);
+        assert_eq!(r2_stats[0].sum, None);
+        assert_eq!(r2_stats[0].mean, None);
+        assert_eq!(r2_stats[0].min, None);
+        assert_eq!(r2_stats[0].max, None);
+        assert_eq!(r2_stats[0].stddev, None);
+
+        // 6. Invalid band index returns Err
+        assert!(zonal_statistics(&r, &[99], &cell_poly).is_err());
+
+        // 7. Invalid exterior ring (< 3 points) with surviving hole does NOT promote hole to exterior
+        let bad_ext_with_hole = vec![vec![
+            vec![Coord::new(0.005, 0.005), Coord::new(0.015, 0.015)], // only 2 points!
+            sq(0.009, 0.009, 0.011, 0.011),                           // 5 points (hole)
+        ]];
+        let bad_ext_stats = zonal_statistics(&r, &[0], &bad_ext_with_hole).unwrap();
+        assert_eq!(bad_ext_stats[0], ZoneStats::default());
+        assert!(polygon_mask(&r, &bad_ext_with_hole).is_none());
+
+        // 8. Exterior ring with non-finite/unprojectable coordinate is rejected
+        let nan_ext_with_hole = vec![vec![
+            vec![Coord::new(f64::NAN, 0.005), Coord::new(0.015, 0.005), Coord::new(0.015, 0.015), Coord::new(f64::NAN, 0.005)],
+            sq(0.009, 0.009, 0.011, 0.011),
+        ]];
+        let nan_ext_stats = zonal_statistics(&r, &[0], &nan_ext_with_hole).unwrap();
+        assert_eq!(nan_ext_stats[0], ZoneStats::default());
+        assert!(polygon_mask(&r, &nan_ext_with_hole).is_none());
+
+        // 9. Polygon with defective hole (< 3 vertices or non-finite) is rejected rather than silently discarding the hole
+        let bad_hole_poly = vec![vec![
+            sq(0.005, 0.005, 0.015, 0.015),                           // valid exterior (5 points)
+            vec![Coord::new(0.009, 0.009), Coord::new(0.011, 0.011)], // invalid hole: only 2 points!
+        ]];
+        assert!(polygon_mask(&r, &bad_hole_poly).is_none(), "polygon with defective hole must be rejected");
+        let bad_hole_stats = zonal_statistics(&r, &[0], &bad_hole_poly).unwrap();
+        assert_eq!(bad_hole_stats[0], ZoneStats::default());
+
+        let nan_hole_poly = vec![vec![
+            sq(0.005, 0.005, 0.015, 0.015),
+            vec![Coord::new(f64::NAN, 0.009), Coord::new(0.011, 0.009), Coord::new(0.011, 0.011), Coord::new(f64::NAN, 0.009)],
+        ]];
+        assert!(polygon_mask(&r, &nan_hole_poly).is_none(), "polygon with non-finite hole must be rejected");
+    }
+
+    #[test]
+    fn polygon_transformation_and_hole_validation_semantics() {
+        let sq = |x0: f64, y0: f64, x1: f64, y1: f64| {
+            vec![Coord::new(x0, y0), Coord::new(x1, y0), Coord::new(x1, y1), Coord::new(x0, y1), Coord::new(x0, y0)]
+        };
+        let ext = sq(0.0, 0.0, 10.0, 10.0);
+        let hole1 = sq(2.0, 2.0, 4.0, 4.0);
+        let hole2 = sq(6.0, 6.0, 8.0, 8.0);
+
+        // Valid exterior and valid holes
+        let valid_poly = vec![ext.clone(), hole1.clone(), hole2.clone()];
+        let res = transform_polygon_to_lonlat(&valid_poly, None).unwrap();
+        assert_eq!(res.len(), 3);
+        assert_eq!(res[0].len(), ext.len());
+        assert_eq!(res[1].len(), hole1.len());
+        assert_eq!(res[2].len(), hole2.len());
+
+        // Empty polygon
+        let empty_poly: Vec<Vec<Coord>> = Vec::new();
+        assert_eq!(transform_polygon_to_lonlat(&empty_poly, None), Err(PolygonTransformError::EmptyPolygon));
+
+        // Invalid exterior: fewer than 3 vertices
+        let bad_ext = vec![vec![Coord::new(0.0, 0.0), Coord::new(1.0, 1.0)]];
+        assert!(matches!(transform_polygon_to_lonlat(&bad_ext, None), Err(PolygonTransformError::InvalidExterior(_))));
+
+        // Invalid exterior: non-finite coordinate
+        let nan_ext = vec![vec![Coord::new(0.0, 0.0), Coord::new(f64::NAN, 1.0), Coord::new(2.0, 2.0), Coord::new(0.0, 0.0)]];
+        assert!(matches!(transform_polygon_to_lonlat(&nan_ext, None), Err(PolygonTransformError::InvalidExterior(_))));
+
+        // Invalid hole 1: fewer than 3 vertices
+        let bad_hole = vec![ext.clone(), vec![Coord::new(2.0, 2.0), Coord::new(3.0, 3.0)]];
+        assert_eq!(
+            transform_polygon_to_lonlat(&bad_hole, None),
+            Err(PolygonTransformError::InvalidHole { index: 1, reason: "fewer than 3 vertices".into() })
+        );
+
+        // Invalid hole 2: non-finite coordinate
+        let nan_hole = vec![ext.clone(), hole1.clone(), vec![Coord::new(6.0, 6.0), Coord::new(f64::INFINITY, 8.0), Coord::new(8.0, 8.0)]];
+        assert!(matches!(transform_polygon_to_lonlat(&nan_hole, None), Err(PolygonTransformError::InvalidHole { index: 2, .. })));
+
+        // Multipart features
+        let multi_valid = vec![vec![ext.clone(), hole1.clone()], vec![sq(20.0, 20.0, 30.0, 30.0)]];
+        let multi_res = transform_feature_polygons_to_lonlat(&multi_valid, None).unwrap();
+        assert_eq!(multi_res.len(), 2);
+        assert_eq!(multi_res[0].len(), 2);
+        assert_eq!(multi_res[1].len(), 1);
+
+        // Multipart with failing part
+        let multi_bad = vec![vec![ext.clone()], vec![vec![Coord::new(1.0, 1.0)]]];
+        assert!(transform_feature_polygons_to_lonlat(&multi_bad, None).is_err());
+    }
+
+    #[test]
+    fn parse_bands_str_validation_and_boundaries() {
+        let total = 3usize;
+
+        // Valid formats
+        assert_eq!(parse_bands_str("", total).unwrap(), vec![0]);
+        assert_eq!(parse_bands_str("   ", total).unwrap(), vec![0]);
+        assert_eq!(parse_bands_str("all", total).unwrap(), vec![0, 1, 2]);
+        assert_eq!(parse_bands_str("ALL", total).unwrap(), vec![0, 1, 2]);
+        assert_eq!(parse_bands_str("1", total).unwrap(), vec![0]);
+        assert_eq!(parse_bands_str("3", total).unwrap(), vec![2]); // boundary: total_bands
+        assert_eq!(parse_bands_str("1, 2", total).unwrap(), vec![0, 1]);
+        assert_eq!(parse_bands_str("[1, 2, 3]", total).unwrap(), vec![0, 1, 2]);
+        assert_eq!(parse_bands_str("1; 3", total).unwrap(), vec![0, 2]);
+        assert_eq!(parse_bands_str(" 1   2 ", total).unwrap(), vec![0, 1]);
+
+        // De-duplication
+        assert_eq!(parse_bands_str("1, 1, 2, 2, 3", total).unwrap(), vec![0, 1, 2]);
+
+        // Boundary violations and negative
+        assert!(parse_bands_str("0", total).is_err(), "0 is invalid (1-based)");
+        assert!(parse_bands_str("-1", total).is_err(), "negative band is invalid");
+        assert!(parse_bands_str("4", total).is_err(), "out of range band is invalid");
+        assert!(parse_bands_str("999", total).is_err(), "out of range band is invalid");
+
+        // Decimals and non-finite
+        assert!(parse_bands_str("1.5", total).is_err(), "decimal band is invalid");
+        assert_eq!(parse_bands_str("1.0", total).unwrap(), vec![0]); // 1.0 has fract == 0.0
+        assert!(parse_bands_str("NaN", total).is_err(), "NaN is invalid");
+        assert!(parse_bands_str("Infinity", total).is_err(), "Infinity is invalid");
+
+        // Malformed strings
+        assert!(parse_bands_str("abc", total).is_err(), "non-numeric string is invalid");
+        assert!(parse_bands_str("[]", total).is_err(), "empty brackets are invalid");
+        assert!(parse_bands_str("[   ]", total).is_err(), "whitespace brackets are invalid");
+        assert!(parse_bands_str("1,,2", total).is_err(), "empty item between commas is invalid");
+        assert!(parse_bands_str(",", total).is_err(), "standalone comma is invalid");
+        assert!(parse_bands_str(";", total).is_err(), "standalone semicolon is invalid");
+    }
 }
diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs
index 109aaf622..feeaf7054 100644
--- a/crates/analysis/src/tools.rs
+++ b/crates/analysis/src/tools.rs
@@ -457,6 +457,56 @@ pub static TOOLS: &[ToolSpec] = &[
         ],
         run: kriging_interpolation,
     },
+    ToolSpec {
+        id: "zonalStats",
+        name: "Zonal Statistics",
+        toolbox: "Analysis › Raster",
+        summary: "Calculates raster summary statistics for each polygon feature.",
+        params: &[
+            P {
+                name: "raster",
+                label: "Input Raster",
+                kind: K::RasterLayer,
+                required: true,
+                default: "",
+                help: "The raster containing pixel values to summarize.",
+            },
+            P {
+                name: "polygons",
+                label: "Zone Polygons",
+                kind: K::PolygonLayer,
+                required: true,
+                default: "",
+                help: "Polygon layer defining the zones.",
+            },
+            P {
+                name: "bands",
+                label: "Bands",
+                kind: K::Text,
+                required: false,
+                default: "\"1\"",
+                help: "1-based band number(s) (e.g. 1 or [1, 2]); defaults to band 1.",
+            },
+            P {
+                name: "stats",
+                label: "Statistics",
+                kind: K::Text,
+                required: false,
+                default: "\"All\"",
+                help: "Statistical operations: All, Mean, Min, Max, Sum, Count, StdDev (comma-separated or array); defaults to All.",
+            },
+            P {
+                name: "prefix",
+                label: "Field Prefix",
+                kind: K::Text,
+                required: false,
+                default: "",
+                help: "Prefix for generated fields; defaults to raster/band name.",
+            },
+            P { name: "output", label: "Output Name", kind: K::Text, required: false, default: "", help: "Output polygon layer name." },
+        ],
+        run: zonal_stats,
+    },
 ];
 
 // ---------------------------------------------------------------------------------------------
@@ -523,6 +573,33 @@ fn copy_schema(ds: &Dataset, name: String, g: GeometryType) -> Dataset {
     Dataset::new(name, g, ds.fields.clone())
 }
 
+/// Finds the index of a field by case-insensitive name, explicitly detecting collisions.
+///
+/// If multiple fields in `ds.fields` match `name` case-insensitively, returns an error
+/// detailing the collision rather than silently selecting an arbitrary field.
+pub fn find_field_case_insensitive(ds: &Dataset, name: &str) -> Result<usize, String> {
+    if let Some(idx) = ds.fields.iter().position(|f| f.name == name) {
+        return Ok(idx);
+    }
+    let mut matches = Vec::new();
+    for (idx, f) in ds.fields.iter().enumerate() {
+        if f.name.eq_ignore_ascii_case(name) {
+            matches.push(idx);
+        }
+    }
+    match matches.len() {
+        0 => Err(format!("no field named '{name}'")),
+        1 => Ok(matches[0]),
+        _ => {
+            let names: Vec<String> = matches.iter().map(|&i| ds.fields[i].name.clone()).collect();
+            Err(format!(
+                "field name '{name}' is ambiguous: case-insensitive collision between {}",
+                names.iter().map(|n| format!("'{n}'")).collect::<Vec<_>>().join(" and ")
+            ))
+        }
+    }
+}
+
 /// Union of polygons; rings are oriented first because the overlay infers its fill rule from
 /// the first ring's winding, and input files mix orientations.
 fn union_all<'a>(polys: impl IntoIterator<Item = &'a geo::Polygon<f64>>) -> MultiPolygon<f64> {
@@ -632,7 +709,7 @@ fn dissolve(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
         return Err("Dissolve needs a polygon layer".into());
     }
     let field = s(p, "field").trim();
-    let fi = if field.is_empty() { None } else { Some(ds.field_index(field).ok_or_else(|| format!("no field named '{field}'"))?) };
+    let fi = if field.is_empty() { None } else { Some(find_field_case_insensitive(ds, field)?) };
     let mut groups: Vec<(Value, Vec<MultiPolygon<f64>>)> = Vec::new();
     for r in rows(inp, p, "input", ds) {
         let Some(mp) = ds.features.get(r).and_then(|f| f.geometry.as_ref()).and_then(to_multipolygon) else { continue };
@@ -745,7 +822,7 @@ fn xy_table_to_point(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
     let ds = input(inp, p, "input")?;
     let field = |k: &str| {
         let name = s(p, k).trim();
-        ds.field_index(name).ok_or_else(|| format!("no field named '{name}'"))
+        find_field_case_insensitive(ds, name)
     };
     let (xi, yi) = (field("x")?, field("y")?);
     let crs = match s(p, "crs").trim() {
@@ -863,7 +940,7 @@ fn select_by_attributes(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String
 fn idw_interpolation(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
     let ds = input(inp, p, "input")?;
     let field = s(p, "field");
-    let fi = ds.field_index(field).ok_or_else(|| format!("no field named '{field}'"))?;
+    let fi = find_field_case_insensitive(ds, field)?;
     let mut pts: Vec<(Coord, f64)> = Vec::new();
     for r in rows(inp, p, "input", ds) {
         let Some(g) = ds.features.get(r).and_then(|f| f.geometry.as_ref()) else { continue };
@@ -935,7 +1012,7 @@ fn kriging_interpolation(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, Strin
         return Err("kriging interpolation expects coordinates in degrees (WGS 84 / EPSG:4326)".into());
     }
     let field = s(p, "field");
-    let fi = ds.field_index(field).ok_or_else(|| format!("no field named '{field}'"))?;
+    let fi = find_field_case_insensitive(ds, field)?;
     let mut raw_pts: Vec<(Coord, f64)> = Vec::new();
     for r in rows(inp, p, "input", ds) {
         let Some(g) = ds.features.get(r).and_then(|f| f.geometry.as_ref()) else { continue };
@@ -1155,9 +1232,9 @@ fn calculate_geometry(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String>
 fn summary_statistics(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
     let ds = input(inp, p, "input")?;
     let field = s(p, "field");
-    let fi = ds.field_index(field).ok_or_else(|| format!("no field named '{field}'"))?;
+    let fi = find_field_case_insensitive(ds, field)?;
     let case = s(p, "case").trim();
-    let ci = if case.is_empty() { None } else { Some(ds.field_index(case).ok_or_else(|| format!("no field named '{case}'"))?) };
+    let ci = if case.is_empty() { None } else { Some(find_field_case_insensitive(ds, case)?) };
     let mut groups: Vec<(Value, Vec<f64>)> = Vec::new();
     for r in rows(inp, p, "input", ds) {
         let key = ci.map(|c| ds.value(r, c).clone()).unwrap_or(Value::Null);
@@ -1421,9 +1498,320 @@ fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
     result(out, format!("Spatial Join: {written_features} features joined ({total_matches} matches)"))
 }
 
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+enum StatOp {
+    Mean,
+    Min,
+    Max,
+    Sum,
+    Count,
+    StdDev,
+}
+
+impl StatOp {
+    pub const ALL: [StatOp; 6] = [StatOp::Mean, StatOp::Min, StatOp::Max, StatOp::Sum, StatOp::Count, StatOp::StdDev];
+
+    fn name(self) -> &'static str {
+        match self {
+            StatOp::Mean => "mean",
+            StatOp::Min => "min",
+            StatOp::Max => "max",
+            StatOp::Sum => "sum",
+            StatOp::Count => "count",
+            StatOp::StdDev => "stddev",
+        }
+    }
+
+    fn field_type(self) -> FieldType {
+        match self {
+            StatOp::Count => FieldType::Integer,
+            _ => FieldType::Double,
+        }
+    }
+
+    fn parse(raw: &str) -> Result<Vec<StatOp>, String> {
+        let clean = raw.trim().to_ascii_lowercase();
+        match clean.as_str() {
+            "all" => Ok(Self::ALL.to_vec()),
+            "mean" => Ok(vec![StatOp::Mean]),
+            "min" => Ok(vec![StatOp::Min]),
+            "max" => Ok(vec![StatOp::Max]),
+            "sum" => Ok(vec![StatOp::Sum]),
+            "count" => Ok(vec![StatOp::Count]),
+            "stddev" | "std" | "std_dev" => Ok(vec![StatOp::StdDev]),
+            "" => Err("empty statistic name: choose from mean, min, max, sum, count, stddev, or all".into()),
+            other => Err(format!("unknown statistic '{other}': choose from mean, min, max, sum, count, stddev, or all")),
+        }
+    }
+}
+
+/// Parses the `bands` parameter for zonal statistics.
+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])
+        }
+        Some(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);
+                            }
+                        }
+                    }
+                    _ => return Err("band items must be numbers or strings".into()),
+                }
+            }
+            Ok(out)
+        }
+        _ => Err("bands must be a number, string, or array of numbers/strings".into()),
+    }
+}
+
+fn parse_stats(p: &Json) -> Result<Vec<StatOp>, String> {
+    match p.get("stats") {
+        None | Some(Json::Null) => Ok(StatOp::ALL.to_vec()),
+        Some(Json::Array(arr)) => {
+            if arr.is_empty() {
+                return Err("stats array cannot be empty".into());
+            }
+            let mut out = Vec::new();
+            for item in arr {
+                let txt = item.as_str().ok_or_else(|| "statistic names must be strings".to_string())?;
+                if txt.trim().is_empty() {
+                    return Err("statistic name cannot be empty".into());
+                }
+                let ops = StatOp::parse(txt)?;
+                for op in ops {
+                    if !out.contains(&op) {
+                        out.push(op);
+                    }
+                }
+            }
+            Ok(out)
+        }
+        Some(Json::String(txt)) => {
+            let trimmed = txt.trim();
+            if trimmed.is_empty() {
+                return Err("stats parameter cannot be empty string".into());
+            }
+            let clean = trimmed.trim_matches(|c| c == '[' || c == ']');
+            if clean.trim().is_empty() {
+                return Err("stats parameter cannot be empty brackets".into());
+            }
+            let mut out = Vec::new();
+            for part in clean.split([',', ';']) {
+                let part = part.trim().trim_matches('"').trim();
+                if part.is_empty() {
+                    return Err("stats selection contains an empty statistic entry".into());
+                }
+                let ops = StatOp::parse(part)?;
+                for op in ops {
+                    if !out.contains(&op) {
+                        out.push(op);
+                    }
+                }
+            }
+            Ok(out)
+        }
+        _ => Err("stats parameter must be a string or array of strings".into()),
+    }
+}
+
+fn zonal_stats(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
+    let r_name = {
+        let n = s(p, "raster");
+        if n.is_empty() { s(p, "layer") } else { n }
+    };
+    if r_name.is_empty() {
+        return Err("parameter 'raster' is required".into());
+    }
+    let r_input = inp.raster(r_name).ok_or_else(|| format!("no raster layer named '{r_name}'"))?;
+    let r = r_input.raster;
+    r.require_coordinates()?;
+
+    let poly_param = if !s(p, "polygons").is_empty() { "polygons" } else { "input" };
+    let poly_ds = input(inp, p, poly_param)?;
+    if poly_ds.geometry_type != GeometryType::Polygon {
+        return Err(format!(
+            "layer '{}' is {} data; zonal statistics needs polygon geometries",
+            poly_ds.name,
+            poly_ds.geometry_type.label()
+        ));
+    }
+    let poly_tr = match poly_ds.crs {
+        Some(crs) if crs != rustgis_data::raster::Crs::LonLat && crs != rustgis_data::raster::Crs::Epsg(4326) => {
+            Some(crs.transformer().map_err(|e| format!("polygon layer coordinate system isn't supported: {e}"))?)
+        }
+        _ => None,
+    };
+
+    let bands = parse_bands(p, r.band_count())?;
+    let stats = parse_stats(p)?;
+
+    let mut out = copy_schema(poly_ds, out_name(p, format!("{}_ZonalStats", poly_ds.name)), poly_ds.geometry_type);
+    out.source = None;
+    out.crs = poly_ds.crs;
+
+    let user_prefix = s(p, "prefix").trim();
+    let multi_bands = bands.len() > 1;
+    let mut out_fields = Vec::new();
+    let mut new_field_names: Vec<String> = Vec::new();
+
+    for (b_idx, &band) in bands.iter().enumerate() {
+        let band_num = band + 1;
+        let band_name = r.bands.get(band).map(|b| b.name.as_str()).unwrap_or("");
+        let has_named_band = !band_name.is_empty() && !band_name.starts_with("Band ") && band_name != "Gray";
+        let prefix = if !user_prefix.is_empty() {
+            if multi_bands { format!("{user_prefix}_b{band_num}") } else { user_prefix.to_string() }
+        } else if multi_bands {
+            if has_named_band { format!("{}_{}", r.name, band_name) } else { format!("{}_b{band_num}", r.name) }
+        } else if has_named_band {
+            band_name.to_string()
+        } else {
+            r.name.to_string()
+        };
+
+        for &stat in &stats {
+            let base_name = format!("{prefix}_{}", stat.name());
+            let mut final_name = base_name.clone();
+            let mut counter = 1usize;
+            while out.field_index(&final_name).is_some() || new_field_names.iter().any(|n| n.eq_ignore_ascii_case(&final_name)) {
+                counter += 1;
+                final_name = format!("{base_name}_{counter}");
+            }
+            new_field_names.push(final_name.clone());
+            out_fields.push((b_idx, band, stat, final_name, stat.field_type()));
+        }
+    }
+
+    for (_, _, _, name, ty) in &out_fields {
+        out.fields.push(Field::new(name.clone(), *ty));
+    }
+
+    let target_rows = rows(inp, p, poly_param, poly_ds);
+    let mut skipped_features = 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)?;
+                        (st, true)
+                    }
+                    Err(e) => {
+                        failure_details.push(format!("feature {}: {e}", f.id));
+                        (vec![crate::rasterops::ZoneStats::default(); bands.len()], false)
+                    }
+                }
+            }
+            _ => {
+                failure_details.push(format!("feature {}: missing or empty geometry", f.id));
+                (vec![crate::rasterops::ZoneStats::default(); bands.len()], false)
+            }
+        };
+
+        if !feature_valid {
+            skipped_features += 1;
+        }
+
+        let mut attrs = f.attrs.clone();
+        attrs.resize(poly_ds.fields.len(), Value::Null);
+        for (b_idx, _, stat, _, _) in &out_fields {
+            let b_stats = stats_by_band.get(*b_idx);
+            let val = match (b_stats, stat) {
+                (Some(s), StatOp::Count) => Value::Int(s.count as i64),
+                (Some(s), StatOp::Mean) => s.mean.map(Value::Float).unwrap_or(Value::Null),
+                (Some(s), StatOp::Min) => s.min.map(Value::Float).unwrap_or(Value::Null),
+                (Some(s), StatOp::Max) => s.max.map(Value::Float).unwrap_or(Value::Null),
+                (Some(s), StatOp::Sum) => s.sum.map(Value::Float).unwrap_or(Value::Null),
+                (Some(s), StatOp::StdDev) => s.stddev.map(Value::Float).unwrap_or(Value::Null),
+                _ => Value::Null,
+            };
+            attrs.push(val);
+        }
+        out.features.push(rustgis_data::Feature { id: f.id, geometry: f.geometry.clone(), attrs });
+    }
+
+    let n_features = out.features.len();
+    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"));
+        for detail in failure_details.iter().take(5) {
+            messages.push(detail.clone());
+        }
+        if failure_details.len() > 5 {
+            messages.push(format!("... and {} more feature failure(s)", failure_details.len() - 5));
+        }
+    }
+    Ok(ToolResult { output: Output::Layer(out), messages })
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
+    use rustgis_data::Feature;
     use std::collections::HashMap;
 
     struct Map(HashMap<String, Dataset>);
@@ -1436,6 +1824,24 @@ mod tests {
         }
     }
 
+    struct RasterMap<'a> {
+        datasets: HashMap<String, Dataset>,
+        rasters: HashMap<String, &'a rustgis_data::raster::Raster>,
+    }
+
+    impl<'a> Inputs for RasterMap<'a> {
+        fn dataset(&self, layer: &str) -> Option<&Dataset> {
+            self.datasets.get(layer)
+        }
+        fn selection(&self, _: &str) -> &[usize] {
+            &[]
+        }
+        fn raster(&self, layer: &str) -> Option<crate::RasterInput<'_>> {
+            let r = self.rasters.get(layer)?;
+            Some(crate::RasterInput { raster: r, path: None, embedded: None })
+        }
+    }
+
     fn sq(x: f64, y: f64, s: f64) -> Geometry {
         Geometry::Polygon(vec![vec![vec![
             Coord::new(x, y),
@@ -2215,4 +2621,532 @@ mod tests {
         let err = crate::run("spatialJoin", &m_crs, &serde_json::json!({"target": "pts", "join": "polys"})).unwrap_err();
         assert!(err.contains("coordinate reference system"));
     }
+
+    #[test]
+    fn zonal_stats_tool_cases() {
+        let (w, h) = (10u32, 10u32);
+        let mut b1_data: Vec<f32> = Vec::with_capacity((w * h) as usize);
+        let mut b2_data: Vec<f32> = Vec::with_capacity((w * h) as usize);
+        for row in 0..h {
+            for col in 0..w {
+                b1_data.push((col * 10 + row) as f32);
+                let v2 = 100.0 - (col * 10 + row) as f32;
+                if col == 1 && row == 1 {
+                    b2_data.push(-9999.0);
+                } else {
+                    b2_data.push(v2);
+                }
+            }
+        }
+        let gt = rustgis_data::raster::GeoTransform { origin_x: 0.0, origin_y: 10.0, pixel_width: 1.0, pixel_height: -1.0 };
+        let raster = rustgis_data::raster::Raster::new(
+            "elev",
+            w,
+            h,
+            vec![b1_data, b2_data],
+            gt,
+            rustgis_data::raster::Crs::LonLat,
+            Some(-9999.0),
+            rustgis_data::raster::ColorHint::Gray,
+            "Float32",
+        )
+        .unwrap();
+
+        let mut poly_ds = Dataset::new(
+            "zones",
+            GeometryType::Polygon,
+            vec![Field::new("zone_id", FieldType::Integer), Field::new("elev_mean", FieldType::Double)],
+        );
+
+        let sq_geom = |x0: f64, y0: f64, x1: f64, y1: f64| {
+            Geometry::Polygon(vec![vec![vec![
+                Coord::new(x0, y0),
+                Coord::new(x1, y0),
+                Coord::new(x1, y1),
+                Coord::new(x0, y1),
+                Coord::new(x0, y0),
+            ]]])
+        };
+
+        poly_ds.push(Some(sq_geom(1.0, 7.0, 3.0, 9.0)), vec![Value::Int(101), Value::Float(999.0)]);
+        poly_ds.push(Some(sq_geom(20.0, 20.0, 25.0, 25.0)), vec![Value::Int(102), Value::Float(888.0)]);
+
+        let hole_geom = Geometry::Polygon(vec![vec![
+            vec![Coord::new(0.0, 7.0), Coord::new(3.0, 7.0), Coord::new(3.0, 10.0), Coord::new(0.0, 10.0), Coord::new(0.0, 7.0)],
+            vec![Coord::new(1.0, 8.0), Coord::new(2.0, 8.0), Coord::new(2.0, 9.0), Coord::new(1.0, 9.0), Coord::new(1.0, 8.0)],
+        ]]);
+        poly_ds.push(Some(hole_geom), vec![Value::Int(103), Value::Float(777.0)]);
+        poly_ds.push(None, vec![Value::Int(104), Value::Float(666.0)]);
+
+        let input_clone = poly_ds.clone();
+        let m = RasterMap {
+            datasets: HashMap::from([("zones".to_string(), poly_ds)]),
+            rasters: HashMap::from([("elev".to_string(), &raster)]),
+        };
+
+        let res = crate::run(
+            "zonalStats",
+            &m,
+            &serde_json::json!({
+                "raster": "elev",
+                "polygons": "zones",
+                "bands": [1],
+            }),
+        )
+        .unwrap();
+
+        let out = layer(res);
+        assert_eq!(out.features.len(), 4);
+
+        assert_eq!(out.fields[0].name, "zone_id");
+        assert_eq!(out.fields[1].name, "elev_mean");
+        assert_eq!(out.value(0, 0), &Value::Int(101));
+        assert_eq!(out.value(0, 1), &Value::Float(999.0));
+
+        let mean_idx = out.field_index("elev_mean_2").expect("collision should produce elev_mean_2");
+        let min_idx = out.field_index("elev_min").unwrap();
+        let max_idx = out.field_index("elev_max").unwrap();
+        let sum_idx = out.field_index("elev_sum").unwrap();
+        let count_idx = out.field_index("elev_count").unwrap();
+        let stddev_idx = out.field_index("elev_stddev").unwrap();
+
+        assert_eq!(out.value(0, count_idx), &Value::Int(4));
+        assert_eq!(out.value(0, sum_idx), &Value::Float(66.0));
+        assert_eq!(out.value(0, min_idx), &Value::Float(11.0));
+        assert_eq!(out.value(0, max_idx), &Value::Float(22.0));
+        assert_eq!(out.value(0, mean_idx), &Value::Float(16.5));
+        let std_val = out.value(0, stddev_idx).as_f64().unwrap();
+        assert!((std_val - 25.25f64.sqrt()).abs() < 1e-10);
+
+        assert_eq!(out.value(1, count_idx), &Value::Int(0));
+        assert_eq!(out.value(1, sum_idx), &Value::Null);
+        assert_eq!(out.value(1, mean_idx), &Value::Null);
+        assert_eq!(out.value(1, min_idx), &Value::Null);
+        assert_eq!(out.value(1, max_idx), &Value::Null);
+        assert_eq!(out.value(1, stddev_idx), &Value::Null);
+
+        assert_eq!(out.value(2, count_idx), &Value::Int(8));
+
+        assert_eq!(out.value(3, count_idx), &Value::Int(0));
+        assert_eq!(out.value(3, mean_idx), &Value::Null);
+
+        let res_mb = crate::run(
+            "zonalStats",
+            &m,
+            &serde_json::json!({
+                "raster": "elev",
+                "polygons": "zones",
+                "bands": [1, 2],
+                "stats": ["mean", "count"],
+            }),
+        )
+        .unwrap();
+        let out_mb = layer(res_mb);
+        let b1_count = out_mb.field_index("elev_b1_count").unwrap();
+        let b2_count = out_mb.field_index("elev_b2_count").unwrap();
+        let b2_mean = out_mb.field_index("elev_b2_mean").unwrap();
+
+        assert_eq!(out_mb.value(0, b1_count), &Value::Int(4));
+        assert_eq!(out_mb.value(0, b2_count), &Value::Int(3));
+        let b2_m = out_mb.value(0, b2_mean).as_f64().unwrap();
+        assert!((b2_m - 245.0 / 3.0).abs() < 1e-10);
+
+        let mut pt_ds = Dataset::new("points", GeometryType::Point, vec![]);
+        pt_ds.push(Some(Geometry::Point(Coord::new(1.0, 8.0))), vec![]);
+        let m_pts =
+            RasterMap { datasets: HashMap::from([("points".to_string(), pt_ds)]), rasters: HashMap::from([("elev".to_string(), &raster)]) };
+        let err_pt = crate::run(
+            "zonalStats",
+            &m_pts,
+            &serde_json::json!({
+                "raster": "elev",
+                "polygons": "points",
+            }),
+        )
+        .unwrap_err();
+        assert!(err_pt.contains("Point data"));
+
+        let err_band = crate::run(
+            "zonalStats",
+            &m,
+            &serde_json::json!({
+                "raster": "elev",
+                "polygons": "zones",
+                "bands": [5],
+            }),
+        )
+        .unwrap_err();
+        assert!(err_band.contains("band index 5 is invalid"));
+
+        let empty_ds = Dataset::new("empty", GeometryType::Polygon, vec![Field::new("id", FieldType::Integer)]);
+        let m_empty = RasterMap {
+            datasets: HashMap::from([("empty".to_string(), empty_ds)]),
+            rasters: HashMap::from([("elev".to_string(), &raster)]),
+        };
+        let res_empty = crate::run(
+            "zonalStats",
+            &m_empty,
+            &serde_json::json!({
+                "raster": "elev",
+                "polygons": "empty",
+            }),
+        )
+        .unwrap();
+        let out_empty = layer(res_empty);
+        assert_eq!(out_empty.features.len(), 0);
+        assert!(out_empty.field_index("elev_mean").is_some());
+
+        // 9. Bands validation: "all", integral float, decimal rejection, zero rejection, negative rejection
+        let res_all_bands = crate::run(
+            "zonalStats",
+            &m,
+            &serde_json::json!({
+                "raster": "elev",
+                "polygons": "zones",
+                "bands": "all",
+            }),
+        )
+        .unwrap();
+        let out_all = layer(res_all_bands);
+        assert!(out_all.field_index("elev_b1_mean").is_some());
+        assert!(out_all.field_index("elev_b2_mean").is_some());
+
+        let res_float_band = crate::run(
+            "zonalStats",
+            &m,
+            &serde_json::json!({
+                "raster": "elev",
+                "polygons": "zones",
+                "bands": [1.0],
+            }),
+        )
+        .unwrap();
+        let out_float = layer(res_float_band);
+        assert!(out_float.field_index("elev_mean").is_some());
+
+        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": [0]})).is_err());
+        assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "bands": [-1]})).is_err());
+
+        // 10. Projected CRS reprojection: polygon in EPSG:3857 (WebMercator)
+        let merc_tr = rustgis_data::raster::Crs::WebMercator.transformer().unwrap();
+        let p0 = merc_tr.from_lonlat(Coord::new(1.0, 7.0)).unwrap();
+        let p1 = merc_tr.from_lonlat(Coord::new(3.0, 9.0)).unwrap();
+        let merc_geom = Geometry::Polygon(vec![vec![vec![p0, Coord::new(p1.x, p0.y), p1, Coord::new(p0.x, p1.y), p0]]]);
+        let mut merc_poly_ds = Dataset::new("merc_zones", GeometryType::Polygon, vec![Field::new("id", FieldType::Integer)]);
+        merc_poly_ds.crs = Some(rustgis_data::raster::Crs::WebMercator);
+        merc_poly_ds.push(Some(merc_geom.clone()), vec![Value::Int(1)]);
+
+        let m_merc = RasterMap {
+            datasets: HashMap::from([("merc_zones".to_string(), merc_poly_ds)]),
+            rasters: HashMap::from([("elev".to_string(), &raster)]),
+        };
+        let res_merc = crate::run(
+            "zonalStats",
+            &m_merc,
+            &serde_json::json!({
+                "raster": "elev",
+                "polygons": "merc_zones",
+                "bands": [1],
+            }),
+        )
+        .unwrap();
+        let out_merc = layer(res_merc);
+        assert_eq!(out_merc.crs, Some(rustgis_data::raster::Crs::WebMercator));
+        assert_eq!(out_merc.features[0].geometry, Some(merc_geom));
+        let merc_count = out_merc.field_index("elev_count").unwrap();
+        let merc_mean = out_merc.field_index("elev_mean").unwrap();
+        assert_eq!(out_merc.value(0, merc_count), &Value::Int(4));
+        assert_eq!(out_merc.value(0, merc_mean), &Value::Float(16.5));
+
+        // 11. Unsupported polygon CRS returns error
+        let mut bad_crs_ds = Dataset::new("bad_crs", GeometryType::Polygon, vec![]);
+        bad_crs_ds.crs = Some(rustgis_data::raster::Crs::Epsg(60000));
+        bad_crs_ds.push(Some(sq_geom(1.0, 7.0, 3.0, 9.0)), vec![]);
+        let m_bad_crs = RasterMap {
+            datasets: HashMap::from([("bad_crs".to_string(), bad_crs_ds)]),
+            rasters: HashMap::from([("elev".to_string(), &raster)]),
+        };
+        assert!(crate::run("zonalStats", &m_bad_crs, &serde_json::json!({"raster": "elev", "polygons": "bad_crs"})).is_err());
+
+        // 12. Exterior ring dropped during reprojection (does not promote hole)
+        let hole_p0 = merc_tr.from_lonlat(Coord::new(1.0, 8.0)).unwrap();
+        let hole_p1 = merc_tr.from_lonlat(Coord::new(2.0, 9.0)).unwrap();
+        let bad_ext_geom = Geometry::Polygon(vec![vec![
+            vec![p0, p1], // Only 2 points (< 3 points, invalid exterior ring)
+            vec![hole_p0, Coord::new(hole_p1.x, hole_p0.y), hole_p1, Coord::new(hole_p0.x, hole_p1.y), hole_p0], // valid hole
+        ]]);
+        let mut drop_ext_ds = Dataset::new("drop_ext", GeometryType::Polygon, vec![Field::new("id", FieldType::Integer)]);
+        drop_ext_ds.crs = Some(rustgis_data::raster::Crs::WebMercator);
+        drop_ext_ds.push(Some(bad_ext_geom), vec![Value::Int(1)]);
+        let m_drop_ext = RasterMap {
+            datasets: HashMap::from([("drop_ext".to_string(), drop_ext_ds)]),
+            rasters: HashMap::from([("elev".to_string(), &raster)]),
+        };
+        let res_drop_ext = crate::run("zonalStats", &m_drop_ext, &serde_json::json!({"raster": "elev", "polygons": "drop_ext"})).unwrap();
+        assert!(res_drop_ext.messages.iter().any(|m| m.contains("failed coordinate transformation") || m.contains("invalid geometry")));
+        let out_drop_ext = layer(res_drop_ext);
+        let drop_count_idx = out_drop_ext.field_index("elev_count").unwrap();
+        let drop_mean_idx = out_drop_ext.field_index("elev_mean").unwrap();
+        assert_eq!(out_drop_ext.value(0, drop_count_idx), &Value::Int(0));
+        assert_eq!(out_drop_ext.value(0, drop_mean_idx), &Value::Null);
+
+        // 13. Silent vertex loss prevention during reprojection:
+        // A polygon where one exterior ring vertex has non-finite coords (cannot be reprojected to lon/lat).
+        let nan_ext_geom = Geometry::Polygon(vec![vec![vec![
+            Coord::new(f64::NAN, p0.y),
+            Coord::new(p1.x, p0.y),
+            p1,
+            Coord::new(p0.x, p1.y),
+            Coord::new(f64::NAN, p0.y),
+        ]]]);
+        let mut nan_ext_ds = Dataset::new("nan_ext", GeometryType::Polygon, vec![Field::new("id", FieldType::Integer)]);
+        nan_ext_ds.crs = Some(rustgis_data::raster::Crs::WebMercator);
+        nan_ext_ds.push(Some(nan_ext_geom), vec![Value::Int(1)]);
+        let m_nan_ext = RasterMap {
+            datasets: HashMap::from([("nan_ext".to_string(), nan_ext_ds)]),
+            rasters: HashMap::from([("elev".to_string(), &raster)]),
+        };
+        let res_nan_ext = crate::run("zonalStats", &m_nan_ext, &serde_json::json!({"raster": "elev", "polygons": "nan_ext"})).unwrap();
+        assert!(res_nan_ext.messages.iter().any(|m| m.contains("failed coordinate transformation") || m.contains("invalid geometry")));
+        let out_nan_ext = layer(res_nan_ext);
+        let nan_count_idx = out_nan_ext.field_index("elev_count").unwrap();
+        let nan_mean_idx = out_nan_ext.field_index("elev_mean").unwrap();
+        assert_eq!(out_nan_ext.value(0, nan_count_idx), &Value::Int(0));
+        assert_eq!(out_nan_ext.value(0, nan_mean_idx), &Value::Null);
+
+        // 14. Valid projected polygon with holes reprojects both and subtracts holes correctly
+        let merc_outer0 = merc_tr.from_lonlat(Coord::new(0.0, 7.0)).unwrap();
+        let merc_outer1 = merc_tr.from_lonlat(Coord::new(3.0, 10.0)).unwrap();
+        let merc_hole0 = merc_tr.from_lonlat(Coord::new(1.0, 8.0)).unwrap();
+        let merc_hole1 = merc_tr.from_lonlat(Coord::new(2.0, 9.0)).unwrap();
+        let valid_hole_geom = Geometry::Polygon(vec![vec![
+            vec![merc_outer0, Coord::new(merc_outer1.x, merc_outer0.y), merc_outer1, Coord::new(merc_outer0.x, merc_outer1.y), merc_outer0],
+            vec![merc_hole0, Coord::new(merc_hole1.x, merc_hole0.y), merc_hole1, Coord::new(merc_hole0.x, merc_hole1.y), merc_hole0],
+        ]]);
+        let mut hole_ds = Dataset::new("merc_hole", GeometryType::Polygon, vec![Field::new("id", FieldType::Integer)]);
+        hole_ds.crs = Some(rustgis_data::raster::Crs::WebMercator);
+        hole_ds.push(Some(valid_hole_geom), vec![Value::Int(1)]);
+        let m_hole = RasterMap {
+            datasets: HashMap::from([("merc_hole".to_string(), hole_ds)]),
+            rasters: HashMap::from([("elev".to_string(), &raster)]),
+        };
+        let res_hole = crate::run("zonalStats", &m_hole, &serde_json::json!({"raster": "elev", "polygons": "merc_hole"})).unwrap();
+        let out_hole = layer(res_hole);
+        let hole_cnt_idx = out_hole.field_index("elev_count").unwrap();
+        // 9 outer cells - 1 hole cell = 8 cells
+        assert_eq!(out_hole.value(0, hole_cnt_idx), &Value::Int(8));
+
+        assert_eq!(m.datasets["zones"], input_clone);
+    }
+
+    #[test]
+    fn zonal_stats_attribute_padding_and_alignment_regression() {
+        use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster};
+
+        let (w, h) = (4u32, 4u32);
+        let data: Vec<f32> = (1..=16).map(|v| v as f32).collect();
+        let gt = GeoTransform { origin_x: 0.0, origin_y: 4.0, pixel_width: 1.0, pixel_height: -1.0 };
+        let raster = Raster::new("elev", w, h, vec![data], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap();
+
+        let fields =
+            vec![Field::new("name", FieldType::Text), Field::new("category", FieldType::Text), Field::new("score", FieldType::Double)];
+        let mut ds = Dataset::new("zones_pad", GeometryType::Polygon, fields);
+        let sq_geom = |x0: f64, y0: f64, x1: f64, y1: f64| {
+            Geometry::Polygon(vec![vec![vec![
+                Coord::new(x0, y0),
+                Coord::new(x1, y0),
+                Coord::new(x1, y1),
+                Coord::new(x0, y1),
+                Coord::new(x0, y0),
+            ]]])
+        };
+        let geom = Some(sq_geom(0.0, 3.0, 2.0, 4.0)); // 2 cells: 1.0, 2.0
+
+        // Feature 1: shorter attribute array (only 1 attribute instead of 3)
+        ds.features.push(Feature { id: 1, geometry: geom.clone(), attrs: vec![Value::Text("FeatureA".into())] });
+        // Feature 2: full attribute array (3 attributes)
+        ds.features.push(Feature {
+            id: 2,
+            geometry: geom.clone(),
+            attrs: vec![Value::Text("FeatureB".into()), Value::Text("Cat1".into()), Value::Float(88.5)],
+        });
+        // Feature 3: empty attribute array (0 attributes)
+        ds.features.push(Feature { id: 3, geometry: geom, attrs: vec![] });
+
+        let m =
+            RasterMap { datasets: HashMap::from([("zones_pad".to_string(), ds)]), rasters: HashMap::from([("elev".to_string(), &raster)]) };
+
+        let res = crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones_pad"})).unwrap();
+        let out = layer(res);
+
+        // Schema must contain 3 original fields + 6 stat fields
+        assert_eq!(out.fields.len(), 9);
+        assert_eq!(out.fields[0].name, "name");
+        assert_eq!(out.fields[1].name, "category");
+        assert_eq!(out.fields[2].name, "score");
+        assert_eq!(out.fields[3].name, "elev_mean");
+        assert_eq!(out.fields[4].name, "elev_min");
+        assert_eq!(out.fields[5].name, "elev_max");
+        assert_eq!(out.fields[6].name, "elev_sum");
+        assert_eq!(out.fields[7].name, "elev_count");
+        assert_eq!(out.fields[8].name, "elev_stddev");
+
+        // Feature 1: attrs must be padded to 3 before stats fields!
+        let f1 = &out.features[0];
+        assert_eq!(f1.attrs.len(), 9);
+        assert_eq!(f1.attrs[0], Value::Text("FeatureA".into()));
+        assert_eq!(f1.attrs[1], Value::Null, "missing category attribute must be padded with Null");
+        assert_eq!(f1.attrs[2], Value::Null, "missing score attribute must be padded with Null");
+        assert_eq!(f1.attrs[3], Value::Float(1.5), "elev_mean must be at index 3, not in column 1!");
+        assert_eq!(f1.attrs[7], Value::Int(2), "elev_count must be at index 7");
+
+        // Feature 2: complete attributes
+        let f2 = &out.features[1];
+        assert_eq!(f2.attrs.len(), 9);
+        assert_eq!(f2.attrs[0], Value::Text("FeatureB".into()));
+        assert_eq!(f2.attrs[1], Value::Text("Cat1".into()));
+        assert_eq!(f2.attrs[2], Value::Float(88.5));
+        assert_eq!(f2.attrs[3], Value::Float(1.5));
+
+        // Feature 3: empty attributes padded completely
+        let f3 = &out.features[2];
+        assert_eq!(f3.attrs.len(), 9);
+        assert_eq!(f3.attrs[0], Value::Null);
+        assert_eq!(f3.attrs[1], Value::Null);
+        assert_eq!(f3.attrs[2], Value::Null);
+        assert_eq!(f3.attrs[3], Value::Float(1.5));
+    }
+
+    #[test]
+    fn field_name_case_insensitive_and_collision_regression() {
+        let fields = vec![Field::new("Population", FieldType::Double), Field::new("Area_KM2", FieldType::Double)];
+        let mut ds = Dataset::new("cities", GeometryType::None, fields);
+        ds.push(None, vec![Value::Float(100.0), Value::Float(10.0)]);
+        ds.push(None, vec![Value::Float(200.0), Value::Float(20.0)]);
+
+        // 1. Case-insensitive lookup succeeds when casing differs
+        assert_eq!(find_field_case_insensitive(&ds, "population").unwrap(), 0);
+        assert_eq!(find_field_case_insensitive(&ds, "POPULATION").unwrap(), 0);
+        assert_eq!(find_field_case_insensitive(&ds, "area_km2").unwrap(), 1);
+        assert_eq!(find_field_case_insensitive(&ds, "AREA_km2").unwrap(), 1);
+
+        // 2. Missing field
+        assert!(find_field_case_insensitive(&ds, "unknown").is_err());
+
+        // 3. Summary statistics with lowercase lookup preserves original casing in output fields
+        let m = Map(HashMap::from([("cities".to_string(), ds.clone())]));
+        let res = crate::run("summaryStatistics", &m, &serde_json::json!({"input": "cities", "field": "population"})).unwrap();
+        let out = layer(res);
+        assert!(out.field_index("COUNT_Population").is_some(), "output field must preserve original casing");
+        assert!(out.field_index("SUM_Population").is_some());
+        assert!(out.field_index("MEAN_Population").is_some());
+
+        // 4. Exact match precedence and case-insensitive collision detection
+        let colliding_fields = vec![Field::new("Population", FieldType::Double), Field::new("population", FieldType::Double)];
+        let ds_collision = Dataset::new("colliding", GeometryType::None, colliding_fields);
+        // Exact matches must resolve directly:
+        assert_eq!(find_field_case_insensitive(&ds_collision, "population").unwrap(), 1);
+        assert_eq!(find_field_case_insensitive(&ds_collision, "Population").unwrap(), 0);
+        // Inexact matches with multiple candidates are ambiguous:
+        let err = find_field_case_insensitive(&ds_collision, "POPULATION").unwrap_err();
+        assert!(err.contains("ambiguous") && err.contains("collision"), "collision must be explicitly reported: {err}");
+        let err2 = find_field_case_insensitive(&ds_collision, "pOpUlAtIoN").unwrap_err();
+        assert!(err2.contains("collision"));
+    }
+
+    #[test]
+    fn zonal_stats_partial_feature_failures_regression() {
+        use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster};
+
+        let (w, h) = (4u32, 4u32);
+        let data: Vec<f32> = (1..=16).map(|v| v as f32).collect();
+        let gt = GeoTransform { origin_x: 0.0, origin_y: 4.0, pixel_width: 1.0, pixel_height: -1.0 };
+        let raster = Raster::new("elev", w, h, vec![data], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap();
+
+        let fields = vec![Field::new("zone_id", FieldType::Integer)];
+        let mut ds = Dataset::new("mixed_zones", GeometryType::Polygon, fields);
+
+        let sq_ring = |x0: f64, y0: f64, x1: f64, y1: f64| {
+            vec![Coord::new(x0, y0), Coord::new(x1, y0), Coord::new(x1, y1), Coord::new(x0, y1), Coord::new(x0, y0)]
+        };
+
+        // Feature 1: Valid polygon with hole (outer square 0..3, inner hole 1..2)
+        let f1_geom = Geometry::Polygon(vec![vec![
+            sq_ring(0.0, 1.0, 3.0, 4.0), // outer 3x3 = 9 cells
+            sq_ring(1.0, 2.0, 2.0, 3.0), // hole 1x1 = 1 cell
+        ]]);
+        ds.features.push(Feature { id: 101, geometry: Some(f1_geom), attrs: vec![Value::Int(1)] });
+
+        // Feature 2: Defective polygon: hole has fewer than 3 vertices!
+        let f2_geom = Geometry::Polygon(vec![vec![
+            sq_ring(0.0, 1.0, 3.0, 4.0),
+            vec![Coord::new(1.0, 2.0), Coord::new(2.0, 3.0)], // only 2 points!
+        ]]);
+        ds.features.push(Feature { id: 102, geometry: Some(f2_geom), attrs: vec![Value::Int(2)] });
+
+        // Feature 3: Valid simple polygon
+        let f3_geom = Geometry::Polygon(vec![vec![sq_ring(0.0, 3.0, 1.0, 4.0)]]); // 1 cell
+        ds.features.push(Feature { id: 103, geometry: Some(f3_geom), attrs: vec![Value::Int(3)] });
+
+        let m = RasterMap {
+            datasets: HashMap::from([("mixed_zones".to_string(), ds)]),
+            rasters: HashMap::from([("elev".to_string(), &raster)]),
+        };
+
+        let res = crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "mixed_zones"})).unwrap();
+        let messages = res.messages.clone();
+        let out = layer(res);
+
+        assert_eq!(out.features.len(), 3);
+        let count_idx = out.field_index("elev_count").unwrap();
+        let mean_idx = out.field_index("elev_mean").unwrap();
+
+        // Feature 1 succeeded (9 - 1 = 8 cells)
+        assert_eq!(out.value(0, count_idx), &Value::Int(8));
+        assert!(out.value(0, mean_idx).as_f64().is_some());
+
+        // Feature 2 failed (hole was defective -> polygon rejected -> default stats)
+        assert_eq!(out.value(1, count_idx), &Value::Int(0));
+        assert_eq!(out.value(1, mean_idx), &Value::Null);
+
+        // Feature 3 succeeded (1 cell)
+        assert_eq!(out.value(2, count_idx), &Value::Int(1));
+        assert!(out.value(2, mean_idx).as_f64().is_some());
+
+        // Messages must explicitly report feature 102 failure
+        assert!(messages.iter().any(|m| m.contains("1 feature(s) had invalid geometry")));
+        assert!(messages.iter().any(|m| m.contains("feature 102") && m.contains("fewer than 3 vertices")));
+    }
+
+    #[test]
+    fn parse_bands_and_stats_validation_regression() {
+        use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster};
+
+        let (w, h) = (2u32, 2u32);
+        let data: Vec<f32> = vec![1.0, 2.0, 3.0, 4.0];
+        let gt = GeoTransform { origin_x: 0.0, origin_y: 2.0, pixel_width: 1.0, pixel_height: -1.0 };
+        let raster = Raster::new("elev", w, h, vec![data], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap();
+
+        let mut ds = Dataset::new("zones", GeometryType::Polygon, vec![Field::new("id", FieldType::Integer)]);
+        let sq_ring = vec![Coord::new(0.0, 0.0), Coord::new(2.0, 0.0), Coord::new(2.0, 2.0), Coord::new(0.0, 2.0), Coord::new(0.0, 0.0)];
+        ds.push(Some(Geometry::Polygon(vec![vec![sq_ring]])), vec![Value::Int(1)]);
+
+        let m = RasterMap { datasets: HashMap::from([("zones".to_string(), ds)]), rasters: HashMap::from([("elev".to_string(), &raster)]) };
+
+        // 1. Invalid stats inputs
+        assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "stats": []})).is_err());
+        assert!(crate::run("zonalStats", &m, &serde_json::json!({"raster": "elev", "polygons": "zones", "stats": ""})).is_err());
+        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
+        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": 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());
+        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());
+    }
 }
diff --git a/crates/engine/src/cmd/gp.rs b/crates/engine/src/cmd/gp.rs
index 4276822c3..f7c0e8095 100644
--- a/crates/engine/src/cmd/gp.rs
+++ b/crates/engine/src/cmd/gp.rs
@@ -15,6 +15,11 @@ pub fn commands() -> Vec<C> {
             .params(
                 r#"{"target": layer, "join": layer, "match"?: "Intersects" (default) | "Contains" | "Within", "operation"?: "Join one to one" (default) | "Join one to many", "keepAll"?: "Keep all" (default) | "Keep only matching", "prefix"?: "join_", "output"?: text}"#,
             ),
+        C::new("gp.zonalStats", "Zonal Statistics", "Analyze › Raster", zonal_stats)
+            .when(need_layers)
+            .params(
+                r#"{"raster": raster, "polygons": polygon layer, "bands"?: [1-based band numbers] or single band (default 1), "stats"?: ["mean", "min", "max", "sum", "count", "stddev"] or "all" (default "all"), "prefix"?: text, "output"?: text}"#,
+            ),
         C::new("gp.bufferPoint", "Buffer A Point", "Analyze › Geoprocessing", buffer_point).params(
             r#"{"lon": number, "lat": number, "distance": number, "units"?: "Kilometers" (default) | "Meters" | "Miles" | "Feet", "name"?: output}"#,
         ),
@@ -701,6 +706,11 @@ fn spatial_join(s: &mut Session, p: &Value) -> CmdResult {
     run_tool(s, "spatialJoin", p)
 }
 
+/// Calculates raster summary statistics for polygon zones.
+fn zonal_stats(s: &mut Session, p: &Value) -> CmdResult {
+    run_tool(s, "zonalStats", p)
+}
+
 /// Runs a (built-in) tool and applies its result, as `gp.run` does.
 pub fn run_tool(s: &mut Session, tool: &str, params: &Value) -> CmdResult {
     let tool = tool_id(tool).unwrap_or(tool).to_string();
diff --git a/crates/engine/src/cmd/imagery.rs b/crates/engine/src/cmd/imagery.rs
index fecf2d1e2..4fb18210b 100644
--- a/crates/engine/src/cmd/imagery.rs
+++ b/crates/engine/src/cmd/imagery.rs
@@ -70,7 +70,7 @@ pub fn commands() -> Vec<C> {
             r#"{"layer": raster, "index": "ndvi"|"ndwi"|"savi"|"nbr"|"ndbi"|"evi", "red"?, "green"?, "blue"?, "nir"?, "swir"?: band numbers, "name"?: output}"#,
         ),
         C::new("raster.process", "Process", "Analyze › Imagery Tools", process)
-            .params(r#"{"layer": raster, "operation": "clip"|"resample", "extent"?: [minLon, minLat, maxLon, maxLat] (clip; the view when omitted), "polygon"?: polygon layer (clip to its selected or all polygons; cells outside become no-data), "factor"?: n (resample: n×n pixels become one), "name"?: output}"#),
+            .params(r#"{"layer": raster, "operation": "clip"|"resample", "extent"?: [minLon, minLat, maxLon, maxLat] (clip; the view when omitted), "polygon"?: polygon layer (clip to its selected or all polygons; cells outside become no-data), "selectedOnly"?: bool (clip to selected polygons when true, all when false), "factor"?: n (resample: n×n pixels become one), "name"?: output}"#),
         C::new("raster.mosaic", "Mosaic Rasters", "Analyze › Imagery Tools", raster_mosaic).params(
             r#"{"inputs": [raster layer or file path, …] (two or more, in one coordinate system and cell size; where they overlap the first wins), "name"?: output, "path"?: GeoTIFF to write the mosaic to}"#,
         ),
@@ -519,6 +519,28 @@ fn index(s: &mut Session, p: &Value) -> CmdResult {
     add(s, rasterops::derive(&r, &name(p, format!("{}_{idx}", r.name)), vec![out])?, false)
 }
 
+/// Processes a raster via clipping or resampling.
+///
+/// ### Polygon Clipping Semantics (`operation: "clip"`, `"polygon": "<layer>"`)
+/// - **Selection Propagation (`selectedOnly`)**:
+///   - When `Some(true)`: Clips strictly using features currently selected in the polygon layer.
+///     If the selection is empty, returns an actionable error (`"<layer>: no features are selected"`).
+///   - When `Some(false)`: Clips using all features in the polygon layer, ignoring any selection.
+///   - When omitted (`None` or `null`): For backward compatibility, clips to the selection if non-empty,
+///     otherwise clips using all features.
+///   - Non-boolean values are rejected with an explicit validation error.
+/// - **Per-Zone and Overlap Semantics**:
+///   - Every polygon geometry (including multipart geometries) is transformed to WGS 84 lon/lat.
+///   - Polygon exterior rings and interior holes are strictly validated; invalid or failed holes reject
+///     the containing polygon rather than distorting the boundary.
+///   - Polygons are combined into a composite scanline mask over the raster. Holes cut out cells using the
+///     even-odd rule within each polygon. Where multiple polygons overlap, their coverage is combined
+///     (union: any cell falling inside at least one polygon is retained).
+///   - The raster is cropped to the tight bounding box of all selected polygons, and cells outside the
+///     polygons receive NoData.
+/// - **Partial Failure Reporting**:
+///   - Features with invalid geometry or unprojectable coordinates are skipped and logged with a warning.
+///   - If all features fail, the operation fails and reports the specific failure cause of the first feature.
 fn process(s: &mut Session, p: &Value) -> CmdResult {
     let r = raster(s, p)?;
     match ps(p, "operation").ok_or("operation is required")? {
@@ -531,16 +553,59 @@ fn process(s: &mut Session, p: &Value) -> CmdResult {
                 .layer_by_name(pl)
                 .filter(|l| l.geometry_type() == GeometryType::Polygon)
                 .ok_or_else(|| format!("no polygon layer '{pl}'"))?;
-            // The selected polygons, or all of them.
-            let rows: Vec<usize> = if l.selection.is_empty() { (0..l.dataset.features.len()).collect() } else { l.selection.clone() };
-            let polygons: Vec<Vec<Vec<Coord>>> = rows
-                .iter()
-                .filter_map(|row| match l.dataset.features.get(*row).and_then(|f| f.geometry.as_ref()) {
-                    Some(Geometry::Polygon(polys)) => Some(polys.iter().cloned()),
-                    _ => None,
-                })
-                .flatten()
-                .collect();
+            let rows: Vec<usize> = match p.get("selectedOnly") {
+                Some(Value::Bool(true)) => {
+                    if l.selection.is_empty() {
+                        return Err(format!("{pl}: no features are selected"));
+                    }
+                    l.selection.clone()
+                }
+                Some(Value::Bool(false)) => (0..l.dataset.features.len()).collect(),
+                None | Some(Value::Null) => {
+                    // Backward-compatible default: if features are selected, clip to selection; otherwise clip all.
+                    if l.selection.is_empty() { (0..l.dataset.features.len()).collect() } else { l.selection.clone() }
+                }
+                Some(_) => {
+                    return Err("parameter 'selectedOnly' must be a boolean (true or false)".into());
+                }
+            };
+            if rows.is_empty() {
+                return Err(format!("{pl}: has no features to clip by"));
+            }
+            let poly_tr = match l.dataset.crs {
+                Some(crs) if crs != Crs::LonLat && crs != Crs::Epsg(4326) => {
+                    Some(crs.transformer().map_err(|e| format!("{pl} coordinate system isn't supported: {e}"))?)
+                }
+                _ => None,
+            };
+            let mut polygons: Vec<Vec<Vec<Coord>>> = Vec::new();
+            let mut failed_features = Vec::new();
+            for &row in &rows {
+                let Some(f) = l.dataset.features.get(row) else { continue };
+                match &f.geometry {
+                    Some(Geometry::Polygon(polys)) => match rasterops::transform_feature_polygons_to_lonlat(polys, poly_tr.as_ref()) {
+                        Ok(transformed) => {
+                            polygons.extend(transformed);
+                        }
+                        Err(e) => {
+                            failed_features.push((f.id, e));
+                        }
+                    },
+                    _ => {
+                        failed_features
+                            .push((f.id, rasterops::PolygonTransformError::InvalidExterior("missing or non-polygon geometry".into())));
+                    }
+                }
+            }
+            if polygons.is_empty() {
+                if let Some((first_id, first_err)) = failed_features.first() {
+                    return Err(format!("{pl}: has no valid polygon geometries (feature {first_id}: {first_err})"));
+                }
+                return Err(format!("{pl}: has no valid polygon geometries"));
+            }
+            if !failed_features.is_empty() {
+                log::warn!("{pl}: {} feature(s) failed polygon transformation during clipping", failed_features.len());
+            }
             let (cols, rows, mask) =
                 rasterops::polygon_mask(&r, &polygons).ok_or_else(|| format!("no cell of {} falls inside {pl}", r.name))?;
             let out = clip_window(&r, cols, rows, Some(&mask), name(p, format!("{}_clip", r.name)))?;
diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs
index 12795133f..223bf8319 100644
--- a/crates/engine/src/tests.rs
+++ b/crates/engine/src/tests.rs
@@ -11306,7 +11306,14 @@ fn literal_keys(body: &str, prefix: &str) -> Vec<String> {
 #[test]
 fn params_name_every_key_a_command_reads() {
     // Read by a shared helper but of no use to the command.
-    let unused = [("gp.spatialJoin", "input"), ("profile.alongSelected", "row"), ("sources.remove", "name"), ("sources.remove", "extra")];
+    let unused = [
+        ("gp.spatialJoin", "input"),
+        ("gp.zonalStats", "input"),
+        ("gp.zonalStats", "target"),
+        ("profile.alongSelected", "row"),
+        ("sources.remove", "name"),
+        ("sources.remove", "extra"),
+    ];
     let s = Session::new(Project::default());
     let mut missing = Vec::new();
     let found = keys_read_by_commands();
@@ -14249,6 +14256,49 @@ fn dem_tiles_mosaic_clip_by_polygon_share_a_stretch_and_report_overviews() {
     assert!(e.to_string().contains("no cell"), "{e}");
     assert!(s.run("raster.process", &json!({"layer": "dem", "operation": "clip", "polygon": "west"})).is_err());
 
+    // selectedOnly: true with no selected features fails with a descriptive error.
+    s.map_mut().layers.iter_mut().find(|l| l.name == "City").unwrap().selection.clear();
+    let e_sel_none =
+        s.run("raster.process", &json!({"layer": "dem", "operation": "clip", "polygon": "City", "selectedOnly": true})).unwrap_err();
+    assert!(e_sel_none.to_string().contains("no features are selected"), "{e_sel_none}");
+
+    // Add a second feature to City further east in dem's extent.
+    let ring2 = [ll(361_500.0, 5_700_000.0), ll(362_000.0, 5_700_000.0), ll(362_000.0, 5_699_000.0), ll(361_500.0, 5_700_000.0)];
+    s.run("edit.addFeature", &json!({"layer": "City", "geometry": {"type": "Polygon", "coordinates": [ring2]}})).unwrap();
+    assert_eq!(s.map().layer_by_name("City").unwrap().dataset.features.len(), 2);
+
+    // With feature 0 selected: selectedOnly: true clips only feature 0.
+    s.map_mut().layers.iter_mut().find(|l| l.name == "City").unwrap().selection = vec![0];
+    let v_sel = s.run("raster.process", &json!({"layer": "dem", "operation": "clip", "polygon": "City", "selectedOnly": true})).unwrap();
+    let c_sel = raster(&s, &v_sel);
+    assert_eq!(c_sel.bands[0].stats.count, 25);
+
+    // selectedOnly: false clips all features even with selection active.
+    let v_all = s.run("raster.process", &json!({"layer": "dem", "operation": "clip", "polygon": "City", "selectedOnly": false})).unwrap();
+    let c_all = raster(&s, &v_all);
+    assert!(c_all.bands[0].stats.count > 25);
+
+    // An empty polygon layer fails validation.
+    s.run("layer.template", &json!({"template": "polygons", "name": "EmptyPoly"})).unwrap();
+    let e_empty = s.run("raster.process", &json!({"layer": "dem", "operation": "clip", "polygon": "EmptyPoly"})).unwrap_err();
+    assert!(e_empty.to_string().contains("has no features"), "{e_empty}");
+
+    // A polygon layer in a projected CRS (UTM) reprojects and clips accurately.
+    let mut proj_ds = rustgis_data::Dataset::new("CityProj", rustgis_data::GeometryType::Polygon, vec![]);
+    proj_ds.crs = Some(utm);
+    let utm_ring = [
+        rustgis_geom::Coord::new(360_500.0, 5_700_000.0),
+        rustgis_geom::Coord::new(361_500.0, 5_700_000.0),
+        rustgis_geom::Coord::new(361_500.0, 5_699_500.0),
+        rustgis_geom::Coord::new(360_500.0, 5_700_000.0),
+    ];
+    proj_ds.push(Some(rustgis_data::Geometry::Polygon(vec![vec![utm_ring.to_vec()]])), vec![]);
+    s.map_mut().add_layer(rustgis_map::Layer::new(0, proj_ds));
+    let v_proj = s.run("raster.process", &json!({"layer": "dem", "operation": "clip", "polygon": "CityProj"})).unwrap();
+    let c_proj = raster(&s, &v_proj);
+    assert_eq!((c_proj.width, c_proj.height, c_proj.crs), (9, 5, utm));
+    assert_eq!(c_proj.bands[0].stats.count, 25);
+
     // One stretch shared across tiles, as fixed bounds.
     s.run("raster.stretch", &json!({"layer": "west", "type": "custom", "min": 100, "max": 300})).unwrap();
     s.run("raster.stretch", &json!({"layer": "east", "like": "west"})).unwrap();
@@ -15204,6 +15254,332 @@ fn scans_have_no_coordinates_to_give_away() {
     let _ = std::fs::remove_dir_all(&dir);
 }
 
+#[test]
+fn gp_zonal_stats_end_to_end_and_undo() {
+    use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster};
+    use rustgis_data::{Dataset, Field, FieldType, Geometry, GeometryType, Value};
+    use rustgis_geom::Coord;
+
+    let mut s = world();
+
+    let (w, h) = (10u32, 10u32);
+    let mut data = Vec::with_capacity((w * h) as usize);
+    for row in 0..h {
+        for col in 0..w {
+            data.push((col * 10 + row) as f32);
+        }
+    }
+    let gt = GeoTransform { origin_x: 0.0, origin_y: 10.0, pixel_width: 1.0, pixel_height: -1.0 };
+    let raster = Raster::new("dem", w, h, vec![data], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap();
+    let raster_layer = rustgis_map::RasterLayer::from_path(raster, "dem.tif");
+    s.map_mut().add_raster("dem", raster_layer);
+
+    let mut poly_ds = Dataset::new("districts", GeometryType::Polygon, vec![Field::new("code", FieldType::Text)]);
+    let sq = |x0: f64, y0: f64, x1: f64, y1: f64| {
+        Geometry::Polygon(vec![vec![vec![
+            Coord::new(x0, y0),
+            Coord::new(x1, y0),
+            Coord::new(x1, y1),
+            Coord::new(x0, y1),
+            Coord::new(x0, y0),
+        ]]])
+    };
+    poly_ds.push(Some(sq(1.0, 7.0, 3.0, 9.0)), vec![Value::Text("D1".into())]);
+    poly_ds.push(Some(sq(50.0, 50.0, 52.0, 52.0)), vec![Value::Text("D2".into())]);
+    s.map_mut().add_layer(rustgis_map::Layer::new(0, poly_ds));
+
+    let layers_before = s.map().layers.len();
+
+    let res = s.run(
+        "gp.zonalStats",
+        &json!({
+            "raster": "dem",
+            "polygons": "districts",
+            "output": "zonal_out",
+        }),
+    );
+    assert!(res.is_ok(), "gp.zonalStats failed: {res:?}");
+    assert_eq!(s.map().layers.len(), layers_before + 1);
+
+    let out_layer = s.map().layer_by_name("zonal_out").expect("output layer should exist");
+    assert_eq!(out_layer.dataset.features.len(), 2);
+
+    let mean_idx = out_layer.dataset.field_index("dem_mean").expect("dem_mean field");
+    let count_idx = out_layer.dataset.field_index("dem_count").expect("dem_count field");
+    let min_idx = out_layer.dataset.field_index("dem_min").expect("dem_min field");
+    let max_idx = out_layer.dataset.field_index("dem_max").expect("dem_max field");
+    let sum_idx = out_layer.dataset.field_index("dem_sum").expect("dem_sum field");
+    let stddev_idx = out_layer.dataset.field_index("dem_stddev").expect("dem_stddev field");
+
+    assert_eq!(out_layer.dataset.value(0, 0), &Value::Text("D1".into()));
+    assert_eq!(out_layer.dataset.value(0, count_idx), &Value::Int(4));
+    assert_eq!(out_layer.dataset.value(0, mean_idx), &Value::Float(16.5));
+    assert_eq!(out_layer.dataset.value(0, min_idx), &Value::Float(11.0));
+    assert_eq!(out_layer.dataset.value(0, max_idx), &Value::Float(22.0));
+    assert_eq!(out_layer.dataset.value(0, sum_idx), &Value::Float(66.0));
+    let std = out_layer.dataset.value(0, stddev_idx).as_f64().unwrap();
+    assert!((std - 25.25f64.sqrt()).abs() < 1e-10);
+
+    assert_eq!(out_layer.dataset.value(1, 0), &Value::Text("D2".into()));
+    assert_eq!(out_layer.dataset.value(1, count_idx), &Value::Int(0));
+    assert_eq!(out_layer.dataset.value(1, mean_idx), &Value::Null);
+    assert_eq!(out_layer.dataset.value(1, min_idx), &Value::Null);
+    assert_eq!(out_layer.dataset.value(1, max_idx), &Value::Null);
+    assert_eq!(out_layer.dataset.value(1, sum_idx), &Value::Null);
+    assert_eq!(out_layer.dataset.value(1, stddev_idx), &Value::Null);
+
+    assert!(s.undo());
+    assert_eq!(s.map().layers.len(), layers_before);
+    assert!(s.map().layer_by_name("zonal_out").is_none());
+
+    assert!(s.redo());
+    assert_eq!(s.map().layers.len(), layers_before + 1);
+    assert!(s.map().layer_by_name("zonal_out").is_some());
+
+    let run_res = s.run(
+        "gp.run",
+        &json!({
+            "tool": "zonalStats",
+            "params": {
+                "raster": "dem",
+                "polygons": "districts",
+                "stats": ["mean", "count"],
+                "output": "via_gp_run",
+            }
+        }),
+    );
+    assert!(run_res.is_ok(), "gp.run with zonalStats failed: {run_res:?}");
+    let run_layer = s.map().layer_by_name("via_gp_run").expect("via_gp_run layer should exist");
+    assert!(run_layer.dataset.field_index("dem_mean").is_some());
+    assert!(run_layer.dataset.field_index("dem_count").is_some());
+    assert!(run_layer.dataset.field_index("dem_min").is_none());
+
+    // Projected CRS (WebMercator) reprojection end-to-end
+    let merc_tr = Crs::WebMercator.transformer().unwrap();
+    let mp0 = merc_tr.from_lonlat(Coord::new(1.0, 7.0)).unwrap();
+    let mp1 = merc_tr.from_lonlat(Coord::new(3.0, 9.0)).unwrap();
+    let mut merc_ds = Dataset::new("merc_districts", GeometryType::Polygon, vec![Field::new("code", FieldType::Text)]);
+    merc_ds.crs = Some(Crs::WebMercator);
+    merc_ds.push(Some(sq(mp0.x, mp0.y, mp1.x, mp1.y)), vec![Value::Text("M1".into())]);
+    s.map_mut().add_layer(rustgis_map::Layer::new(0, merc_ds));
+
+    let merc_res = s.run(
+        "gp.zonalStats",
+        &json!({
+            "raster": "dem",
+            "polygons": "merc_districts",
+            "output": "merc_out",
+            "bands": [1.0],
+        }),
+    );
+    assert!(merc_res.is_ok(), "gp.zonalStats with WebMercator failed: {merc_res:?}");
+    let merc_layer = s.map().layer_by_name("merc_out").expect("merc_out should exist");
+    assert_eq!(merc_layer.dataset.crs, Some(Crs::WebMercator));
+    let m_count_idx = merc_layer.dataset.field_index("dem_count").unwrap();
+    let m_mean_idx = merc_layer.dataset.field_index("dem_mean").unwrap();
+    assert_eq!(merc_layer.dataset.value(0, m_count_idx), &Value::Int(4));
+    assert_eq!(merc_layer.dataset.value(0, m_mean_idx), &Value::Float(16.5));
+
+    // "bands": "all"
+    let all_res = s.run(
+        "gp.zonalStats",
+        &json!({
+            "raster": "dem",
+            "polygons": "districts",
+            "output": "all_bands_out",
+            "bands": "all",
+        }),
+    );
+    assert!(all_res.is_ok(), "gp.zonalStats with bands: all failed: {all_res:?}");
+
+    assert!(s.run("gp.zonalStats", &json!({})).is_err());
+    assert!(s.run("gp.zonalStats", &json!({"raster": "nonexistent", "polygons": "districts"})).is_err());
+    assert!(s.run("gp.zonalStats", &json!({"raster": "dem", "polygons": "Cities"})).is_err());
+    assert!(s.run("gp.zonalStats", &json!({"raster": "dem", "polygons": "districts", "bands": [999]})).is_err());
+    assert!(s.run("gp.zonalStats", &json!({"raster": "dem", "polygons": "districts", "bands": [0]})).is_err());
+    assert!(s.run("gp.zonalStats", &json!({"raster": "dem", "polygons": "districts", "bands": [-1]})).is_err());
+    assert!(s.run("gp.zonalStats", &json!({"raster": "dem", "polygons": "districts", "bands": [1.5]})).is_err());
+}
+
+#[test]
+fn raster_clip_polygon_selection_semantics_and_partial_failure() {
+    use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster};
+    use rustgis_data::{Dataset, Field, FieldType, Geometry, GeometryType, Value};
+    use rustgis_geom::Coord;
+
+    let mut s = world();
+
+    // 10x10 raster: origin (0.0, 10.0), cell size (1.0, -1.0)
+    let (w, h) = (10u32, 10u32);
+    let mut data = Vec::with_capacity((w * h) as usize);
+    for row in 0..h {
+        for col in 0..w {
+            data.push((col * 10 + row) as f32);
+        }
+    }
+    let gt = GeoTransform { origin_x: 0.0, origin_y: 10.0, pixel_width: 1.0, pixel_height: -1.0 };
+    let raster = Raster::new("dem", w, h, vec![data], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap();
+    let raster_layer = rustgis_map::RasterLayer::from_path(raster, "dem.tif");
+    s.map_mut().add_raster("dem", raster_layer);
+
+    let sq = |x0: f64, y0: f64, x1: f64, y1: f64| {
+        Geometry::Polygon(vec![vec![vec![
+            Coord::new(x0, y0),
+            Coord::new(x1, y0),
+            Coord::new(x1, y1),
+            Coord::new(x0, y1),
+            Coord::new(x0, y0),
+        ]]])
+    };
+
+    let mut poly_ds = Dataset::new("clip_zones", GeometryType::Polygon, vec![Field::new("name", FieldType::Text)]);
+    // Feature 0: [1.0, 7.0] to [3.0, 9.0]
+    poly_ds.push(Some(sq(1.0, 7.0, 3.0, 9.0)), vec![Value::Text("Z1".into())]);
+    // Feature 1: [2.0, 6.0] to [4.0, 8.0] (overlaps with Z1)
+    poly_ds.push(Some(sq(2.0, 6.0, 4.0, 8.0)), vec![Value::Text("Z2".into())]);
+    // Feature 2: [6.0, 1.0] to [8.0, 3.0] (separate area)
+    poly_ds.push(Some(sq(6.0, 1.0, 8.0, 3.0)), vec![Value::Text("Z3".into())]);
+    // Feature 3: Degenerate hole (< 3 coords) -> should trigger rejection of this feature
+    let bad_hole_poly = Geometry::Polygon(vec![vec![
+        vec![Coord::new(4.0, 4.0), Coord::new(6.0, 4.0), Coord::new(6.0, 6.0), Coord::new(4.0, 6.0), Coord::new(4.0, 4.0)],
+        vec![Coord::new(4.5, 4.5), Coord::new(5.0, 5.0)], // invalid hole with only 2 points
+    ]]);
+    poly_ds.push(Some(bad_hole_poly), vec![Value::Text("Z4_bad".into())]);
+
+    let poly_layer_id = s.map_mut().add_layer(rustgis_map::Layer::new(0, poly_ds));
+
+    // Test 1: selectedOnly: "invalid" or non-boolean returns validation error
+    let err_type = s.run(
+        "raster.process",
+        &json!({
+            "layer": "dem",
+            "operation": "clip",
+            "polygon": "clip_zones",
+            "selectedOnly": "yes",
+        }),
+    );
+    assert!(err_type.is_err());
+    assert!(err_type.unwrap_err().contains("must be a boolean"));
+
+    // Test 2: selectedOnly: true with empty selection returns error
+    let err_empty_sel = s.run(
+        "raster.process",
+        &json!({
+            "layer": "dem",
+            "operation": "clip",
+            "polygon": "clip_zones",
+            "selectedOnly": true,
+        }),
+    );
+    assert!(err_empty_sel.is_err());
+    assert!(err_empty_sel.unwrap_err().contains("no features are selected"));
+
+    // Test 3: selectedOnly: true with Feature 0 selected
+    {
+        let l = s.map_mut().layer_mut(poly_layer_id).unwrap();
+        l.selection = vec![0];
+    }
+    let res_sel = s.run(
+        "raster.process",
+        &json!({
+            "layer": "dem",
+            "operation": "clip",
+            "polygon": "clip_zones",
+            "selectedOnly": true,
+            "name": "dem_sel_0",
+        }),
+    );
+    assert!(res_sel.is_ok(), "clip selected failed: {res_sel:?}");
+    let sel_layer = s.map().layer_by_name("dem_sel_0").expect("dem_sel_0 exists");
+    let sel_raster = sel_layer.raster.as_ref().unwrap().data.as_ref().unwrap();
+    // Feature 0 covers x: [1.0, 3.0], y: [7.0, 9.0] -> 2 cols x 2 rows
+    assert_eq!(sel_raster.width, 2);
+    assert_eq!(sel_raster.height, 2);
+
+    // Test 4: Overlapping polygons combined coverage (Features 0 and 1)
+    {
+        let l = s.map_mut().layer_mut(poly_layer_id).unwrap();
+        l.selection = vec![0, 1];
+    }
+    let res_overlap = s.run(
+        "raster.process",
+        &json!({
+            "layer": "dem",
+            "operation": "clip",
+            "polygon": "clip_zones",
+            "selectedOnly": true,
+            "name": "dem_overlap",
+        }),
+    );
+    assert!(res_overlap.is_ok(), "clip overlap failed: {res_overlap:?}");
+    let ov_layer = s.map().layer_by_name("dem_overlap").expect("dem_overlap exists");
+    let ov_raster = ov_layer.raster.as_ref().unwrap().data.as_ref().unwrap();
+    // Features 0 & 1 union bounding box: x: [1.0, 4.0], y: [6.0, 9.0] -> 3 cols x 3 rows
+    assert_eq!(ov_raster.width, 3);
+    assert_eq!(ov_raster.height, 3);
+
+    // Test 5: selectedOnly: false ignores active selection and clips using all valid features
+    let res_all = s.run(
+        "raster.process",
+        &json!({
+            "layer": "dem",
+            "operation": "clip",
+            "polygon": "clip_zones",
+            "selectedOnly": false,
+            "name": "dem_all",
+        }),
+    );
+    assert!(res_all.is_ok(), "clip all failed: {res_all:?}");
+    let all_layer = s.map().layer_by_name("dem_all").expect("dem_all exists");
+    let all_raster = all_layer.raster.as_ref().unwrap().data.as_ref().unwrap();
+    // Features 0, 1, 2 bounding box spans x: [1.0, 8.0], y: [1.0, 9.0] -> 7 cols x 8 rows
+    assert_eq!(all_raster.width, 7);
+    assert_eq!(all_raster.height, 8);
+
+    // Test 6: Partial failure reporting: Feature 0 (valid) + Feature 3 (invalid hole)
+    // Feature 3 should be skipped with a warning, Feature 0 should succeed
+    {
+        let l = s.map_mut().layer_mut(poly_layer_id).unwrap();
+        l.selection = vec![0, 3];
+    }
+    let res_partial = s.run(
+        "raster.process",
+        &json!({
+            "layer": "dem",
+            "operation": "clip",
+            "polygon": "clip_zones",
+            "selectedOnly": true,
+            "name": "dem_partial",
+        }),
+    );
+    assert!(res_partial.is_ok(), "clip partial should succeed for valid feature: {res_partial:?}");
+    let partial_layer = s.map().layer_by_name("dem_partial").expect("dem_partial exists");
+    let partial_raster = partial_layer.raster.as_ref().unwrap().data.as_ref().unwrap();
+    assert_eq!(partial_raster.width, 2);
+    assert_eq!(partial_raster.height, 2);
+
+    // Test 7: Total failure when all selected features fail validation (only Feature 3 selected)
+    {
+        let l = s.map_mut().layer_mut(poly_layer_id).unwrap();
+        l.selection = vec![3];
+    }
+    let res_bad = s.run(
+        "raster.process",
+        &json!({
+            "layer": "dem",
+            "operation": "clip",
+            "polygon": "clip_zones",
+            "selectedOnly": true,
+            "name": "dem_bad",
+        }),
+    );
+    assert!(res_bad.is_err());
+    let err_msg = res_bad.unwrap_err();
+    assert!(err_msg.contains("has no valid polygon geometries"));
+    assert!(err_msg.contains("feature 4"));
+    assert!(err_msg.contains("hole 1 is invalid"));
+}
+
 #[test]
 fn layer_selection_replaces_toggles_and_extends() {
     let mut s = world();
diff --git a/crates/ui-egui/src/analysispane.rs b/crates/ui-egui/src/analysispane.rs
index bcab7947b..22d02bcb0 100644
--- a/crates/ui-egui/src/analysispane.rs
+++ b/crates/ui-egui/src/analysispane.rs
@@ -36,6 +36,16 @@ pub struct Workflow {
     /// The Script pane's last run, shown above its results.
     pub script_last: Option<ScriptRun>,
     pub trend: TrendForm,
+    pub selected_only: Option<bool>,
+    pub zonal_stats: ZonalStatsForm,
+}
+
+/// The Zonal Statistics section's inputs.
+#[derive(Default)]
+pub struct ZonalStatsForm {
+    pub bands: String,
+    pub prefix: String,
+    pub output: String,
 }
 
 /// A Script pane run, as the line above its results tells it: its number this session (so a run
@@ -124,6 +134,7 @@ pub const SECTIONS: &[(&str, &str)] = &[
     ("rasterFunctions", "Raster Functions"),
     ("functionEditor", "Function Editor"),
     ("suitability", "Suitability Modeler"),
+    ("zonalStats", "Zonal Statistics"),
     ("visibility", "Visibility Analysis"),
     ("lineOfSight", "Exploratory 3D Analysis"),
     ("interpolate", "Geostatistical Wizard"),
@@ -206,6 +217,7 @@ pub fn workflows(app: &mut RustGisApp, ui: &mut Ui, t: &Tokens) {
     let inputs = match w.section {
         "rasterFunctions" | "functionEditor" | "suitability" | "visibility" | "lineOfSight" | "georeference" | "change" | "spectral"
         | "classify" | "supervised" | "deepLearning" | "hosted" | "indices" | "process" | "pixels" | "mensuration" | "reality"
+        | "zonalStats"
         // The trend's one-layer mode takes a raster (a Zarr variable, or a raster of several bands).
         | "trend" => &rasters,
         "interpolate" => &points,
@@ -216,6 +228,9 @@ pub fn workflows(app: &mut RustGisApp, ui: &mut Ui, t: &Tokens) {
     if !inputs.contains(&w.layer) {
         w.layer = inputs.first().cloned().unwrap_or_default();
     }
+    if (w.section == "zonalStats" || w.section == "process") && !polygons.contains(&w.layer2) {
+        w.layer2 = polygons.first().cloned().unwrap_or_default();
+    }
     ui.add_space(6.0);
     match app.workflow.section {
         "rasterFunctions" | "functionEditor" => {
@@ -301,6 +316,66 @@ pub fn workflows(app: &mut RustGisApp, ui: &mut Ui, t: &Tokens) {
                 run(app, "analysis.suitability", json!({"criteria": [c1, c2]}));
             }
         }
+        "zonalStats" => {
+            heading(ui, "Calculate raster statistics per polygon");
+            ui.label(egui::RichText::new(tl!("Raster layer")).small().weak());
+            combo(ui, "wf_zs_raster", "Raster layer", &mut app.workflow.layer, &rasters);
+            ui.label(egui::RichText::new(tl!("Polygon layer")).small().weak());
+            combo(ui, "wf_zs_poly", "Polygon layer", &mut app.workflow.layer2, &polygons);
+            ui.horizontal(|ui| {
+                let band_label = ui.label(tl!("Band"));
+                ui.add(egui::TextEdit::singleline(&mut app.workflow.zonal_stats.bands).hint_text("1").desired_width(40.0))
+                    .labelled_by(band_label.id);
+                let prefix_label = ui.label(tl!("Prefix"));
+                ui.add(egui::TextEdit::singleline(&mut app.workflow.zonal_stats.prefix).hint_text("elev").desired_width(70.0))
+                    .labelled_by(prefix_label.id);
+            });
+            ui.horizontal(|ui| {
+                let output_label = ui.label(tl!("Output"));
+                ui.add(
+                    egui::TextEdit::singleline(&mut app.workflow.zonal_stats.output)
+                        .hint_text(tl!("output layer name"))
+                        .desired_width(120.0),
+                )
+                .labelled_by(output_label.id);
+            });
+            let band_str = app.workflow.zonal_stats.bands.trim();
+            let total_bands = app
+                .session
+                .map()
+                .layer_by_name(&app.workflow.layer)
+                .and_then(|l| l.raster.as_ref())
+                .and_then(|r| r.data.as_ref())
+                .map(|d| d.bands.len())
+                .unwrap_or(1);
+            let band_res =
+                if band_str.is_empty() { Ok(vec![0]) } else { rustgis_engine::analysis::rasterops::parse_bands_str(band_str, total_bands) };
+            let band_valid = band_res.is_ok();
+            let band_err = band_res.err();
+            let layers_ready = !app.workflow.layer.is_empty() && !app.workflow.layer2.is_empty();
+            let can_run = layers_ready && band_valid;
+            if !layers_ready {
+                ui.label(egui::RichText::new(tl!("Add a raster and a polygon layer to calculate zonal statistics.")).small().weak());
+            } else if let Some(err) = band_err {
+                ui.label(egui::RichText::new(err).small().color(t.amber));
+            }
+            if ui.add_enabled(can_run, egui::Button::new(tl!("Calculate Statistics"))).clicked() {
+                let mut params = json!({
+                    "raster": app.workflow.layer,
+                    "polygons": app.workflow.layer2,
+                });
+                if !band_str.is_empty() {
+                    params["bands"] = json!(band_str);
+                }
+                if !app.workflow.zonal_stats.prefix.trim().is_empty() {
+                    params["prefix"] = json!(app.workflow.zonal_stats.prefix.trim());
+                }
+                if !app.workflow.zonal_stats.output.trim().is_empty() {
+                    params["output"] = json!(app.workflow.zonal_stats.output.trim());
+                }
+                run(app, "gp.zonalStats", params);
+            }
+        }
         "visibility" => {
             heading(ui, "Viewshed from a point");
             combo(ui, "wf_v", "Layer", &mut app.workflow.layer, &rasters);
@@ -611,13 +686,69 @@ pub fn workflows(app: &mut RustGisApp, ui: &mut Ui, t: &Tokens) {
         }
         "process" => {
             heading(ui, "Clip or resample");
-            combo(ui, "wf_pr", "Layer", &mut app.workflow.layer, &rasters);
+            combo(ui, "wf_pr", "Raster layer", &mut app.workflow.layer, &rasters);
+
+            ui.add_space(4.0);
+            ui.label(egui::RichText::new(tl!("Clip to polygon layer")).small().weak());
+            let has_polys = !polygons.is_empty();
+            if has_polys {
+                let layer2_before = app.workflow.layer2.clone();
+                combo(ui, "wf_pr_poly", "Polygon layer", &mut app.workflow.layer2, &polygons);
+                if app.workflow.layer2 != layer2_before {
+                    app.workflow.selected_only = None;
+                }
+                let poly_layer = app.session.map().layer_by_name(&app.workflow.layer2);
+                let sel_count = poly_layer.as_ref().map_or(0, |l| l.selection.len());
+                let total_features = poly_layer.as_ref().map_or(0, |l| l.dataset.features.len());
+
+                let effective_selected_only = app.workflow.selected_only.unwrap_or(sel_count > 0);
+                let mut checked = effective_selected_only;
+                ui.horizontal(|ui| {
+                    if ui.checkbox(&mut checked, tl!("Selected only")).changed() {
+                        app.workflow.selected_only = Some(checked);
+                    }
+                    if sel_count > 0 {
+                        let text = crate::i18n::fmt(
+                            "({selected} of {total})",
+                            &[("selected", &sel_count.to_string()), ("total", &total_features.to_string())],
+                        );
+                        ui.label(egui::RichText::new(text).small().weak());
+                    }
+                });
+
+                let sel_invalid = effective_selected_only && sel_count == 0;
+                let can_clip_poly = !app.workflow.layer.is_empty() && !app.workflow.layer2.is_empty() && !sel_invalid;
+
+                if sel_invalid {
+                    ui.label(egui::RichText::new(tl!("No features are selected in this layer.")).small().color(t.amber));
+                }
+
+                if ui.add_enabled(can_clip_poly, egui::Button::new(tl!("Clip To Polygon"))).clicked() {
+                    let mut p = json!({
+                        "layer": app.workflow.layer,
+                        "operation": "clip",
+                        "polygon": app.workflow.layer2,
+                    });
+                    p["selectedOnly"] = json!(app.workflow.selected_only);
+                    run(app, "raster.process", p);
+                }
+            } else {
+                ui.label(egui::RichText::new(tl!("No polygon layer available in the map.")).small().weak());
+            }
+
+            ui.add_space(4.0);
             ui.horizontal(|ui| {
-                if ui.button(tl!("Clip To Map View")).clicked() {
+                if ui.add_enabled(!app.workflow.layer.is_empty(), egui::Button::new(tl!("Clip To Map View"))).clicked() {
                     run(app, "raster.process", json!({"layer": app.workflow.layer, "operation": "clip"}));
                 }
-                ui.add(egui::TextEdit::singleline(&mut app.workflow.number).hint_text(tl!("factor")).desired_width(50.0));
-                if ui.button(tl!("Resample")).clicked() {
+            });
+
+            ui.add_space(4.0);
+            ui.horizontal(|ui| {
+                let factor_label = ui.label(tl!("factor"));
+                ui.add(egui::TextEdit::singleline(&mut app.workflow.number).hint_text(tl!("factor")).desired_width(50.0))
+                    .labelled_by(factor_label.id);
+                if ui.add_enabled(!app.workflow.layer.is_empty(), egui::Button::new(tl!("Resample"))).clicked() {
                     run(
                         app,
                         "raster.process",
@@ -1245,6 +1376,11 @@ mod tests {
                     assert!(h.query_by_role_and_label(egui::accesskit::Role::ComboBox, label).is_some());
                 }
             }
+            if section == "zonalStats" {
+                for label in ["Band", "Prefix", "Output"] {
+                    assert!(h.query_by_role_and_label(egui::accesskit::Role::TextInput, label).is_some());
+                }
+            }
         }
     }
 }
diff --git a/crates/ui-egui/tests/smoke.rs b/crates/ui-egui/tests/smoke.rs
index 022f12757..f520ea392 100644
--- a/crates/ui-egui/tests/smoke.rs
+++ b/crates/ui-egui/tests/smoke.rs
@@ -235,6 +235,30 @@ fn command_search_workflow_collects_inputs_before_georeferencing() {
     assert_eq!(h.state().ui.right, Some(rustgis_ui_egui::RightPane::Workflows));
 }
 
+/// Clipping follows a selection until the user explicitly opts out, even after changing it.
+#[test]
+fn polygon_clip_defaults_to_selection_and_honors_opt_out() {
+    let mut h = harness();
+    h.state_mut().workflow.layer2 = "Countries".into();
+    h.state_mut().run("ui.pane.workflows", json!(null)).unwrap();
+    h.state_mut().workflow.section = "process";
+    h.state_mut().session.map_mut().layers.iter_mut().find(|l| l.name == "Countries").unwrap().selection.push(0);
+    h.run_steps(4);
+    assert_eq!(h.state().workflow.selected_only, None);
+    let total = h.state().session.map().layer_by_name("Countries").unwrap().dataset.features.len();
+    assert!(h.query_by_label(&format!("(1 of {total})")).is_some());
+    h.get_by_label("Selected only").click();
+    h.run_steps(3);
+    assert_eq!(h.state().workflow.selected_only, Some(false));
+    h.state_mut().session.map_mut().layers.iter_mut().find(|l| l.name == "Countries").unwrap().selection.push(1);
+    h.run_steps(3);
+    assert_eq!(h.state().workflow.selected_only, Some(false));
+    assert!(h.query_by_label(&format!("(2 of {total})")).is_some());
+    h.get_by_label("Selected only").click();
+    h.run_steps(3);
+    assert_eq!(h.state().workflow.selected_only, Some(true));
+}
+
 #[cfg(not(target_arch = "wasm32"))]
 #[test]
 fn mcp_relay_early_start_wakes_queued_request_and_updates_ui() {
diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md
index 04e0c314c..9c7f52bac 100644
--- a/docs/user-guide/analysis.md
+++ b/docs/user-guide/analysis.md
@@ -533,6 +533,20 @@ rustgis-cli run --data cobertura_1985.tif --data cobertura_1986.tif --data cober
 rustgis-cli run --cmd 'cube.reduce={"url":"https://data.inpe.br/bdc/stac/v1","collection":"mod13q1-6.1","asset":"NDVI","bbox":[-44.8,-7.5,-44.3,-7.0],"start":"2001","end":"2025","cellSize":1000,"stat":"mean","by":"year","minValid":0.5,"name":"NDVI"}' --save urucui.rgis
 ```
 
+## Zonal statistics and polygon clipping
+
+Zonal Statistics uses the selected polygon features when the layer has a selection. It keeps
+original geometries and fields, and reports features with missing or rejected polygon parts.
+If any polygon part fails validation or transformation, the feature receives count 0 and null
+statistics, so results never summarize only the surviving parts. Large
+overlapping zones each require a separate pass over their raster window.
+
+For `raster.process` polygon clipping, `selectedOnly: true` requires a selection,
+`selectedOnly: false` uses every polygon, and omitting it uses the selection when present.
+The Process pane initially checks Selected only when the chosen layer has a selection, and sends
+the choice explicitly. Clipping and Zonal Statistics share
+polygon transformation rules so a failed hole cannot silently enlarge the sampled area.
+
 ## Deep learning and hosted imagery
 
 **Analyze › Deep Learning Tools** trains a neural network (a multilayer perceptron) on the pixels