fix(analysis): use normalized polygon coordinates for zonal statistics

crates/analysis/src/tools.rs

@@ -1037,11 +1037,6 @@ fn idw_auto_cell(extent: &Extent, n: usize) -> f64 {
 /// Ordinary Kriging interpolation: builds prediction and error surfaces from spatial autocorrelation.
 fn kriging_interpolation(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
     let ds = input(inp, p, "input")?;
-    if let Some(crs) = ds.crs
-        && crs != rustgis_data::raster::Crs::LonLat
-    {
-        return Err("kriging interpolation expects coordinates in degrees (WGS 84 / EPSG:4326)".into());
-    }
     let field = s(p, "field");
     let fi = find_field_case_insensitive(ds, field)?;
     let mut raw_pts: Vec<(Coord, f64)> = Vec::new();
@@ -1906,12 +1901,6 @@ fn zonal_stats(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
             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)?;
@@ -1964,7 +1953,8 @@ fn zonal_stats(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
         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()) {
+                // Dataset geometry is already lon/lat; crs records the source system.
+                match crate::rasterops::transform_feature_polygons_to_lonlat(polys, None) {
                     Ok(polys_lonlat) => {
                         let mut failed_parts = Vec::new();
                         let st = crate::rasterops::zonal_statistics_with_failures(r, &bands, &polys_lonlat, &mut failed_parts)?;
@@ -2012,7 +2002,7 @@ fn zonal_stats(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
     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"));
+        messages.push(format!("{skipped_features} feature(s) had invalid geometry"));
     }
     if skipped_parts > 0 {
         messages.push(format!("{skipped_parts} polygon part(s) could not be projected onto the raster grid; valid parts were retained"));
@@ -2542,13 +2532,22 @@ mod tests {
         assert!(crate::run("kriging", &m, &serde_json::json!({"input": "pts", "field": "v", "range": -5})).is_err());
         assert!(crate::run("kriging", &m, &serde_json::json!({"input": "pts", "field": "v", "range": "none"})).is_err());
 
-        // Projected layer with non-degrees CRS returns error
+        // Source CRS metadata does not change the normalized longitude/latitude geometry.
         let Some(pts_layer) = m.0.get("pts") else { panic!("missing pts") };
         let mut proj_pts = pts_layer.clone();
         proj_pts.crs = Some(rustgis_data::raster::Crs::WebMercator);
         let m_proj = Map(HashMap::from([("proj_pts".to_string(), proj_pts)]));
-        let err = crate::run("kriging", &m_proj, &serde_json::json!({"input": "proj_pts", "field": "v"})).unwrap_err();
-        assert!(err.contains("degrees"));
+        let projected =
+            crate::run("kriging", &m_proj, &serde_json::json!({"input": "proj_pts", "field": "v", "extent": "0,0,2,2", "pixelSize": 0.5}))
+                .unwrap();
+        let Output::Layers(projected_layers) = projected.output else { panic!("expected layers") };
+        assert_eq!(projected_layers.len(), out.len());
+        for (actual, expected) in projected_layers.iter().zip(&out) {
+            let (OutputLayer::Raster { bytes: actual, .. }, OutputLayer::Raster { bytes: expected, .. }) = (actual, expected) else {
+                panic!("expected raster layers")
+            };
+            assert_eq!(actual, expected);
+        }
 
         // Datasets with fewer than two points are cleanly rejected
         let mut single_pt = Dataset::new("single", GeometryType::Point, vec![Field::new("v", FieldType::Double)]);
@@ -3437,10 +3436,9 @@ mod tests {
         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();
+        // 10. A loaded projected layer keeps its source CRS, but its geometry is lon/lat.
+        let p0 = Coord::new(1.0, 7.0);
+        let p1 = Coord::new(3.0, 9.0);
         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);
@@ -3468,19 +3466,20 @@ mod tests {
         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
+        // 11. Source CRS metadata does not change already-normalized geometry.
         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());
+        let out = layer(crate::run("zonalStats", &m_bad_crs, &serde_json::json!({"raster": "elev", "polygons": "bad_crs"})).unwrap());
+        assert_eq!(out.value(0, out.field_index("elev_count").unwrap()), &Value::Int(4));
 
-        // 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();
+        // 12. An invalid exterior does not promote a valid hole.
+        let hole_p0 = Coord::new(1.0, 8.0);
+        let hole_p1 = Coord::new(2.0, 9.0);
         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
@@ -3500,8 +3499,7 @@ mod tests {
         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).
+        // 13. Reject a non-finite exterior vertex without silently dropping it.
         let nan_ext_geom = Geometry::Polygon(vec![vec![vec![
             Coord::new(f64::NAN, p0.y),
             Coord::new(p1.x, p0.y),
@@ -3524,11 +3522,11 @@ mod tests {
         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();
+        // 14. A loaded projected polygon retains its hole.
+        let merc_outer0 = Coord::new(0.0, 7.0);
+        let merc_outer1 = Coord::new(3.0, 10.0);
+        let merc_hole0 = Coord::new(1.0, 8.0);
+        let merc_hole1 = Coord::new(2.0, 9.0);
         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],

crates/data/src/lib.rs

@@ -247,6 +247,8 @@ impl<'de> Deserialize<'de> for Features {
 }
 
 /// A feature class: a named table of features sharing a schema and geometry type.
+/// Geometry coordinates are always WGS 84 longitude/latitude, including datasets created by
+/// tools and plugins. Projected source coordinates must be normalized before insertion.
 #[derive(Clone, Debug, PartialEq, Deserialize)]
 pub struct Dataset {
     pub name: String,

crates/engine/src/cmd/imagery.rs

@@ -530,7 +530,7 @@ fn index(s: &mut Session, p: &Value) -> CmdResult {
 ///     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.
+///   - Every polygon geometry (including multipart geometries) is already 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
@@ -539,7 +539,8 @@ fn index(s: &mut Session, p: &Value) -> CmdResult {
 ///   - The raster is cropped to the tight bounding box of valid selected polygons, and cells outside the
 ///     polygons receive NoData.
 /// - **Partial Failure Reporting**:
-///   - Features with invalid geometry or unprojectable coordinates are skipped and reported in the result messages.
+///   - Features with invalid geometry are skipped and reported in the result messages.
+///   - Polygon parts that cannot be projected onto the raster grid are reported separately.
 ///   - 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)?;
@@ -572,18 +573,12 @@ fn process(s: &mut Session, p: &Value) -> CmdResult {
             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()) {
+                    Some(Geometry::Polygon(polys)) => match rasterops::transform_feature_polygons_to_lonlat(polys, None) {
                         Ok(transformed) => {
                             polygons.extend(transformed);
                         }
@@ -604,7 +599,7 @@ fn process(s: &mut Session, p: &Value) -> CmdResult {
                 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());
+                log::warn!("{pl}: {} feature(s) failed polygon validation during clipping", failed_features.len());
             }
             let mut failed_parts = Vec::new();
             let (cols, rows, mask) = rasterops::polygon_mask_with_failures(&r, &polygons, &mut failed_parts).ok_or_else(|| {
@@ -614,10 +609,7 @@ fn process(s: &mut Session, p: &Value) -> CmdResult {
             let mut result = add(s, out, false)?;
             let mut messages = Vec::new();
             if !failed_features.is_empty() {
-                messages.push(format!(
-                    "{pl}: skipped {} feature(s) with invalid polygon geometry or failed coordinate transformation",
-                    failed_features.len()
-                ));
+                messages.push(format!("{pl}: skipped {} feature(s) with invalid polygon geometry", failed_features.len()));
             }
             if !failed_parts.is_empty() {
                 messages.push(format!(

crates/engine/src/tests.rs

@@ -15207,7 +15207,7 @@ fn dem_tiles_mosaic_clip_by_polygon_share_a_stretch_and_report_overviews() {
     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.
+    // A polygon layer with UTM source metadata still stores lon/lat geometry.
     let mut proj_ds = rustgis_data::Dataset::new("CityProj", rustgis_data::GeometryType::Polygon, vec![]);
     proj_ds.crs = Some(utm);
     let utm_ring = [
@@ -15216,7 +15216,8 @@ fn dem_tiles_mosaic_clip_by_polygon_share_a_stretch_and_report_overviews() {
         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![]);
+    let lonlat_ring = utm_ring.into_iter().map(|c| tr.to_lonlat(c).unwrap()).collect();
+    proj_ds.push(Some(rustgis_data::Geometry::Polygon(vec![vec![lonlat_ring]])), 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);
@@ -16178,6 +16179,71 @@ fn scans_have_no_coordinates_to_give_away() {
     let _ = std::fs::remove_dir_all(&dir);
 }
 
+#[test]
+fn raster_tools_read_projected_zone_files_without_reprojecting_twice() {
+    use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster};
+    use rustgis_data::{Dataset, Geometry, GeometryType, Value};
+    let dir = std::env::temp_dir().join(format!("rustgis-zonal-crs-{}", std::process::id()));
+    std::fs::create_dir_all(&dir).unwrap();
+    let mut s = Session::new(Project::default());
+    let gt = GeoTransform { origin_x: 622_000.0, origin_y: 5_164_000.0, pixel_width: 100.0, pixel_height: -100.0 };
+    let tr = Crs::Epsg(25832).transformer().unwrap();
+    let mut zones = Dataset::new("zones", GeometryType::Polygon, vec![]);
+    let ring = [(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]
+        .into_iter()
+        .map(|(x, y)| tr.to_lonlat(gt.to_world(x, y)).unwrap())
+        .collect();
+    zones.push(Some(Geometry::Polygon(vec![vec![ring]])), vec![]);
+    s.map_mut().add_layer(rustgis_map::Layer::new(0, zones));
+    let nw = tr.to_lonlat(gt.to_world(0.0, 0.0)).unwrap();
+    let se = tr.to_lonlat(gt.to_world(4.0, 4.0)).unwrap();
+    let geographic = GeoTransform { origin_x: nw.x, origin_y: nw.y, pixel_width: (se.x - nw.x) / 4.0, pixel_height: (se.y - nw.y) / 4.0 };
+    for (name, transform, crs) in [("utm", gt, Crs::Epsg(25832)), ("geographic", geographic, Crs::LonLat)] {
+        let raster =
+            Raster::new(name, 4, 4, vec![(1..=16).map(|v| v as f32).collect()], transform, crs, None, ColorHint::Gray, "Float32").unwrap();
+        s.map_mut().add_raster(name, rustgis_map::RasterLayer::from_path(raster, format!("{name}.tif")));
+    }
+    for epsg in [4326, 25832, 32632, 3003, 3857] {
+        for ext in ["gpkg", "shp"] {
+            let path = dir.join(format!("zones_{epsg}.{ext}")).display().to_string();
+            s.run("share.exportLayer", &json!({"layer": "zones", "path": path, "crs": format!("EPSG:{epsg}")})).unwrap();
+            let loaded = s.run("layer.addData", &json!({"path": path, "name": "loaded_zones"})).unwrap()["layer"].clone();
+            let input = s.map().layer(loaded.as_u64().unwrap() as u32).unwrap().dataset.clone();
+            assert_eq!(input.crs.map(|c| c.label()).unwrap_or_else(|| "EPSG:4326".into()), format!("EPSG:{epsg}"), "{ext}");
+            assert!(input.extent().contains(tr.to_lonlat(gt.to_world(2.0, 2.0)).unwrap()), "{epsg}, {ext}: {:?}", input.extent());
+            for raster in ["utm", "geographic"] {
+                let depth = s.undo_depth();
+                s.run(
+                    "gp.run",
+                    &json!({"tool": "zonalStats", "params": {"raster": raster, "polygons": "loaded_zones", "output": "stats"}}),
+                )
+                .unwrap();
+                let out = &s.map().layer_by_name("stats").unwrap().dataset;
+                for (stat, expected) in [("count", Value::Int(4)), ("mean", Value::Float(8.5)), ("sum", Value::Float(34.0))] {
+                    let field = out.field_index(&format!("{raster}_{stat}")).unwrap();
+                    assert_eq!(out.value(0, field), &expected, "EPSG:{epsg}, {ext}, {raster}, {stat}");
+                }
+                assert_eq!(out.features[0].geometry, input.features[0].geometry);
+                assert_eq!(out.crs, input.crs);
+                assert_eq!(s.undo_depth(), depth + 1);
+                assert!(s.undo());
+                assert!(s.map().layer_by_name("stats").is_none());
+                assert!(s.redo());
+                assert!(s.map().layer_by_name("stats").is_some());
+                assert!(s.undo());
+                s.run("raster.process", &json!({"layer": raster, "operation": "clip", "polygon": "loaded_zones", "name": "clipped"}))
+                    .unwrap();
+                let clipped = s.map().layer_by_name("clipped").unwrap().raster.as_ref().unwrap().data.as_ref().unwrap();
+                assert_eq!((clipped.width, clipped.height, clipped.bands[0].stats.count), (2, 2, 4), "{epsg}, {ext}, {raster}");
+                assert_eq!(clipped.bands[0].data.iter().sum::<f32>(), 34.0, "{epsg}, {ext}, {raster}");
+                assert!(s.undo());
+            }
+            s.run("layer.remove", &json!({"layer": loaded})).unwrap();
+        }
+    }
+    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};
@@ -16278,13 +16344,10 @@ fn gp_zonal_stats_end_to_end_and_undo() {
     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();
+    // Loaded WebMercator data already uses longitude/latitude geometry.
     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())]);
+    merc_ds.push(Some(sq(1.0, 7.0, 3.0, 9.0)), vec![Value::Text("M1".into())]);
     s.map_mut().add_layer(rustgis_map::Layer::new(0, merc_ds));
 
     let merc_res = s.run(