ST_PointOnSurface, ST_ConcaveHull

rust/geoarrow/src/algorithm/geo/concave_hull.rs

@@ -0,0 +1,119 @@
+use crate::algorithm::broadcasting::BroadcastablePrimitive;
+use crate::array::*;
+use crate::datatypes::{Dimension, NativeType};
+use crate::error::{GeoArrowError, Result};
+use crate::trait_::ArrayAccessor;
+use crate::NativeArray;
+use arrow::datatypes::Float64Type;
+use geo::algorithm::ConcaveHull as _;
+use geo::Polygon;
+
+/// Returns a polygon which covers a geometry. Unlike convex hulls, which also cover
+/// their geometry, a concave hull does so while trying to further minimize its area by
+/// constructing edges such that the exterior of the polygon incorporates points that would
+/// be interior points in a convex hull.
+///
+/// This implementation is inspired by <https://github.com/mapbox/concaveman>
+/// and also uses ideas from the following paper:
+/// www.iis.sinica.edu.tw/page/jise/2012/201205_10.pdf
+pub trait ConcaveHull {
+    type Output;
+
+    fn concave_hull(&self, concavity: &BroadcastablePrimitive<Float64Type>) -> Self::Output;
+}
+
+/// Implementation that iterates over geo objects
+macro_rules! iter_geo_impl {
+    ($type:ty) => {
+        impl ConcaveHull for $type {
+            type Output = PolygonArray;
+
+            fn concave_hull(
+                &self,
+                concavity: &BroadcastablePrimitive<Float64Type>,
+            ) -> Self::Output {
+                let output_geoms: Vec<Option<Polygon>> = self
+                    .iter_geo()
+                    .zip(concavity)
+                    .map(|(maybe_g, concavity)| {
+                        if let (Some(geom), Some(concavity)) = (maybe_g, concavity) {
+                            Some(geom.concave_hull(concavity))
+                        } else {
+                            None
+                        }
+                    })
+                    .collect();
+
+                PolygonBuilder::from_nullable_polygons(
+                    output_geoms.as_slice(),
+                    Dimension::XY,
+                    self.coord_type(),
+                    self.metadata().clone(),
+                )
+                .finish()
+            }
+        }
+    };
+}
+
+iter_geo_impl!(LineStringArray);
+iter_geo_impl!(PolygonArray);
+iter_geo_impl!(MultiPointArray);
+iter_geo_impl!(MultiLineStringArray);
+iter_geo_impl!(MultiPolygonArray);
+
+impl ConcaveHull for GeometryArray {
+    type Output = Result<PolygonArray>;
+
+    fn concave_hull(&self, concavity: &BroadcastablePrimitive<Float64Type>) -> Self::Output {
+        let output_geoms: Vec<Option<Polygon>> = self
+            .iter_geo()
+            .zip(concavity)
+            .map(|(maybe_g, concavity)| {
+                if let (Some(geom), Some(concavity)) = (maybe_g, concavity) {
+                    let out = match geom {
+                        geo::Geometry::LineString(g) => g.concave_hull(concavity),
+                        geo::Geometry::Polygon(g) => g.concave_hull(concavity),
+                        geo::Geometry::MultiLineString(g) => g.concave_hull(concavity),
+                        geo::Geometry::MultiPolygon(g) => g.concave_hull(concavity),
+                        _ => {
+                            return Err(GeoArrowError::IncorrectType(
+                                "incorrect type in concave_hull".into(),
+                            ))
+                        }
+                    };
+                    Ok(Some(out))
+                } else {
+                    Ok(None)
+                }
+            })
+            .collect::<Result<_>>()?;
+
+        Ok(PolygonBuilder::from_nullable_polygons(
+            output_geoms.as_slice(),
+            Dimension::XY,
+            self.coord_type(),
+            self.metadata().clone(),
+        )
+        .finish())
+    }
+}
+
+impl ConcaveHull for &dyn NativeArray {
+    type Output = Result<PolygonArray>;
+
+    fn concave_hull(&self, concavity: &BroadcastablePrimitive<Float64Type>) -> Self::Output {
+        use NativeType::*;
+
+        let result = match self.data_type() {
+            LineString(_, _) => self.as_line_string().concave_hull(concavity),
+            Polygon(_, _) => self.as_polygon().concave_hull(concavity),
+            MultiPoint(_, _) => self.as_multi_point().concave_hull(concavity),
+            MultiLineString(_, _) => self.as_multi_line_string().concave_hull(concavity),
+            MultiPolygon(_, _) => self.as_multi_polygon().concave_hull(concavity),
+            Geometry(_) => self.as_geometry().concave_hull(concavity)?,
+            _ => return Err(GeoArrowError::IncorrectType("".into())),
+        };
+        Ok(result)
+    }
+}

