Merge #366

366: Exposed spatial predicates over `Geometry` r=metasim a=metasim

- [x] I agree to follow the project's [code of conduct](https://github.com/georust/gdal/blob/master/CODE_OF_CONDUCT.md).
- [X] I added an entry to `CHANGES.md` if knowledge of this change could be valuable to users.
---

Partially addresses #306 

cc: `@thomas001`

Co-authored-by: Simeon H.K. Fitch <fitch@astraea.io>

CHANGES.md

@@ -46,6 +46,10 @@
 
   - <https://github.com/georust/gdal/pull/367>
 
+- Exposed spatial predicates over `Geometry`: `intersects`, `contains`, `disjoint`, `touches`, `crosses`, `within`, and `overlaps`.
+
+  - <https://github.com/georust/gdal/pull/366>
+
 ## 0.14
 
 - Added new content to `README.md` and the root docs.

src/dataset.rs

@@ -529,7 +529,7 @@ impl Dataset {
     /// * `overviews` - list of overview decimation factors, e.g. `&[2, 4, 8, 16, 32]`
     /// * `bands` - list of bands to build the overviews for, or empty for all bands
     ///
-    /// [`GDALBuildOverviews`]: https://gdal.org/doxygen/gdal_8h.html#a767f4456a6249594ee18ea53f68b7e80
+    /// [`GDALBuildOverviews`]: https://gdal.org/api/raster_c_api.html#_CPPv418GDALBuildOverviews12GDALDatasetHPKciPKiiPKi16GDALProgressFuncPv
     pub fn build_overviews(
         &mut self,
         resampling: &str,

src/lib.rs

@@ -25,7 +25,7 @@
 //! ## Usage
 //!
 //! This crate provides high-level, idiomatic Rust bindings for GDAL.
-//! To do that, it uses [`gdal-sys`](gdal-sys) internally, a low-level interface to the GDAL C library,
+//! To do that, it uses [`gdal_sys`] internally, a low-level interface to the GDAL C library,
 //! which is generated using [`bindgen`](https://rust-lang.github.io/rust-bindgen/).
 //! Using the `gdal-sys` crate directly is normally not needed, but it can be useful in order to call APIs that have not yet been exposed in this crate.
 //!

src/raster/mdarray.rs

@@ -199,12 +199,12 @@ impl<'a> MDArray<'a> {
         Ok(())
     }
 
-    /// Read a 'Buffer<T>' from this band. T implements 'GdalType'
+    /// Read a [`Vec<T>`] from this band, where `T` implements [`GdalType`].
     ///
     /// # Arguments
-    /// * array_start_index - Values representing the starting index to read in each dimension (in `[0, aoDims[i].GetSize()-1]` range).
+    /// * `array_start_index` - Values representing the starting index to read in each dimension (in `[0, aoDims[i].GetSize()-1]` range).
     ///   Array of `GetDimensionCount()` values. Must not be empty, unless for a zero-dimensional array.
-    /// * count - Values representing the number of values to extract in each dimension. Array of `GetDimensionCount()` values.
+    /// * `count` - Values representing the number of values to extract in each dimension. Array of `GetDimensionCount()` values.
     ///   Must not be empty, unless for a zero-dimensional array.
     ///
     pub fn read_as<T: Copy + GdalType>(

src/raster/rasterband.rs

@@ -222,7 +222,7 @@ impl<'a> RasterBand<'a> {
         (self.x_size(), self.y_size())
     }
 
-    /// Read data from this band into a slice, where `T` implements ['GdalType']
+    /// Read data from this band into a slice, where `T` implements [`GdalType`]
     ///
     /// # Arguments
     ///
@@ -293,7 +293,7 @@ impl<'a> RasterBand<'a> {
         Ok(())
     }
 
-    /// Read a 'Buffer<T>' from this band, where `T` implements ['GdalType'].
+    /// Read a [`Buffer<T>`] from this band, where `T` implements [`GdalType`].
     ///
     /// # Arguments
     ///
@@ -373,7 +373,7 @@ impl<'a> RasterBand<'a> {
 
     #[cfg(feature = "ndarray")]
     #[cfg_attr(docsrs, doc(cfg(feature = "array")))]
-    /// Read a 'Array2<T>' from this band, where `T` implements ['GdalType'].
+    /// Read a [`Array2<T>`] from this band, where `T` implements [`GdalType`].
     ///
     /// # Arguments
     ///
@@ -400,15 +400,15 @@ impl<'a> RasterBand<'a> {
         )?)
     }
 
-    /// Read the full band as a 'Buffer<T>', where `T` implements ['GdalType'].
+    /// Read the full band as a [`Buffer<T>`], where `T` implements [`GdalType`].
     pub fn read_band_as<T: Copy + GdalType>(&self) -> Result<Buffer<T>> {
         let size = self.size();
         self.read_as::<T>((0, 0), (size.0, size.1), (size.0, size.1), None)
     }
 
     #[cfg(feature = "ndarray")]
     #[cfg_attr(docsrs, doc(cfg(feature = "array")))]
-    /// Read a 'Array2<T>' from a 'Dataset' block, where `T` implements ['GdalType'].
+    /// Read a [`Array2<T>`] from a [`Dataset`] block, where `T` implements [`GdalType`].
     ///
     /// # Arguments
     /// * `block_index` - the block index
@@ -440,7 +440,7 @@ impl<'a> RasterBand<'a> {
         Array2::from_shape_vec((size.1, size.0), data).map_err(Into::into)
     }
 
-    /// Write a 'Buffer<T>' into a 'Dataset'.
+    /// Write a [`Buffer<T>`] into a [`Dataset`].
     ///
     /// # Arguments
     ///

src/vector/geometry.rs

@@ -375,7 +375,7 @@ impl Geometry {
         Ok(())
     }
 
-    // Transform the geometry inplace (when we own the Geometry)
+    /// Transform the geometry inplace (when we own the Geometry)
     pub fn transform_inplace(&mut self, htransform: &CoordTransform) -> Result<()> {
         let rv = unsafe { gdal_sys::OGR_G_Transform(self.c_geometry(), htransform.to_c_hct()) };
         if rv != OGRErr::OGRERR_NONE {
@@ -387,7 +387,7 @@ impl Geometry {
         Ok(())
     }
 
-    // Return a new transformed geometry (when the Geometry is owned by a Feature)
+    /// Return a new transformed geometry (when the Geometry is owned by a Feature)
     pub fn transform(&self, htransform: &CoordTransform) -> Result<Geometry> {
         let new_c_geom = unsafe { gdal_sys::OGR_G_Clone(self.c_geometry()) };
         let rv = unsafe { gdal_sys::OGR_G_Transform(new_c_geom, htransform.to_c_hct()) };
@@ -411,7 +411,7 @@ impl Geometry {
         Ok(())
     }
 
-    /// Transforms a geometrys coordinates into another SpatialRef
+    /// Transforms a geometry's coordinates into another SpatialRef
     pub fn transform_to(&self, spatial_ref: &SpatialRef) -> Result<Geometry> {
         let new_c_geom = unsafe { gdal_sys::OGR_G_Clone(self.c_geometry()) };
         let rv = unsafe { gdal_sys::OGR_G_TransformTo(new_c_geom, spatial_ref.to_c_hsrs()) };
@@ -432,7 +432,7 @@ impl Geometry {
     ///
     /// Returns `Some(SpatialRef)`, or `None` if one isn't defined.
     ///
-    /// Refer: [OGR_G_GetSpatialReference](https://gdal.org/doxygen/ogr__api_8h.html#abc393e40282eec3801fb4a4abc9e25bf)
+    /// See: [OGR_G_GetSpatialReference](https://gdal.org/doxygen/ogr__api_8h.html#abc393e40282eec3801fb4a4abc9e25bf)
     pub fn spatial_ref(&self) -> Option<SpatialRef> {
         let c_spatial_ref = unsafe { gdal_sys::OGR_G_GetSpatialReference(self.c_geometry()) };
 
@@ -504,6 +504,20 @@ impl Geometry {
             Ok(unsafe { Geometry::with_c_geometry(c_geom, true) })
         }
     }
+
+    /// Test if the geometry is valid.
+    ///
+    /// This function requires the GEOS library.
+    /// If OGR is built without the GEOS library, this function will always return `false`.
+    /// Check with [`VersionInfo::has_geos`][has_geos].
+    ///
+    /// See: [`OGR_G_IsValid`](https://gdal.org/api/vector_c_api.html#_CPPv413OGR_G_IsValid12OGRGeometryH)
+    ///
+    /// [has_geos]: crate::version::VersionInfo::has_geos
+    pub fn is_valid(&self) -> bool {
+        let p = unsafe { gdal_sys::OGR_G_IsValid(self.c_geometry()) };
+        p != 0
+    }
 }
 
 impl Drop for Geometry {
@@ -689,13 +703,15 @@ mod tests {
         let src = Geometry::from_wkt("POINT (0 0)").unwrap();
         let dst = src.make_valid(&CslStringList::new());
         assert!(dst.is_ok());
+        assert!(dst.unwrap().is_valid());
     }
 
     #[test]
     /// Un-repairable geometry case
     pub fn test_make_valid_invalid() {
         let _nolog = SuppressGDALErrorLog::new();
         let src = Geometry::from_wkt("LINESTRING (0 0)").unwrap();
+        assert!(!src.is_valid());
         let dst = src.make_valid(&CslStringList::new());
         assert!(dst.is_err());
     }
@@ -704,8 +720,10 @@ mod tests {
     /// Repairable case (self-intersecting)
     pub fn test_make_valid_repairable() {
         let src = Geometry::from_wkt("POLYGON ((0 0, 10 10, 0 10, 10 0, 0 0))").unwrap();
+        assert!(!src.is_valid());
         let dst = src.make_valid(&CslStringList::new());
         assert!(dst.is_ok());
+        assert!(dst.unwrap().is_valid());
     }
 
     #[cfg(all(major_ge_3, minor_ge_4))]
@@ -718,5 +736,6 @@ mod tests {
         let opts = CslStringList::try_from(&[("STRUCTURE", "LINEWORK")]).unwrap();
         let dst = src.make_valid(&opts);
         assert!(dst.is_ok(), "{dst:?}");
+        assert!(dst.unwrap().is_valid());
     }
 }

src/vector/mod.rs

@@ -1,5 +1,7 @@
 //! GDAL Vector Data API
 //!
+//! ## Example
+//!
 //! This example opens a vector [`Dataset`](crate::Dataset) and iterates over the various levels of structure within it.
 //! The GDAL vector data model is quite sophisticated, so please refer to the GDAL
 //! [Vector Data Model](https://gdal.org/user/vector_data_model.html) document for specifics.

src/vector/ops/mod.rs

@@ -1,3 +1,4 @@
 mod geometry;
+mod predicates;
 
 pub use geometry::intersection::Intersection as GeometryIntersection;

src/vector/ops/predicates.rs

@@ -0,0 +1,189 @@
+use crate::vector::Geometry;
+
+/// These methods define common [spatial relations](https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates) between
+/// two geometries.
+impl Geometry {
+    /// Tests if two geometries [_intersect_][DE-9IM];
+    /// `self` and `other` have at least one point in common.
+    ///
+    /// If GEOS is enabled, then this is done in rigorous fashion, otherwise `true` is returned
+    /// if the envelopes (bounding boxes) of the two geometries overlap. Check with [`VersionInfo::has_geos`][has_geos].
+    ///
+    /// See: [`OGR_G_Intersects`][OGR_G_Intersects]
+    ///
+    /// [DE-9IM]: https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates
+    /// [OGR_G_Intersects]: https://gdal.org/api/vector_c_api.html#ogr__api_8h_1acaed6926b75cd33a42b284c10def6e87
+    /// [has_geos]: crate::version::VersionInfo::has_geos
+    pub fn intersects(&self, other: &Self) -> bool {
+        let p = unsafe { gdal_sys::OGR_G_Intersects(self.c_geometry(), other.c_geometry()) };
+        p != 0
+    }
+
+    /// Tests if this geometry [_contains_][DE-9IM] the other geometry;
+    /// `other` lies in `self`, and the interiors intersect.
+    ///
+    /// # Notes
+    /// * Geometry validity is not checked, and invalid geometry will generate unpredictable results.
+    /// Use [`Geometry::is_valid`] if validity might be in question.
+    /// * If GEOS is *not* enabled, this function will always return `false`. Check with [`VersionInfo::has_geos`][has_geos].
+    ///
+    /// See: [`OGR_G_Contains`][OGR_G_Contains]
+    ///
+    /// [DE-9IM]: https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates
+    /// [OGR_G_Contains]: https://gdal.org/api/vector_c_api.html#_CPPv414OGR_G_Contains12OGRGeometryH12OGRGeometryH
+    /// [has_geos]: crate::version::VersionInfo::has_geos
+    pub fn contains(&self, other: &Self) -> bool {
+        let p = unsafe { gdal_sys::OGR_G_Contains(self.c_geometry(), other.c_geometry()) };
+        p != 0
+    }
+
+    /// Tests if this geometry and the other geometry are [_disjoint_][DE-9IM];
+    /// `self` and `other` form a set of disconnected geometries.
+    ///
+    /// # Notes
+    /// * Geometry validity is not checked, and invalid geometry will generate unpredictable results.
+    /// Use [`Geometry::is_valid`] if validity might be in question.
+    /// * If GEOS is *not* enabled, this function will always return `false`. Check with [`VersionInfo::has_geos`][has_geos].
+    ///
+    /// See: [`OGR_G_Disjoint`][OGR_G_Disjoint]
+    ///
+    /// [DE-9IM]: https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates
+    /// [OGR_G_Disjoint]: https://gdal.org/api/vector_c_api.html#_CPPv414OGR_G_Disjoint12OGRGeometryH12OGRGeometryH
+    /// [has_geos]: crate::version::VersionInfo::has_geos
+    pub fn disjoint(&self, other: &Self) -> bool {
+        let p = unsafe { gdal_sys::OGR_G_Disjoint(self.c_geometry(), other.c_geometry()) };
+        p != 0
+    }
+
+    /// Tests if this geometry and the other geometry are [_touching_][DE-9IM];
+    /// `self` and `other` have at least one point in common, but their interiors do not intersect.
+    ///
+    /// # Notes
+    /// * Geometry validity is not checked, and invalid geometry will generate unpredictable results.
+    /// Use [`Geometry::is_valid`] if validity might be in question.
+    /// * If GEOS is *not* enabled, this function will always return `false`. Check with [`VersionInfo::has_geos`][has_geos].
+    ///
+    /// See: [`OGR_G_Touches`][OGR_G_Touches]
+    ///
+    /// [DE-9IM]: https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates
+    /// [OGR_G_Touches]: https://gdal.org/api/ogrgeometry_cpp.html#_CPPv4NK11OGRGeometry7TouchesEPK11OGRGeometry
+    /// [has_geos]: crate::version::VersionInfo::has_geos
+    pub fn touches(&self, other: &Self) -> bool {
+        let p = unsafe { gdal_sys::OGR_G_Touches(self.c_geometry(), other.c_geometry()) };
+        p != 0
+    }
+
+    /// Tests if this geometry and the other geometry are [_crossing_][DE-9IM];
+    /// `self` and `other` have some but not all interior points in common, and the dimension of
+    /// the intersection is less than that of at least one of them.
+    ///
+    /// # Notes
+    /// * Geometry validity is not checked, and invalid geometry will generate unpredictable results.
+    /// Use [`Geometry::is_valid`] if validity might be in question.
+    /// * If GEOS is *not* enabled, this function will always return `false`. Check with [`VersionInfo::has_geos`][has_geos].
+    ///
+    /// See: [`OGR_G_Crosses`][OGR_G_Crosses]
+    ///
+    /// [DE-9IM]: https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates
+    /// [OGR_G_Crosses]: https://gdal.org/api/ogrgeometry_cpp.html#_CPPv4NK11OGRGeometry7CrossesEPK11OGRGeometry
+    /// [has_geos]: crate::version::VersionInfo::has_geos
+    pub fn crosses(&self, other: &Self) -> bool {
+        let p = unsafe { gdal_sys::OGR_G_Crosses(self.c_geometry(), other.c_geometry()) };
+        p != 0
+    }
+
+    /// Tests if this geometry is [_within_][DE-9IM] the other;
+    /// `self` lies fully in the interior of `other`.
+    ///
+    /// # Notes
+    /// * Geometry validity is not checked, and invalid geometry will generate unpredictable results.
+    /// Use [`Geometry::is_valid`] if validity might be in question.
+    /// * If GEOS is *not* enabled, this function will always return `false`. Check with [`VersionInfo::has_geos`][has_geos].
+    ///
+    /// See: [`OGR_G_Within`][OGR_G_Within]
+    ///
+    /// [DE-9IM]: https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates
+    /// [OGR_G_Within]: https://gdal.org/api/ogrgeometry_cpp.html#_CPPv4NK11OGRGeometry6WithinEPK11OGRGeometry
+    /// [has_geos]: crate::version::VersionInfo::has_geos
+    pub fn within(&self, other: &Self) -> bool {
+        let p = unsafe { gdal_sys::OGR_G_Within(self.c_geometry(), other.c_geometry()) };
+        p != 0
+    }
+
+    /// Tests if this geometry and the other [_overlap_][DE-9IM];
+    /// `self` and `other` they have some but not all points in common,
+    /// they have the same dimension,
+    /// and the intersection of the interiors has the same dimension as the geometries.
+    ///
+    /// # Notes
+    /// * Geometry validity is not checked, and invalid geometry will generate unpredictable results.
+    /// Use [`Geometry::is_valid`] if validity might be in question.
+    /// * If GEOS is *not* enabled, this function will always return `false`. Check with [`VersionInfo::has_geos`][has_geos].
+    ///
+    /// See: [`OGR_G_Overlaps`][OGR_G_Overlaps]
+    ///
+    /// [DE-9IM]: https://en.wikipedia.org/wiki/DE-9IM#Spatial_predicates
+    /// [OGR_G_Overlaps]: https://gdal.org/api/ogrgeometry_cpp.html#_CPPv4NK11OGRGeometry8OverlapsEPK11OGRGeometry
+    /// [has_geos]: crate::version::VersionInfo::has_geos
+    pub fn overlaps(&self, other: &Self) -> bool {
+        let p = unsafe { gdal_sys::OGR_G_Overlaps(self.c_geometry(), other.c_geometry()) };
+        p != 0
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    #[test]
+    fn test_intersects() {
+        let poly = Geometry::from_wkt("POLYGON((0 0,5 5,10 0,0 0))").unwrap();
+        let point = Geometry::from_wkt("POINT(10 0)").unwrap();
+        assert!(poly.intersects(&point));
+    }
+
+    #[test]
+    fn test_contains() {
+        let poly = Geometry::from_wkt("POLYGON((0 0,5 5,10 0,0 0))").unwrap();
+        let point = Geometry::from_wkt("POINT(10 0)").unwrap();
+        assert!(!poly.contains(&point));
+        let point = Geometry::from_wkt("POINT(0.1 0.01)").unwrap();
+        assert!(poly.contains(&point));
+    }
+
+    #[test]
+    fn test_disjoint() {
+        let poly = Geometry::from_wkt("POLYGON((0 0,5 5,10 0,0 0))").unwrap();
+        let line = Geometry::from_wkt("LINESTRING(-1 -1, -2 -2)").unwrap();
+        assert!(poly.disjoint(&line));
+    }
+
+    #[test]
+    fn test_touches() {
+        let line1 = Geometry::from_wkt("LINESTRING(0 0, 10 10)").unwrap();
+        let line2 = Geometry::from_wkt("LINESTRING(0 0, 0 10)").unwrap();
+        assert!(line1.touches(&line2));
+    }
+
+    #[test]
+    fn test_crosses() {
+        let line1 = Geometry::from_wkt("LINESTRING(0 0, 10 10)").unwrap();
+        let line2 = Geometry::from_wkt("LINESTRING(10 0, 0 10)").unwrap();
+        assert!(line1.crosses(&line2));
+    }
+
+    #[test]
+    fn test_within() {
+        let poly1 = Geometry::from_wkt("POLYGON((0 0, 10 10, 10 0, 0 0))").unwrap();
+        let poly2 = Geometry::from_wkt("POLYGON((-90 -90, -90 90, 190 -90, -90 -90))").unwrap();
+        assert!(poly1.within(&poly2));
+        assert!(!poly2.within(&poly1));
+    }
+
+    #[test]
+    fn test_overlaps() {
+        let poly1 = Geometry::from_wkt("POLYGON((0 0, 10 10, 10 0, 0 0))").unwrap();
+        let poly2 = Geometry::from_wkt("POLYGON((0 -5,10 5,10 -5,0 -5))").unwrap();
+        assert!(poly1.overlaps(&poly2));
+    }
+}