fix(analysis): Minimum Bounding Geometry and Kernel Density say what they leave out

crates/data/src/export.rs

@@ -116,6 +116,7 @@ fn write_chunks_to(
     let left_out = match ext.as_str() {
         #[cfg(not(target_arch = "wasm32"))]
         "gpkg" => {
+            // Every feature is written: one of another kind makes the column GEOMETRY.
             table = Some(crate::source::gpkg::write_features_in(&layer.name, geometry_type, &layer.fields, path, crs, gpkg, &mut feed)?);
             0
         }

crates/data/src/export/tests.rs

@@ -370,3 +370,138 @@ fn geoparquet_counts_values_written_as_null() {
     assert_eq!(note, format!("{want}, '{}…'", "x".repeat(40)));
     let _ = std::fs::remove_dir_all(d);
 }
+
+/// A streamed GeoPackage export of a layer holding another kind of geometry than its own writes
+/// every feature in a GEOMETRY column, as the export of the whole layer does, even when the
+/// other kind comes only in a late chunk (the export starts again for it): squares with points
+/// among them, a few features a chunk, give the same table.
+#[test]
+fn a_streamed_geopackage_of_a_mixed_layer_keeps_every_kind() {
+    let d = folder("mixed-kinds");
+    let mut mixed = squares();
+    mixed.name = "mixed".into();
+    for i in [0usize, 5, 6, 299] {
+        if let Some(f) = mixed.features.get_mut(i) {
+            f.geometry = Some(Geometry::Point(Coord::new(i as f64, 1.0)));
+        }
+    }
+    let layer = Dataset::new(mixed.name.clone(), mixed.geometry_type, mixed.fields.clone());
+    let (a, b) = (d.join("streamed.gpkg"), d.join("whole.gpkg"));
+    let streamed = write_chunks_in(&Runs::new(&mixed, 4), &layer, &a, Crs::LonLat, 1).unwrap();
+    assert_eq!((streamed.left_out, streamed.written), (0, 301));
+    assert_eq!(crate::io::write_path(&mixed, &b).unwrap(), 0);
+    let table = |p: &Path| {
+        let c = rusqlite::Connection::open(p).unwrap();
+        let kind = c.query_row("SELECT geometry_type_name FROM gpkg_geometry_columns", [], |r| r.get::<_, String>(0)).unwrap();
+        let blobs = c
+            .query_row("SELECT count(*) FROM mixed WHERE geom IS NOT NULL AND substr(hex(geom), 1, 4) = '4750'", [], |r| r.get::<_, i64>(0))
+            .unwrap();
+        (kind, blobs)
+    };
+    // Every geometry row is a GeoPackage blob: the 296 squares and the 4 points.
+    assert_eq!(table(&a), ("GEOMETRY".to_string(), 300));
+    assert_eq!(table(&b), ("GEOMETRY".to_string(), 300));
+    assert_eq!(read_back(&a).features, read_back(&b).features);
+    let points = read_back(&a).features.iter().filter(|f| matches!(f.geometry, Some(Geometry::Point(_)))).count();
+    assert_eq!(points, 4);
+    // A multipoint among them comes back a multipoint, with its members, streamed and whole.
+    let pair = Geometry::MultiPoint(vec![Coord::new(1.5, 2.5), Coord::new(3.5, 4.5)]);
+    if let Some(f) = mixed.features.get_mut(10) {
+        f.geometry = Some(pair.clone());
+    }
+    let (c, e) = (d.join("streamed-multi.gpkg"), d.join("whole-multi.gpkg"));
+    write_chunks_in(&Runs::new(&mixed, 4), &layer, &c, Crs::LonLat, 1).unwrap();
+    crate::io::write_path(&mixed, &e).unwrap();
+    for p in [&c, &e] {
+        let back = read_back(p);
+        let multis: Vec<_> =
+            back.features.iter().filter_map(|f| f.geometry.as_ref()).filter(|g| matches!(g, Geometry::MultiPoint(_))).collect();
+        assert_eq!(multis, [&pair], "{}", p.display());
+    }
+}
+
+/// A layer that declares no geometry type but carries geometry exports to the same GeoPackage
+/// streamed as whole: a table of its geometry's kind when it holds one, a GEOMETRY table when it
+/// holds several, nothing left out either way.
+#[test]
+fn a_typeless_layer_exports_the_same_geopackage_streamed_or_whole() {
+    let d = folder("typeless");
+    let mut ds = squares();
+    ds.name = "typeless".into();
+    ds.geometry_type = GeometryType::None;
+    let kind = |p: &Path| {
+        let c = rusqlite::Connection::open(p).unwrap();
+        c.query_row("SELECT geometry_type_name FROM gpkg_geometry_columns", [], |r| r.get::<_, String>(0)).unwrap()
+    };
+    let layer = Dataset::new(ds.name.clone(), ds.geometry_type, ds.fields.clone());
+    for (with_point, want) in [(false, "MULTIPOLYGON"), (true, "GEOMETRY")] {
+        if with_point && let Some(f) = ds.features.get_mut(3) {
+            f.geometry = Some(Geometry::Point(Coord::new(3.0, 1.0)));
+        }
+        let (a, b) = (d.join(format!("streamed-{with_point}.gpkg")), d.join(format!("whole-{with_point}.gpkg")));
+        assert_eq!(write_chunks_in(&Runs::new(&ds, 4), &layer, &a, Crs::LonLat, 1).unwrap().left_out, 0);
+        assert_eq!(crate::io::write_path(&ds, &b).unwrap(), 0);
+        assert_eq!((kind(&a), kind(&b)), (want.to_string(), want.to_string()));
+        assert_eq!(read_back(&a).features, read_back(&b).features);
+        assert_eq!(read_back(&a).features.len(), 300);
+    }
+}
+
+/// A projected streamed GeoPackage export of a mixed layer writes every feature, projected, in a
+/// GEOMETRY column, as the whole export does.
+#[test]
+fn a_projected_streamed_geopackage_of_a_mixed_layer_keeps_every_kind() {
+    let d = folder("mixed-projected");
+    let mut mixed = squares();
+    mixed.name = "mixed".into();
+    for i in [2usize, 7] {
+        if let Some(f) = mixed.features.get_mut(i) {
+            f.geometry = Some(Geometry::Point(Coord::new(3.0, 1.0)));
+        }
+    }
+    let utm = Crs::from_epsg(32631).unwrap();
+    let layer = Dataset::new(mixed.name.clone(), mixed.geometry_type, mixed.fields.clone());
+    let (a, b) = (d.join("streamed.gpkg"), d.join("whole.gpkg"));
+    let streamed = write_chunks_in(&Runs::new(&mixed, 4), &layer, &a, utm, 1).unwrap();
+    assert_eq!((streamed.left_out, streamed.written), (0, 301));
+    assert_eq!(crate::io::write_path_in(&mixed, &b, utm).unwrap(), 0);
+    assert_eq!(read_back(&a).features, read_back(&b).features);
+    assert_eq!(read_back(&a).features.len(), 300);
+}
+
+/// A mixed layer added as a table to an existing GeoPackage, streamed or whole, gets a GEOMETRY
+/// column holding every feature, the file's other table left as it was.
+#[test]
+fn a_mixed_layer_added_to_a_geopackage_keeps_every_kind() {
+    let d = folder("mixed-added");
+    let mut mixed = squares();
+    mixed.name = "mixed".into();
+    for i in [1usize, 298] {
+        if let Some(f) = mixed.features.get_mut(i) {
+            f.geometry = Some(Geometry::Point(Coord::new(i as f64, 1.0)));
+        }
+    }
+    let layer = Dataset::new(mixed.name.clone(), mixed.geometry_type, mixed.fields.clone());
+    let add = GpkgTable { mode: crate::io::GpkgMode::AddTable, name: None };
+    let mut other = squares();
+    other.name = "other".into();
+    other.features.truncate(3);
+    for (i, streamed) in [false, true].into_iter().enumerate() {
+        let path = d.join(format!("county-{i}.gpkg"));
+        crate::io::write_path(&other, &path).unwrap();
+        if streamed {
+            let (e, _) = write_chunks_to(&Runs::new(&mixed, 4), &layer, &path, Crs::LonLat, &add, 1).unwrap();
+            assert_eq!((e.left_out, e.written), (0, 301));
+        } else {
+            assert_eq!(crate::io::write_path_to(&mixed, &path, Crs::LonLat, &add).unwrap().0, 0);
+        }
+        let c = rusqlite::Connection::open(&path).unwrap();
+        let kind = |t: &str| {
+            c.query_row("SELECT geometry_type_name FROM gpkg_geometry_columns WHERE table_name = ?1", [t], |r| r.get::<_, String>(0))
+                .unwrap()
+        };
+        assert_eq!((kind("mixed"), kind("other")), ("GEOMETRY".to_string(), "MULTIPOLYGON".to_string()));
+        let rows: i64 = c.query_row("SELECT count(*) FROM mixed WHERE geom IS NOT NULL", [], |r| r.get(0)).unwrap();
+        assert_eq!(rows, 300);
+    }
+}

crates/data/src/geometry.rs

@@ -334,11 +334,19 @@ impl VertexId {
 /// How many vertices a path has: a closed ring's repeated last point doesn't count.
 fn open_len(path: &[Coord], ring: bool) -> usize {
     match (path.first(), path.last()) {
-        (Some(a), Some(b)) if ring && path.len() > 1 && a == b => path.len() - 1,
+        (Some(a), Some(b)) if ring && path.len() > 1 && same_point(a, b) => path.len() - 1,
         _ => path.len(),
     }
 }
 
+/// Whether two points are the same place. A NaN coordinate equals nothing as a number, so NaN is
+/// matched by its bits: a ring that repeats a NaN point closes on it.
+fn same_point(a: &Coord, b: &Coord) -> bool {
+    // Each coordinate on its own: equal as numbers (so 0 and −0 match), or a NaN by its bits.
+    let same = |p: f64, q: f64| p == q || p.to_bits() == q.to_bits();
+    same(a.x, b.x) && same(a.y, b.y)
+}
+
 impl Geometry {
     /// Calls `f` with each path that vertices lie on: its part and ring, and whether it is a polygon
     /// ring. Nothing is allocated, so a geometry of many parts costs no memory to walk.
@@ -762,6 +770,30 @@ mod tests {
         assert!(polygon_contains(&[square(0.0, 0.0, 10.0), square(4.0, 4.0, 2.0)], lp));
     }
 
+    /// A polygon ring that starts and ends on the same NaN point is closed on it (NaN equals
+    /// nothing as a number, so it's matched by its bits): its last point isn't a vertex of its own,
+    /// so every ring reader sees four vertices and four edges, and deleting the first vertex keeps
+    /// the ring closed on the next one. A line with the same ends isn't a ring: all three count.
+    #[test]
+    fn a_ring_closed_on_a_nan_point_counts_it_once() {
+        let gap = Coord::new(f64::NAN, 0.0);
+        let ring = vec![gap, Coord::new(1.0, 0.0), Coord::new(1.0, 1.0), Coord::new(0.0, 1.0), gap];
+        let g = Geometry::Polygon(vec![vec![ring]]);
+        assert_eq!((g.vertices().len(), g.edges().len()), (4, 4));
+        let Ok(Geometry::Polygon(after)) = g.delete_vertex(VertexId { part: 0, ring: 0, index: 0 }) else { panic!("not deleted") };
+        let ring = after.first().and_then(|rings| rings.first()).cloned().unwrap_or_default();
+        assert_eq!(ring.len(), 4);
+        assert_eq!((ring.first(), ring.last()), (Some(&Coord::new(1.0, 0.0)), Some(&Coord::new(1.0, 0.0))));
+        let line = Geometry::Line(vec![vec![gap, Coord::new(1.0, 0.0), gap]]);
+        assert_eq!(line.vertices().len(), 3);
+        // The same NaN point written once with 0 and once with −0 still closes the ring.
+        let signed = vec![Coord::new(f64::NAN, 0.0), Coord::new(1.0, 0.0), Coord::new(1.0, 1.0), Coord::new(f64::NAN, -0.0)];
+        assert_eq!(Geometry::Polygon(vec![vec![signed]]).vertices().len(), 3);
+        // And the other way round: a signed zero in x, the NaN in y.
+        let mirrored = vec![Coord::new(0.0, f64::NAN), Coord::new(1.0, 0.0), Coord::new(1.0, 1.0), Coord::new(-0.0, f64::NAN)];
+        assert_eq!(Geometry::Polygon(vec![vec![mirrored]]).vertices().len(), 3);
+    }
+
     #[test]
     fn line_hit_and_box() {
         let g = Geometry::Line(vec![vec![Coord::new(0.0, 0.0), Coord::new(10.0, 0.0)]]);

crates/data/src/io.rs

@@ -546,6 +546,8 @@ pub fn write_path_to(ds: &Dataset, path: &Path, crs: Crs, gpkg: &GpkgTable) -> R
             #[cfg(not(target_arch = "wasm32"))]
             {
                 crate::progress::begin(ds.features.len());
+                // Every feature is written: a layer of several kinds of geometry gets a GEOMETRY
+                // column, and one that declares no type takes its geometry's.
                 let table = crate::source::gpkg::write_table_in(ds, path, crs, gpkg)?;
                 Ok((0, Some(table)))
             }
@@ -2516,6 +2518,20 @@ fn dbf_boolean(value: &Value) -> Option<bool> {
     }
 }
 
+/// The notice for the `left_out` features an export to `name` (a file of extension `ext`) left
+/// out, as [`write_path_in`] and [`write_bytes_in`] count them; `None` when none were. A
+/// Shapefile export (`.shp`, or zipped) writes only drawable geometry of the layer's one kind (it
+/// writes no null shapes, though the format has them); FlatGeobuf, KML and GPX exports write only
+/// features with geometry; the other formats leave nothing out.
+pub fn left_out_note(name: &str, ext: &str, left_out: usize) -> Option<String> {
+    (left_out > 0).then(|| match ext {
+        "shp" | "zip" => format!(
+            "{name}: {left_out} features left out of the Shapefile: without drawable geometry, or with geometry of another kind than the layer's"
+        ),
+        _ => format!("{name}: {left_out} features without drawable geometry were left out"),
+    })
+}
+
 /// Whether [`shapefile_parts_in`] writes a feature with geometry `g` to a layer of type `gt`: one
 /// of the layer's kind with something drawable (a finite point, a line part of two finite points,
 /// a polygon whose outer ring has three).
@@ -2682,6 +2698,102 @@ mod tests {
     use super::*;
     use std::io::Write;
 
+    /// A Shapefile export of a polygon layer holding a point and a feature without geometry writes
+    /// the polygon and counts both left out, the point for its kind; a layer of polygons loses none.
+    #[cfg(not(target_arch = "wasm32"))]
+    #[test]
+    fn a_shapefile_counts_geometry_of_another_kind_left_out() -> std::result::Result<(), Box<dyn std::error::Error>> {
+        let geojson = r#"{"type":"FeatureCollection","features":[
+            {"type":"Feature","properties":{"name":"area"},"geometry":{"type":"Polygon","coordinates":[[[0,0],[1,0],[1,1],[0,0]]]}},
+            {"type":"Feature","properties":{"name":"spot"},"geometry":{"type":"Point","coordinates":[5,5]}},
+            {"type":"Feature","properties":{"name":"none"},"geometry":null}]}"#;
+        let (mixed, _) = read_geojson(geojson, "mixed")?;
+        let dir = tempfile::tempdir()?;
+        assert_eq!(write_path(&mixed, &dir.path().join("mixed.shp"))?, 2);
+        let (back, _) = read_path(&dir.path().join("mixed.shp"))?;
+        assert_eq!(back.features.len(), 1);
+        let mut areas = mixed.clone();
+        areas.features.truncate(1);
+        assert_eq!(write_path(&areas, &dir.path().join("areas.shp"))?, 0);
+        Ok(())
+    }
+
+    /// The left-out notice names what each format can't hold: a Shapefile, features without
+    /// drawable geometry or of another kind; the others, features without drawable geometry.
+    /// Nothing left out, no notice.
+    #[test]
+    fn left_out_notes_say_what_the_format_cant_hold() {
+        for ext in ["shp", "zip"] {
+            assert_eq!(
+                left_out_note("out", ext, 2).as_deref(),
+                Some(
+                    "out: 2 features left out of the Shapefile: without drawable geometry, or with geometry of another kind than the layer's"
+                )
+            );
+        }
+        for ext in ["fgb", "kml", "gpx"] {
+            assert_eq!(left_out_note("out", ext, 1).as_deref(), Some("out: 1 features without drawable geometry were left out"));
+        }
+        assert_eq!(left_out_note("out", "shp", 0), None);
+        assert_eq!(left_out_note("out", "gpkg", 0), None);
+    }
+
+    /// A GeoPackage export loses no geometry: a layer read from a GeoJSON of a polygon and a point
+    /// (declared a polygon layer, by its most common kind) gets a GEOMETRY column holding both, and
+    /// a feature without geometry stays, as a table holds it. A layer of one kind keeps its kind's
+    /// column, whether it declares its type or not.
+    #[cfg(not(target_arch = "wasm32"))]
+    #[test]
+    fn a_geopackage_of_a_mixed_layer_keeps_every_kind() -> std::result::Result<(), Box<dyn std::error::Error>> {
+        let geojson = r#"{"type":"FeatureCollection","features":[
+            {"type":"Feature","properties":{"name":"area"},"geometry":{"type":"Polygon","coordinates":[[[0,0],[1,0],[1,1],[0,0]]]}},
+            {"type":"Feature","properties":{"name":"spot"},"geometry":{"type":"Point","coordinates":[5,5]}},
+            {"type":"Feature","properties":{"name":"none"},"geometry":null}]}"#;
+        let (mixed, _) = read_geojson(geojson, "mixed")?;
+        assert_eq!(mixed.geometry_type, GeometryType::Polygon);
+        let dir = tempfile::tempdir()?;
+        // The column's type, and the rows written with whether each has geometry, read with SQL.
+        let table = |path: &Path| -> std::result::Result<(String, Vec<(String, bool)>), Box<dyn std::error::Error>> {
+            let conn = rusqlite::Connection::open(path)?;
+            let kind: String = conn.query_row("SELECT geometry_type_name FROM gpkg_geometry_columns", [], |r| r.get(0))?;
+            let mut rows = conn.prepare(r#"SELECT "name", "geom" IS NOT NULL FROM "mixed" ORDER BY "fid""#)?;
+            let written = rows.query_map([], |r| Ok((r.get(0)?, r.get(1)?)))?.collect::<std::result::Result<_, _>>()?;
+            Ok((kind, written))
+        };
+        let all = vec![("area".to_string(), true), ("spot".to_string(), true), ("none".to_string(), false)];
+        let path = dir.path().join("mixed.gpkg");
+        assert_eq!(write_path(&mixed, &path)?, 0);
+        assert_eq!(table(&path)?, ("GEOMETRY".to_string(), all.clone()));
+        // Read back, the point and the polygon are both there.
+        let (back, _) = read_path(&path)?;
+        let mut kinds: Vec<GeometryType> = back.features.iter().filter_map(|f| f.geometry.as_ref().map(Geometry::kind)).collect();
+        kinds.sort_by_key(|k| format!("{k:?}"));
+        assert_eq!(kinds, [GeometryType::Point, GeometryType::Polygon]);
+        // A layer that declares no geometry type does the same.
+        let mut typeless = mixed.clone();
+        typeless.geometry_type = GeometryType::None;
+        let at = dir.path().join("typeless.gpkg");
+        assert_eq!(write_path(&typeless, &at)?, 0);
+        assert_eq!(table(&at)?, ("GEOMETRY".to_string(), all));
+        // Projected, every feature is written too.
+        let utm = Crs::from_epsg(32631).ok_or("no UTM zone 31N")?;
+        let projected = dir.path().join("projected.gpkg");
+        assert_eq!(write_path_in(&mixed, &projected, utm)?, 0);
+        assert_eq!(table(&projected)?.0, "GEOMETRY");
+        // One kind, declared or not: that kind's column, nothing left out.
+        let mut areas = mixed.clone();
+        areas.features.retain(|f| !matches!(f.geometry, Some(Geometry::Point(_))));
+        let one = vec![("area".to_string(), true), ("none".to_string(), false)];
+        let at = dir.path().join("areas.gpkg");
+        assert_eq!(write_path(&areas, &at)?, 0);
+        assert_eq!(table(&at)?, ("MULTIPOLYGON".to_string(), one.clone()));
+        areas.geometry_type = GeometryType::None;
+        let at = dir.path().join("typeless-areas.gpkg");
+        assert_eq!(write_path(&areas, &at)?, 0);
+        assert_eq!(table(&at)?, ("MULTIPOLYGON".to_string(), one));
+        Ok(())
+    }
+
     /// Exercises caller cleanup through the public Shapefile writer, one fresh process per collision.
     #[cfg(not(target_arch = "wasm32"))]
     #[test]