rust/geoarrow/src/algorithm/geo/interior_point.rs

@@ -0,0 +1,107 @@
+use crate::array::*;
+use crate::datatypes::{Dimension, NativeType};
+use crate::error::Result;
+use crate::trait_::ArrayAccessor;
+use crate::NativeArray;
+use geo::algorithm::interior_point::InteriorPoint as _;
+
+/// Calculation of interior points.
+///
+/// An interior point is a point that's guaranteed to intersect a given geometry, and will be
+/// strictly on the interior of the geometry if possible, or on the edge if the geometry has zero
+/// area. A best effort will additionally be made to locate the point reasonably centrally.
+///
+/// For polygons, this point is located by drawing a line that approximately subdivides the
+/// bounding box around the polygon in half, intersecting it with the polygon, then calculating
+/// the midpoint of the longest line produced by the intersection. For lines, the non-endpoint
+/// vertex closest to the line's centroid is returned if the line has interior points, or an
+/// endpoint is returned otherwise.
+///
+/// For multi-geometries or collections, the interior points of the constituent components are
+/// calculated, and one of those is returned (for MultiPolygons, it's the point that's the midpoint
+/// of the longest intersection of the intersection lines of any of the constituent polygons, as
+/// described above; for all others, the interior point closest to the collection's centroid is
+/// used).
+///
+pub trait InteriorPoint {
+    type Output;
+
+    fn interior_point(&self) -> Self::Output;
+}
+
+impl InteriorPoint for PointArray {
+    type Output = PointArray;
+
+    fn interior_point(&self) -> Self::Output {
+        self.clone()
+    }
+}
+
+impl InteriorPoint for RectArray {
+    type Output = PointArray;
+
+    fn interior_point(&self) -> Self::Output {
+        let mut output_array = PointBuilder::with_capacity_and_options(
+            Dimension::XY,
+            self.len(),
+            self.coord_type(),
+            self.metadata().clone(),
+        );
+        self.iter_geo().for_each(|maybe_g| {
+            output_array.push_point(maybe_g.map(|g| g.interior_point()).as_ref())
+        });
+        output_array.into()
+    }
+}
+
+/// Implementation that iterates over geo objects
+macro_rules! iter_geo_impl {
+    ($type:ty) => {
+        impl InteriorPoint for $type {
+            type Output = PointArray;
+
+            fn interior_point(&self) -> Self::Output {
+                let mut output_array = PointBuilder::with_capacity_and_options(
+                    Dimension::XY,
+                    self.len(),
+                    self.coord_type(),
+                    self.metadata().clone(),
+                );
+                self.iter_geo().for_each(|maybe_g| {
+                    output_array.push_point(maybe_g.and_then(|g| g.interior_point()).as_ref())
+                });
+                output_array.into()
+            }
+        }
+    };
+}
+
+iter_geo_impl!(LineStringArray);
+iter_geo_impl!(PolygonArray);
+iter_geo_impl!(MultiPointArray);
+iter_geo_impl!(MultiLineStringArray);
+iter_geo_impl!(MultiPolygonArray);
+iter_geo_impl!(MixedGeometryArray);
+iter_geo_impl!(GeometryCollectionArray);
+iter_geo_impl!(GeometryArray);
+
+impl InteriorPoint for &dyn NativeArray {
+    type Output = Result<PointArray>;
+
+    fn interior_point(&self) -> Self::Output {
+        use NativeType::*;
+
+        let result = match self.data_type() {
+            Point(_, _) => self.as_point().interior_point(),
+            LineString(_, _) => self.as_line_string().interior_point(),
+            Polygon(_, _) => self.as_polygon().interior_point(),
+            MultiPoint(_, _) => self.as_multi_point().interior_point(),
+            MultiLineString(_, _) => self.as_multi_line_string().interior_point(),
+            MultiPolygon(_, _) => self.as_multi_polygon().interior_point(),
+            GeometryCollection(_, _) => self.as_geometry_collection().interior_point(),
+            Rect(_) => self.as_rect().interior_point(),
+            Geometry(_) => self.as_geometry().interior_point(),
+        };
+        Ok(result)
+    }
+}

