Add M-dim bindings (#600)

* Added ISO wkb & wkt * added ZM functions
georust · Jan 10, 2025 · 7bbccb7 · 7bbccb7
1 parent efb3ab7
commit 7bbccb7
Show file tree

Hide file tree

Showing 6 changed files with 251 additions and 13 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -16,6 +16,7 @@
   - Drop `LayerAccess::create_feature_fields` ([#581](https://github.com/georust/gdal/pull/581))
   - Drop `Feature::geometry_by_name` ([#594](https://github.com/georust/gdal/pull/594))
   - Update `SpatialRef::auth_name`, `SpatialRef::name`, `SpatialRef::angular_units_name`, `SpatialRef::linear_units_name` to return `Option<String>` instead of `Result<String>` ([#589](https://github.com/georust/gdal/pull/589))
+  - Update `Geometry::get_point_vec` to modify a `&mut Vec` as opposed to new allocations, returning an `i32` of points added ([#600](https://github.com/georust/gdal/pull/600))
 
 ### Added
 
@@ -28,10 +29,13 @@
   - Add `Defn::field_index` and `Feature::field_index` ([#581](https://github.com/georust/gdal/pull/581))
   - Add `Defn::geometry_field_index` and `Feature::geometry_field_index` ([#594](https://github.com/georust/gdal/pull/594))
   - Add `Dataset::has_capability` for dataset capability check ([#581](https://github.com/georust/gdal/pull/585))
+  - Add methods `add_point_zm`, `add_point_m`, `set_point_zm`, `set_point_m`, `get_point_zm`, `get_point_vec_zm`, `iso_wkt` and `iso_wkb` to `Geometry` ([#600](https://github.com/georust/gdal/pull/600))
+  - Add functions `geometry_type_flatten`, `geometry_type_set_z`, `geometry_type_set_m`, `geometry_type_set_modifier`, `geometry_type_has_z` and `geometry_type_has_m` to `vector::geometry` ([#600](https://github.com/georust/gdal/pull/600))
 
 ### Fixed
 
   - Fix conversion from `ndarray` when the data is offsetted from the start of the buffer ([#569](https://github.com/georust/gdal/pull/569))
+  - use ISO WKT for the `Debug` implementation of `Geometry`, in order to properly display measure values,  ([#600](https://github.com/georust/gdal/pull/600))
 
 ### Removed
 

diff --git a/src/vector/geometry.rs b/src/vector/geometry.rs
@@ -114,6 +114,33 @@ impl Geometry {
         };
     }
 
+    pub fn set_point_zm(&mut self, i: usize, point: (f64, f64, f64, f64)) {
+        let (x, y, z, m) = point;
+        unsafe {
+            gdal_sys::OGR_G_SetPointZM(
+                self.c_geometry(),
+                i as c_int,
+                x as c_double,
+                y as c_double,
+                z as c_double,
+                m as c_double,
+            );
+        };
+    }
+
+    pub fn set_point_m(&mut self, i: usize, point: (f64, f64, f64)) {
+        let (x, y, m) = point;
+        unsafe {
+            gdal_sys::OGR_G_SetPointM(
+                self.c_geometry(),
+                i as c_int,
+                x as c_double,
+                y as c_double,
+                m as c_double,
+            );
+        };
+    }
+
     pub fn set_point_2d(&mut self, i: usize, p: (f64, f64)) {
         let (x, y) = p;
         unsafe {
@@ -138,6 +165,31 @@ impl Geometry {
         unsafe { gdal_sys::OGR_G_AddPoint_2D(self.c_geometry(), x as c_double, y as c_double) };
     }
 
+    pub fn add_point_zm(&mut self, p: (f64, f64, f64, f64)) {
+        let (x, y, z, m) = p;
+        unsafe {
+            gdal_sys::OGR_G_AddPointZM(
+                self.c_geometry(),
+                x as c_double,
+                y as c_double,
+                z as c_double,
+                m as c_double,
+            )
+        };
+    }
+
+    pub fn add_point_m(&mut self, p: (f64, f64, f64)) {
+        let (x, y, m) = p;
+        unsafe {
+            gdal_sys::OGR_G_AddPointM(
+                self.c_geometry(),
+                x as c_double,
+                y as c_double,
+                m as c_double,
+            )
+        };
+    }
+
     /// Get point coordinates from a line string or a point geometry.
     ///
     /// `index` is the line string vertex index, from 0 to `point_count()-1`, or `0` when a point.
@@ -151,9 +203,42 @@ impl Geometry {
         (x, y, z)
     }
 
-    pub fn get_point_vec(&self) -> Vec<(f64, f64, f64)> {
+    /// Get point coordinates from a line string or a point geometry.
+    ///
+    /// `index` is the line string vertex index, from 0 to `point_count()-1`, or `0` when a point.
+    ///
+    /// See: [`OGR_G_GetPointZM`](https://gdal.org/en/stable/api/vector_c_api.html#_CPPv416OGR_G_GetPointZM12OGRGeometryHiPdPdPdPd)
+    pub fn get_point_zm(&self, index: i32) -> (f64, f64, f64, f64) {
+        let mut x: c_double = 0.;
+        let mut y: c_double = 0.;
+        let mut z: c_double = 0.;
+        let mut m: c_double = 0.;
+        unsafe { gdal_sys::OGR_G_GetPointZM(self.c_geometry(), index, &mut x, &mut y, &mut z, &mut m) };
+        (x, y, z, m)
+    }
+
+    /// Appends all points of a line string to `out_points`. 
+    ///
+    /// Only wkbPoint[X], wkbLineString[X] or wkbCircularString[X] may alter `out_points`. Other geometry types will silently do nothing, see
+    /// [`OGR_G_GetPointCount`](https://gdal.org/en/stable/api/vector_c_api.html#_CPPv419OGR_G_GetPointCount12OGRGeometryH)
+    pub fn get_point_vec(&self, out_points: &mut Vec<(f64, f64, f64)>) -> usize {
+        // Consider replacing logic with
+        // [OGR_G_GetPoints](https://gdal.org/en/stable/api/vector_c_api.html#_CPPv415OGR_G_GetPoints12OGRGeometryHPviPviPvi)
         let length = unsafe { gdal_sys::OGR_G_GetPointCount(self.c_geometry()) };
-        (0..length).map(|i| self.get_point(i)).collect()
+        out_points.extend((0..length).map(|i| self.get_point(i)));
+        length as usize
+    }
+
+    /// Appends all points of a line string to `out_points`. 
+    ///
+    /// Only wkbPoint[X], wkbLineString[X] or wkbCircularString[X] may alter `out_points`. Other geometry types will silently do nothing, see
+    /// [`OGR_G_GetPointCount`](https://gdal.org/en/stable/api/vector_c_api.html#_CPPv419OGR_G_GetPointCount12OGRGeometryH)
+    pub fn get_point_vec_zm(&self, out_points: &mut Vec<(f64, f64, f64, f64)>) -> usize {
+        // Consider replacing logic with
+        // [OGR_G_GetPoints](https://gdal.org/en/stable/api/vector_c_api.html#_CPPv415OGR_G_GetPoints12OGRGeometryHPviPviPvi)
+        let length = unsafe { gdal_sys::OGR_G_GetPointCount(self.c_geometry()) };
+        out_points.extend((0..length).map(|i| self.get_point_zm(i)));
+        length as usize
     }
 
     /// Get the geometry type ordinal
@@ -337,7 +422,9 @@ impl Clone for Geometry {
 
 impl Debug for Geometry {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
-        match self.wkt() {
+        let wkt = self.iso_wkt();
+
+        match wkt {
             Ok(wkt) => f.write_str(wkt.as_str()),
             Err(_) => Err(fmt::Error),
         }
@@ -357,6 +444,37 @@ pub fn geometry_type_to_name(ty: OGRwkbGeometryType::Type) -> String {
     _string(rv).unwrap_or_default()
 }
 
+
+/// Returns the 2D geometry type corresponding to the passed geometry type. 
+pub fn geometry_type_flatten(ty: OGRwkbGeometryType::Type) -> OGRwkbGeometryType::Type {
+    unsafe { gdal_sys::OGR_GT_Flatten(ty) }
+}
+
+/// Returns the 3D geometry type corresponding to the passed geometry type. 
+pub fn geometry_type_set_z(ty: OGRwkbGeometryType::Type) -> OGRwkbGeometryType::Type {
+    unsafe { gdal_sys::OGR_GT_SetZ(ty) }
+}
+
+/// Returns the measured geometry type corresponding to the passed geometry type.
+pub fn geometry_type_set_m(ty: OGRwkbGeometryType::Type) -> OGRwkbGeometryType::Type {
+    unsafe { gdal_sys::OGR_GT_SetM(ty) }
+}
+
+/// Returns a XY, XYZ, XYM or XYZM geometry type depending on parameter.
+pub fn geometry_type_set_modifier(ty: OGRwkbGeometryType::Type, set_z: bool, set_m: bool) -> OGRwkbGeometryType::Type {
+    unsafe { gdal_sys::OGR_GT_SetModifier(ty, set_z as i32, set_m as i32) }
+}
+
+/// Returns `true` if the geometry type is a 3D geometry type. 
+pub fn geometry_type_has_z(ty: OGRwkbGeometryType::Type) -> bool {
+    unsafe { gdal_sys::OGR_GT_HasZ(ty) != 0 }
+}
+
+/// Returns `true` if the geometry type is a measured type.
+pub fn geometry_type_has_m(ty: OGRwkbGeometryType::Type) -> bool {
+    unsafe { gdal_sys::OGR_GT_HasM(ty) != 0 }
+}
+
 /// Reference to owned geometry
 pub struct GeometryRef<'a> {
     geom: Geometry,
@@ -385,6 +503,8 @@ impl Debug for GeometryRef<'_> {
 
 #[cfg(test)]
 mod tests {
+    use self::OGRwkbGeometryType::{wkbLineStringZM, wkbMultiPointZM, wkbPointZM};
+
     use super::*;
     use crate::spatial_ref::SpatialRef;
     use crate::test_utils::SuppressGDALErrorLog;
@@ -482,6 +602,23 @@ mod tests {
         assert_eq!(geom, expected);
     }
 
+    #[test]
+    pub fn test_create_multipoint_zm() {
+        let mut geom = Geometry::empty(wkbMultiPointZM).unwrap();
+        let mut point = Geometry::empty(wkbPointZM).unwrap();
+        point.add_point_zm((1.0, 2.0, 3.0, 0.0));
+        geom.add_geometry(point).unwrap();
+        let mut point = Geometry::empty(wkbPointZM).unwrap();
+        point.add_point_zm((3.0, 2.0, 1.0, 1.0));
+        assert!(!point.is_empty());
+        point.set_point_zm(0, (4.0, 2.0, 1.0, 1.0));
+        geom.add_geometry(point).unwrap();
+        assert!(!geom.is_empty());
+        let expected = Geometry::from_wkt("MULTIPOINT ZM ((1.0 2.0 3.0 0.0), (4.0 2.0 1.0 1.0))").unwrap();
+        assert_eq!(geom, expected);
+        assert_eq!(geometry_type_has_m(geom.geometry_type()), geometry_type_has_m(expected.geometry_type()))
+    }
+
     #[test]
     pub fn test_spatial_ref() {
         let geom = Geometry::empty(wkbMultiPolygon).unwrap();
@@ -509,18 +646,24 @@ mod tests {
         ring.add_point_2d((1218405.0658121984, 721108.1805541387));
         ring.add_point_2d((1179091.1646903288, 712782.8838459781));
         assert!(!ring.is_empty());
-        assert_eq!(ring.get_point_vec().len(), 6);
+        let mut ring_vec: Vec<(f64, f64, f64)> = Vec::new();
+        ring.get_point_vec(&mut ring_vec);
+        assert_eq!(ring_vec.len(), 6);
         let mut poly = Geometry::empty(wkbPolygon).unwrap();
         poly.add_geometry(ring.to_owned()).unwrap();
+        let mut poly_vec: Vec<(f64, f64, f64)> = Vec::new();
+        poly.get_point_vec(&mut poly_vec);
         // Points are in ring, not containing geometry.
         // NB: In Python SWIG bindings, `GetPoints` is fallible.
-        assert!(poly.get_point_vec().is_empty());
+        assert!(poly_vec.is_empty());
         assert_eq!(poly.geometry_count(), 1);
         let ring_out = poly.get_geometry(0);
+        let mut ring_out_vec: Vec<(f64, f64, f64)> = Vec::new();
+        ring_out.get_point_vec(&mut ring_out_vec);
         // NB: `wkb()` shows it to be a `LINEARRING`, but returned type is LineString
         assert_eq!(ring_out.geometry_type(), wkbLineString);
         assert!(!&ring_out.is_empty());
-        assert_eq!(ring.get_point_vec(), ring_out.get_point_vec());
+        assert_eq!(ring_vec, ring_out_vec);
     }
 
     #[test]
@@ -532,10 +675,23 @@ mod tests {
         assert_eq!(geom.geometry_type(), OGRwkbGeometryType::wkbPolygon);
         assert!(geom.json().unwrap().contains("Polygon"));
         let inner = geom.get_geometry(0);
-        let points = inner.get_point_vec();
+        let mut points: Vec<(f64, f64, f64)> = Vec::new();
+        inner.get_point_vec(&mut points);
         assert!(!points.is_empty());
     }
 
+    #[test]
+    fn test_get_points_zm() {
+        let mut line = Geometry::empty(wkbLineStringZM).unwrap();
+        line.add_point_zm((0.0, 0.0, 0.0, 0.0));
+        line.add_point_zm((1.0, 0.0, 0.25, 0.5));
+        line.add_point_zm((1.0, 1.0, 0.5, 1.0));
+        let mut line_points: Vec<(f64, f64, f64, f64)> = Vec::new();
+        line.get_point_vec_zm(&mut line_points);
+        assert_eq!(line_points.len(), 3);
+        assert_eq!(line_points.get(2), Some(&(1.0, 1.0, 0.5, 1.0)));
+    }
+
     #[test]
     pub fn test_geometry_type_to_name() {
         assert_eq!(geometry_type_to_name(wkbLineString), "Line String");
@@ -558,4 +714,30 @@ mod tests {
             polygon.wkt().unwrap()
         );
     }
+
+    #[test]
+    fn test_geometry_type_modification() {
+        let mut geom_type = OGRwkbGeometryType::wkbPoint;
+        geom_type = geometry_type_set_z(geom_type);
+        assert_eq!(geom_type, OGRwkbGeometryType::wkbPoint25D);
+        geom_type = geometry_type_set_m(geom_type);
+        assert_eq!(geom_type, OGRwkbGeometryType::wkbPointZM);
+        geom_type = geometry_type_set_modifier(geom_type, false, true);
+        assert_eq!(geom_type, OGRwkbGeometryType::wkbPointM);
+        geom_type = geometry_type_flatten(geom_type);
+        assert_eq!(geom_type, OGRwkbGeometryType::wkbPoint);
+    }
+
+    #[test]
+    fn test_geometry_type_has_zm() {
+        let geom = Geometry::from_wkt("POINT(0 1)").unwrap();
+        assert_eq!(geometry_type_has_z(geom.geometry_type()), false);
+        assert_eq!(geometry_type_has_m(geom.geometry_type()), false);
+        let geom = Geometry::from_wkt("POINT(0 1 2)").unwrap();
+        assert_eq!(geometry_type_has_z(geom.geometry_type()), true);
+        assert_eq!(geometry_type_has_m(geom.geometry_type()), false);
+        let geom = Geometry::from_wkt("POINT ZM (0 1 2 3)").unwrap();
+        assert_eq!(geometry_type_has_z(geom.geometry_type()), true);
+        assert_eq!(geometry_type_has_m(geom.geometry_type()), true);
+    }
 }
diff --git a/src/vector/layer.rs b/src/vector/layer.rs
@@ -1270,7 +1270,8 @@ mod tests {
         with_feature("roads.geojson", 236194095, |feature| {
             let geom = feature.geometry().unwrap();
             assert_eq!(geom.geometry_type(), OGRwkbGeometryType::wkbLineString);
-            let coords = geom.get_point_vec();
+            let mut coords: Vec<(f64, f64, f64)> = Vec::new();
+            geom.get_point_vec(&mut coords);
             assert_eq!(
                 coords,
                 [

diff --git a/src/vector/mod.rs b/src/vector/mod.rs
@@ -79,7 +79,10 @@ pub use feature::{
     OwnedFeatureIterator,
 };
 pub use gdal_sys::{OGRFieldType, OGRwkbGeometryType};
-pub use geometry::{geometry_type_to_name, Geometry};
+pub use geometry::{
+    geometry_type_to_name, geometry_type_flatten, geometry_type_set_z, geometry_type_set_m, 
+    geometry_type_set_modifier, geometry_type_has_z, geometry_type_has_m, Geometry
+};
 pub use layer::{FieldDefn, Layer, LayerAccess, LayerCaps, LayerIterator, OwnedLayer};
 pub use options::LayerOptions;
 pub use transaction::Transaction;

diff --git a/src/vector/ops/conversions/formats.rs b/src/vector/ops/conversions/formats.rs
@@ -97,13 +97,30 @@ impl Geometry {
         Ok(wkt)
     }
 
+    /// Serialize the geometry as SFSQL 1.2 / ISO SQL / MM Part 3 WKT.
+    pub fn iso_wkt(&self) -> Result<String> {
+        let mut c_wkt = null_mut();
+        let rv = unsafe { gdal_sys::OGR_G_ExportToIsoWkt(self.c_geometry(), &mut c_wkt) };
+        if rv != OGRErr::OGRERR_NONE {
+            return Err(GdalError::OgrError {
+                err: rv,
+                method_name: "OGR_G_ExportToIsoWkt",
+            });
+        }
+        let wkt = _string(c_wkt).unwrap_or_default();
+        unsafe { gdal_sys::VSIFree(c_wkt as *mut c_void) };
+        Ok(wkt)
+    }
+
     /// Serializes the geometry to
     /// [WKB](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry#Well-known_binary)
     /// (Well-Known Binary) format.
     pub fn wkb(&self) -> Result<Vec<u8>> {
         let wkb_size = unsafe { gdal_sys::OGR_G_WkbSize(self.c_geometry()) as usize };
         // We default to little-endian for now. A WKB string explicitly indicates the byte
         // order, so this is not a problem for interoperability.
+        // 
+        // Consider using `Vec::MaybeUninit` in future.
         let byte_order = gdal_sys::OGRwkbByteOrder::wkbNDR;
         let mut wkb = vec![0; wkb_size];
         let rv =
@@ -117,6 +134,28 @@ impl Geometry {
         Ok(wkb)
     }
 
+    /// Serializes the geometry to SFSQL 1.2 / ISO SQL / MM Part 3
+    /// [WKB](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry#Well-known_binary)
+    /// (Well-Known Binary) format.
+    pub fn iso_wkb(&self) -> Result<Vec<u8>> {
+        let wkb_size = unsafe { gdal_sys::OGR_G_WkbSize(self.c_geometry()) as usize };
+        // We default to little-endian for now. A WKB string explicitly indicates the byte
+        // order, so this is not a problem for interoperability.
+        //
+        // Consider using `Vec::MaybeUninit` in future.
+        let byte_order = gdal_sys::OGRwkbByteOrder::wkbNDR;
+        let mut wkb = vec![0; wkb_size];
+        let rv =
+            unsafe { gdal_sys::OGR_G_ExportToIsoWkb(self.c_geometry(), byte_order, wkb.as_mut_ptr()) };
+        if rv != gdal_sys::OGRErr::OGRERR_NONE {
+            return Err(GdalError::OgrError {
+                err: rv,
+                method_name: "OGR_G_ExportToIsoWkb",
+            });
+        }
+        Ok(wkb)
+    }
+
     /// Serialize the geometry as GeoJSON.
     ///
     /// See: [`OGR_G_ExportToJson`](https://gdal.org/api/vector_c_api.html#_CPPv418OGR_G_ExportToJson12OGRGeometryH)
@@ -141,6 +180,15 @@ mod tests {
         assert_eq!(new_geom, orig_geom);
     }
 
+    #[test]
+    pub fn test_iso_wkb() {
+        let wkt = "POLYGON ZM ((45.0 45.0 0.0 0.0, 45.0 50.0 0.0 0.25, 50.0 50.0 0.0 0.50, 50.0 45.0 0.0 0.75, 45.0 45.0 0.0 1.0))";
+        let orig_geom = Geometry::from_wkt(wkt).unwrap();
+        let wkb = orig_geom.iso_wkb().unwrap();
+        let new_geom = Geometry::from_wkb(&wkb).unwrap();
+        assert_eq!(new_geom, orig_geom);
+    }
+
     #[test]
     pub fn test_geojson() {
         let json = r#"{ "type": "Point", "coordinates": [10, 20] }"#;