rust/geoarrow/src/algorithm/geo/mod.rs

@@ -34,6 +34,10 @@ pub use chaikin_smoothing::ChaikinSmoothing;
 mod chamberlain_duquette_area;
 pub use chamberlain_duquette_area::ChamberlainDuquetteArea;
 
+/// Calculate the concave hull of geometries.
+mod concave_hull;
+pub use concave_hull::ConcaveHull;
+
 /// Determine whether `Geometry` `A` completely encloses `Geometry` `B`.
 mod contains;
 pub use contains::Contains;
@@ -74,6 +78,10 @@ pub use geodesic_length::GeodesicLength;
 mod haversine_length;
 pub use haversine_length::HaversineLength;
 
+/// Calculation of interior points.
+mod interior_point;
+pub use interior_point::InteriorPoint;
+
 /// Determine whether `Geometry` `A` intersects `Geometry` `B`.
 mod intersects;
 pub use intersects::Intersects;

rust/geodatafusion/README.md

@@ -266,7 +266,7 @@ Spatial extensions for [Apache DataFusion](https://datafusion.apache.org/), an e
 | ST_BuildArea                |             | Creates a polygonal geometry formed by the linework of a geometry.                                |
 | ST_Centroid                 | ✅          | Returns the geometric center of a geometry.                                                       |
 | ST_ChaikinSmoothing         |             | Returns a smoothed version of a geometry, using the Chaikin algorithm                             |
-| ST_ConcaveHull              |             | Computes a possibly concave geometry that contains all input geometry vertices                    |
+| ST_ConcaveHull              | ✅          | Computes a possibly concave geometry that contains all input geometry vertices                    |
 | ST_ConvexHull               | ✅          | Computes the convex hull of a geometry.                                                           |
 | ST_DelaunayTriangles        |             | Returns the Delaunay triangulation of the vertices of a geometry.                                 |
 | ST_FilterByM                |             | Removes vertices based on their M value                                                           |
@@ -279,7 +279,7 @@ Spatial extensions for [Apache DataFusion](https://datafusion.apache.org/), an e
 | ST_MinimumBoundingRadius    |             | Returns the center point and radius of the smallest circle that contains a geometry.              |
 | ST_OrientedEnvelope         |             | Returns a minimum-area rectangle containing a geometry.                                           |
 | ST_OffsetCurve              |             | Returns an offset line at a given distance and side from an input line.                           |
-| ST_PointOnSurface           |             | Computes a point guaranteed to lie in a polygon, or on a geometry.                                |
+| ST_PointOnSurface           | ✅          | Computes a point guaranteed to lie in a polygon, or on a geometry.                                |
 | ST_Polygonize               |             | Computes a collection of polygons formed from the linework of a set of geometries.                |
 | ST_ReducePrecision          |             | Returns a valid geometry with points rounded to a grid tolerance.                                 |
 | ST_SharedPaths              |             | Returns a collection containing paths shared by the two input linestrings/multilinestrings.       |

rust/geodatafusion/src/udf/native/processing/concave_hull.rs

@@ -0,0 +1,98 @@
+use std::any::Any;
+use std::sync::OnceLock;
+
+use arrow::array::AsArray;
+use arrow::datatypes::Float64Type;
+use arrow_schema::DataType;
+use datafusion::logical_expr::scalar_doc_sections::DOC_SECTION_OTHER;
+use datafusion::logical_expr::{
+    ColumnarValue, Documentation, ScalarUDFImpl, Signature, Volatility,
+};
+use geoarrow::algorithm::broadcasting::BroadcastablePrimitive;
+use geoarrow::algorithm::geo::ConcaveHull as _;
+use geoarrow::array::GeometryArray;
+use geoarrow::ArrayBase;
+
+use crate::data_types::{parse_to_native_array, GEOMETRY_TYPE, POINT2D_TYPE};
+use crate::error::GeoDataFusionResult;
+
+#[derive(Debug)]
+pub(super) struct ConcaveHull {
+    signature: Signature,
+}
+
+impl ConcaveHull {
+    pub fn new() -> Self {
+        Self {
+            signature: Signature::exact(
+                vec![GEOMETRY_TYPE.into(), DataType::Float64],
+                Volatility::Immutable,
+            ),
+        }
+    }
+}
+
+static DOCUMENTATION: OnceLock<Documentation> = OnceLock::new();
+
+impl ScalarUDFImpl for ConcaveHull {
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+
+    fn name(&self) -> &str {
+        "st_concavehull"
+    }
+
+    fn signature(&self) -> &Signature {
+        &self.signature
+    }
+
+    fn return_type(&self, _arg_types: &[DataType]) -> datafusion::error::Result<DataType> {
+        Ok(POINT2D_TYPE.into())
+    }
+
+    fn invoke(&self, args: &[ColumnarValue]) -> datafusion::error::Result<ColumnarValue> {
+        Ok(concave_hull_impl(args)?)
+    }
+
+    fn documentation(&self) -> Option<&Documentation> {
+        Some(DOCUMENTATION.get_or_init(|| {
+            Documentation::builder(
+                DOC_SECTION_OTHER,
+                "A concave hull is a (usually) concave geometry which contains the input, and whose vertices are a subset of the input vertices. In the general case the concave hull is a Polygon. The concave hull of two or more collinear points is a two-point LineString. The concave hull of one or more identical points is a Point. The polygon will not contain holes unless the optional param_allow_holes argument is specified as true.
+
+One can think of a concave hull as \"shrink-wrapping\" a set of points. This is different to the convex hull, which is more like wrapping a rubber band around the points. A concave hull generally has a smaller area and represents a more natural boundary for the input points.
+
+The param_pctconvex controls the concaveness of the computed hull. A value of 1 produces the convex hull. Values between 1 and 0 produce hulls of increasing concaveness. A value of 0 produces a hull with maximum concaveness (but still a single polygon). Choosing a suitable value depends on the nature of the input data, but often values between 0.3 and 0.1 produce reasonable results.",
+                "ST_ConcaveHull(geometry)",
+            )
+            .with_argument("g1", "geometry")
+            .with_argument("param_pctconvex", "float")
+            .build()
+        }))
+    }
+}
+
+fn concave_hull_impl(args: &[ColumnarValue]) -> GeoDataFusionResult<ColumnarValue> {
+    let array = ColumnarValue::values_to_arrays(&args[..1])?
+        .into_iter()
+        .next()
+        .unwrap();
+    let native_array = parse_to_native_array(array)?;
+    let output = match &args[1] {
+        ColumnarValue::Scalar(concavity) => {
+            let concavity = concavity.to_scalar()?.into_inner();
+            let concavity = concavity.as_primitive::<Float64Type>().value(0);
+            native_array.as_ref().concave_hull(&concavity.into())?
+        }
+        ColumnarValue::Array(concavity) => {
+            native_array
+                .as_ref()
+                .concave_hull(&BroadcastablePrimitive::Array(
+                    concavity.as_primitive().clone(),
+                ))?
+        }
+    };
+
+    Ok(GeometryArray::from(output).to_array_ref().into())
+}

rust/geodatafusion/src/udf/native/processing/mod.rs

@@ -1,6 +1,8 @@
 mod centroid;
 mod chaikin_smoothing;
+mod concave_hull;
 mod convex_hull;
+mod point_on_surface;
 mod simplify;
 mod simplify_preserve_topology;
 mod simplify_vw;
@@ -10,7 +12,9 @@ use datafusion::prelude::SessionContext;
 /// Register all provided [geo] functions for processing geometries
 pub fn register_udfs(ctx: &SessionContext) {
     ctx.register_udf(centroid::Centroid::new().into());
+    ctx.register_udf(concave_hull::ConcaveHull::new().into());
     ctx.register_udf(convex_hull::ConvexHull::new().into());
+    ctx.register_udf(point_on_surface::PointOnSurface::new().into());
     ctx.register_udf(simplify_preserve_topology::SimplifyPreserveTopology::new().into());
     ctx.register_udf(simplify_vw::SimplifyVw::new().into());
     ctx.register_udf(simplify::Simplify::new().into());