Add ST_HausdorffDistance

common/src/main/java/org/apache/sedona/common/Functions.java

@@ -36,6 +36,7 @@
 import org.locationtech.jts.operation.valid.IsValidOp;
 import org.locationtech.jts.precision.GeometryPrecisionReducer;
 import org.locationtech.jts.simplify.TopologyPreservingSimplifier;
+import org.locationtech.jts.algorithm.distance.DiscreteHausdorffDistance;
 import org.opengis.referencing.FactoryException;
 import org.opengis.referencing.NoSuchAuthorityCodeException;
 import org.opengis.referencing.crs.CoordinateReferenceSystem;
@@ -965,4 +966,12 @@ public static Geometry boundingDiagonal(Geometry geometry) {
  }
  }
 
+ public static Double hausdorffDistance(Geometry g1, Geometry g2, double densityFrac) throws Exception {
+ return GeomUtils.getHausdorffDistance(g1, g2, densityFrac);
+ }
+
+ public static Double hausdorffDistance(Geometry g1, Geometry g2) throws Exception{
+ return GeomUtils.getHausdorffDistance(g1, g2, -1);
+ }
+
 }

common/src/main/java/org/apache/sedona/common/utils/GeomUtils.java

@@ -26,6 +26,7 @@
 import org.locationtech.jts.io.WKTWriter;
 import org.locationtech.jts.operation.polygonize.Polygonizer;
 import org.locationtech.jts.operation.union.UnaryUnionOp;
+import org.locationtech.jts.algorithm.distance.DiscreteHausdorffDistance;
 
 import java.awt.*;
 import java.nio.ByteOrder;
@@ -461,4 +462,13 @@ public static void translateGeom(Geometry geometry, double deltaX, double deltaY
  geometry.geometryChanged();
  }
  }
+
+ public static Double getHausdorffDistance(Geometry g1, Geometry g2, double densityFrac) throws Exception {
+ if (g1.isEmpty() || g2.isEmpty()) return 0.0;
+ DiscreteHausdorffDistance hausdorffDistanceObj = new DiscreteHausdorffDistance(g1, g2);
+ if (densityFrac != -1) {
+ hausdorffDistanceObj.setDensifyFraction(densityFrac);
+ }
+ return hausdorffDistanceObj.distance();
+ }
 }

common/src/test/java/org/apache/sedona/common/FunctionsTest.java

@@ -918,4 +918,75 @@ public void boundingDiagonalSingleVertex() {
  String actual = Functions.boundingDiagonal(point).toText();
  assertEquals(expected, actual);
  }
+
+ @Test
+ public void hausdorffDistanceDefaultGeom2D() throws Exception {
+ Polygon polygon1 = GEOMETRY_FACTORY.createPolygon(coordArray3d(1, 0, 1, 1, 1, 2, 2, 1, 5, 2, 0, 1, 1, 0, 1));
+ Polygon polygon2 = GEOMETRY_FACTORY.createPolygon(coordArray3d(4, 0, 4, 6, 1, 4, 6, 4, 9, 6, 1, 3, 4, 0, 4));
+ Double expected = 5.0;
+ Double actual = Functions.hausdorffDistance(polygon1, polygon2);
+ assertEquals(expected, actual);
+ }
+
+ @Test
+ public void hausdorffDistanceGeom2D() throws Exception {
+ Point point = GEOMETRY_FACTORY.createPoint(new Coordinate(10, 34));
+ LineString lineString = GEOMETRY_FACTORY.createLineString(coordArray(1, 2, 1, 5, 2, 6, 1, 2));
+ Double expected = 33.24154027718932;
+ Double actual = Functions.hausdorffDistance(point, lineString, 0.33);
+ assertEquals(expected, actual);
+ }
+
+ @Test
+ public void hausdorffDistanceInvalidDensityFrac() throws Exception {
+ Point point = GEOMETRY_FACTORY.createPoint(new Coordinate(10, 34));
+ LineString lineString = GEOMETRY_FACTORY.createLineString(coordArray(1, 2, 1, 5, 2, 6, 1, 2));
+ Exception e = assertThrows(IllegalArgumentException.class, () -> Functions.hausdorffDistance(point, lineString, 3));
+ String expected = "Fraction is not in range (0.0 - 1.0]";
+ String actual = e.getMessage();
+ assertEquals(expected, actual);
+ }
+
+ @Test
+ public void hausdorffDistanceDefaultGeomCollection() throws Exception {
+ Polygon polygon = GEOMETRY_FACTORY.createPolygon(coordArray(1, 2, 2, 1, 2, 0, 4, 1, 1, 2));
+ Geometry point1 = GEOMETRY_FACTORY.createPoint(new Coordinate(1, 0));
+ Geometry point2 = GEOMETRY_FACTORY.createPoint(new Coordinate(40, 10));
+ Geometry point3 = GEOMETRY_FACTORY.createPoint(new Coordinate(-10, -40));
+ GeometryCollection multiPoint = GEOMETRY_FACTORY.createGeometryCollection(new Geometry[] {point1, point2, point3});
+ Double actual = Functions.hausdorffDistance(polygon, multiPoint);
+ Double expected = 41.7612260356422;
+ assertEquals(expected, actual);
+ }
+
+ @Test
+ public void hausdorffDistanceGeomCollection() throws Exception {
+ Polygon polygon = GEOMETRY_FACTORY.createPolygon(coordArray(1, 2, 2, 1, 2, 0, 4, 1, 1, 2));
+ LineString lineString1 = GEOMETRY_FACTORY.createLineString(coordArray(1, 1, 2, 1, 4, 4, 5, 5));
+ LineString lineString2 = GEOMETRY_FACTORY.createLineString(coordArray(10, 10, 11, 11, 12, 12, 14, 14));
+ LineString lineString3 = GEOMETRY_FACTORY.createLineString(coordArray(-11, -20, -11, -21, -15, -19));
+ MultiLineString multiLineString = GEOMETRY_FACTORY.createMultiLineString(new LineString[] {lineString1, lineString2, lineString3});
+ Double actual = Functions.hausdorffDistance(polygon, multiLineString, 0.0000001);
+ Double expected = 25.495097567963924;
+ assertEquals(expected, actual);
+ }
+
+ @Test
+ public void hausdorffDistanceEmptyGeom() throws Exception {
+ Polygon polygon = GEOMETRY_FACTORY.createPolygon(coordArray(1, 2, 2, 1, 2, 0, 4, 1, 1, 2));
+ LineString emptyLineString = GEOMETRY_FACTORY.createLineString();
+ Double expected = 0.0;
+ Double actual = Functions.hausdorffDistance(polygon, emptyLineString, 0.00001);
+ assertEquals(expected, actual);
+ }
+
+ @Test
+ public void hausdorffDistanceDefaultEmptyGeom() throws Exception {
+ Polygon polygon = GEOMETRY_FACTORY.createPolygon(coordArray(1, 2, 2, 1, 2, 0, 4, 1, 1, 2));
+ LineString emptyLineString = GEOMETRY_FACTORY.createLineString();
+ Double expected = 0.0;
+ Double actual = Functions.hausdorffDistance(polygon, emptyLineString);
+ assertEquals(expected, actual);
+ }
+
 }

docs/api/flink/Function.md

@@ -533,6 +533,46 @@ SELECT ST_GeometryN(ST_GeomFromText('MULTIPOINT((1 2), (3 4), (5 6), (8 9))'), 1
 
 Output: `POINT (3 4)`
 
+## ST_HausdorffDistance
+
+Introduction: Returns a discretized (and hence approximate) [Hausdorff distance](https://en.wikipedia.org/wiki/Hausdorff_distance) between the given 2 geometries. 
+Optionally, a densityFraction parameter can be specified, which gives more accurate results by densifying segments before computing hausdorff distance between them.
+Each segment is broken down into equal-length subsegments whose ratio with segment length is closest to the given density fraction.
+
+Hence, the lower the densityFrac value, the more accurate is the computed hausdorff distance, and the more time it takes to compute it.
+
+If any of the geometry is empty, 0.0 is returned.
+
+
+!!!Note
+ Accepted range of densityFrac is (0.0, 1.0], if any other value is provided, ST_HausdorffDistance throws an IllegalArgumentException
+
+
+!!!Note
+ Even though the function accepts 3D geometry, the z ordinate is ignored and the computed hausdorff distance is equivalent to the geometries not having the z ordinate.
+
+Format: `ST_HausdorffDistance(g1: geometry, g2: geometry, densityFrac)`
+
+Since: `v1.5.0`
+
+Example:
+```sql
+SELECT ST_HausdorffDistance(g1, g2, 0.1)
+```
+
+Input: `g1: POINT (0.0 1.0), g2: LINESTRING (0 0, 1 0, 2 0, 3 0, 4 0, 5 0)`
+
+Output: `5.0990195135927845`
+
+```sql
+SELECT ST_HausdorffDistance(ST_GeomFromText(), ST_GeomFromText())
+```
+
+Input: `g1: POLYGON Z((1 0 1, 1 1 2, 2 1 5, 2 0 1, 1 0 1)), g2: POLYGON Z((4 0 4, 6 1 4, 6 4 9, 6 1 3, 4 0 4))`
+
+Output: `5.0`
+
+
 ## ST_InteriorRingN
 
 Introduction: Returns the Nth interior linestring ring of the polygon geometry. Returns NULL if the geometry is not a polygon or the given N is out of range

docs/api/sql/Function.md

@@ -731,6 +731,44 @@ SELECT ST_GeometryType(polygondf.countyshape)
 FROM polygondf
 ```
 
+## ST_HausdorffDistance
+
+Introduction: Returns a discretized (and hence approximate) [Hausdorff distance](https://en.wikipedia.org/wiki/Hausdorff_distance) between the given 2 geometries.
+Optionally, a densityFraction parameter can be specified, which gives more accurate results by densifying segments before computing hausdorff distance between them.
+Each segment is broken down into equal-length subsegments whose ratio with segment length is closest to the given density fraction.
+
+Hence, the lower the densityFrac value, the more accurate is the computed hausdorff distance, and the more time it takes to compute it.
+
+If any of the geometry is empty, 0.0 is returned.
+
+!!!Note
+ Accepted range of densityFrac is (0.0, 1.0], if any other value is provided, ST_HausdorffDistance throws an IllegalArgumentException
+
+
+!!!Note
+ Even though the function accepts 3D geometry, the z ordinate is ignored and the computed hausdorff distance is equivalent to the geometries not having the z ordinate.
+
+Format: `ST_HausdorffDistance(g1: geometry, g2: geometry, densityFrac)`
+
+Since: `v1.5.0`
+
+Example:
+```sql
+SELECT ST_HausdorffDistance(g1, g2, 0.1)
+```
+
+Input: `g1: POINT (0.0 1.0), g2: LINESTRING (0 0, 1 0, 2 0, 3 0, 4 0, 5 0)`
+
+Output: `5.0990195135927845`
+
+```sql
+SELECT ST_HausdorffDistance(ST_GeomFromText(), ST_GeomFromText())
+```
+
+Input: `g1: POLYGON Z((1 0 1, 1 1 2, 2 1 5, 2 0 1, 1 0 1)), g2: POLYGON Z((4 0 4, 6 1 4, 6 4 9, 6 1 3, 4 0 4))`
+
+Output: `5.0`
+
 ## ST_InteriorRingN
 
 Introduction: Returns the Nth interior linestring ring of the polygon geometry. Returns NULL if the geometry is not a polygon or the given N is out of range

flink/src/main/java/org/apache/sedona/flink/Catalog.java

@@ -100,6 +100,7 @@ public static UserDefinedFunction[] getFuncs() {
  new Functions.ST_NRings(),
  new Functions.ST_Translate(),
  new Functions.ST_BoundingDiagonal(),
+ new Functions.ST_HausdorffDistance(),
  };
  }
 

flink/src/main/java/org/apache/sedona/flink/expressions/Functions.java

@@ -624,12 +624,29 @@ public Geometry eval(@DataTypeHint(value = "RAW", bridgedTo = org.locationtech.j
  }
 
  public static class ST_BoundingDiagonal extends ScalarFunction {
-
  @DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class)
  public Geometry eval(@DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class) Object o) {
  Geometry geometry = (Geometry) o;
  return org.apache.sedona.common.Functions.boundingDiagonal(geometry);
  }
+ }
 
+ public static class ST_HausdorffDistance extends ScalarFunction {
+ @DataTypeHint("Double")
+ public Double eval(@DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class) Object g1,
+ @DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class) Object g2,
+ @DataTypeHint("Double") Double densityFrac) throws Exception {
+ Geometry geom1 = (Geometry) g1;
+ Geometry geom2 = (Geometry) g2;
+ return org.apache.sedona.common.Functions.hausdorffDistance(geom1, geom2, densityFrac);
+ }
+
+ @DataTypeHint("Double")
+ public Double eval(@DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class) Object g1,
+ @DataTypeHint(value = "RAW", bridgedTo = org.locationtech.jts.geom.Geometry.class) Object g2) throws Exception {
+ Geometry geom1 = (Geometry) g1;
+ Geometry geom2 = (Geometry) g2;
+ return org.apache.sedona.common.Functions.hausdorffDistance(geom1, geom2);
+ }
  }
 }

flink/src/test/java/org/apache/sedona/flink/FunctionTest.java

@@ -745,4 +745,17 @@ public void testBoundingDiagonal() {
  assertEquals(expected, actual);
  }
 
+ @Test
+ public void testHausdorffDistance() {
+ Table polyTable = tableEnv.sqlQuery("SELECT ST_GeomFromWKT('POINT (0.0 1.0)') AS g1, ST_GeomFromWKT('LINESTRING (0 0, 1 0, 2 0, 3 0, 4 0, 5 0)') AS g2");
+ Table actualTable = polyTable.select(call(Functions.ST_HausdorffDistance.class.getSimpleName(), $("g1"), $("g2"), 0.4));
+ Table actualTableDefault = polyTable.select(call(Functions.ST_HausdorffDistance.class.getSimpleName(), $("g1"), $("g2")));
+ Double expected = 5.0990195135927845;
+ Double expectedDefault = 5.0990195135927845;
+ Double actual = (Double) first(actualTable).getField(0);
+ Double actualDefault = (Double) first(actualTableDefault).getField(0);
+ assertEquals(expected, actual);
+ assertEquals(expectedDefault, actualDefault);
+ }
+
 }

python/sedona/sql/st_functions.py

@@ -1288,3 +1288,16 @@ def ST_BoundingDiagonal(geometry: ColumnOrName) -> Column:
  """
 
  return _call_st_function("ST_BoundingDiagonal", geometry)
+
+@validate_argument_types
+def ST_HausdorffDistance(g1: ColumnOrName, g2: ColumnOrName, densityFrac: Optional[Union[ColumnOrName, float]] = -1) -> Column:
+ """
+ Returns discretized (and hence approximate) hausdorff distance between two given geometries.
+ Optionally, a distance fraction can also be provided which decreases the gap between actual and discretized hausforff distance
+ :param g1:
+ :param g2:
+ :param densityFrac: Optional
+ :return:
+ """
+ args = (g1, g2, densityFrac)
+ return _call_st_function("ST_HausdorffDistance", args)

python/tests/sql/test_dataframe_api.py

@@ -90,6 +90,7 @@
  (stf.ST_GeometricMedian, ("multipoint",), "multipoint_geom", "", "POINT (22.500002656424286 21.250001168173426)"),
  (stf.ST_GeometryN, ("geom", 0), "multipoint", "", "POINT (0 0)"),
  (stf.ST_GeometryType, ("point",), "point_geom", "", "ST_Point"),
+ (stf.ST_HausdorffDistance, ("point", "line",), "point_and_line", "", 5.0990195135927845),
  (stf.ST_InteriorRingN, ("geom", 0), "geom_with_hole", "", "LINESTRING (1 1, 2 2, 2 1, 1 1)"),
  (stf.ST_Intersection, ("a", "b"), "overlapping_polys", "", "POLYGON ((2 0, 1 0, 1 1, 2 1, 2 0))"),
  (stf.ST_IsClosed, ("geom",), "closed_linestring_geom", "", True),
@@ -391,6 +392,8 @@ def base_df(self, request):
  return TestDataFrameAPI.spark.sql("SELECT ST_GeomFromWKT('LINESTRING (0 0, 2 1)') AS line, ST_GeomFromWKT('POLYGON ((1 0, 2 0, 2 2, 1 2, 1 0))') AS poly")
  elif request.param == "square_geom":
  return TestDataFrameAPI.spark.sql("SELECT ST_GeomFromWKT('POLYGON ((1 0, 1 1, 2 1, 2 0, 1 0))') AS geom")
+ elif request.param == "point_and_line":
+ return TestDataFrameAPI.spark.sql("SELECT ST_GeomFromWKT('POINT (0.0 1.0)') AS point, ST_GeomFromWKT('LINESTRING (0 0, 1 0, 2 0, 3 0, 4 0, 5 0)') AS line")
  raise ValueError(f"Invalid base_df name passed: {request.param}")
 
  def _id_test_configuration(val):

python/tests/sql/test_function.py

@@ -1104,3 +1104,19 @@ def test_boundingDiagonal(self):
  "1 0))')) AS geom")
  actual = actual_df.selectExpr("ST_AsText(geom)").take(1)[0][0]
  assert expected == actual
+
+ def test_hausdorffDistance(self):
+ expected = 5.0
+ actual_df = self.spark.sql("SELECT ST_HausdorffDistance(ST_GeomFromText('POLYGON ((1 0 1, 1 1 2, 2 1 5, "
+ "2 0 1, 1 0 1))'), ST_GeomFromText('POLYGON ((4 0 4, 6 1 4, 6 4 9, 6 1 3, "
+ "4 0 4))'), 0.5)")
+ actual_df_default = self.spark.sql("SELECT ST_HausdorffDistance(ST_GeomFromText('POLYGON ((1 0 1, 1 1 2, "
+ "2 1 5, "
+ "2 0 1, 1 0 1))'), ST_GeomFromText('POLYGON ((4 0 4, 6 1 4, 6 4 9, 6 1 3, "
+ "4 0 4))'))")
+ actual = actual_df.take(1)[0][0]
+ actual_default = actual_df_default.take(1)[0][0]
+ assert expected == actual
+ assert expected == actual_default
+
+

sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala

@@ -152,6 +152,7 @@ object Catalog {
  function[ST_NRings](),
  function[ST_Translate](0.0),
  function[ST_BoundingDiagonal](),
+ function[ST_HausdorffDistance](-1),
  // Expression for rasters
  function[RS_NormalizedDifference](),
  function[RS_Mean](),

sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/Functions.scala

@@ -1016,3 +1016,10 @@ case class ST_BoundingDiagonal(inputExpressions: Seq[Expression])
  copy(inputExpressions = newChildren)
  }
 }
+
+case class ST_HausdorffDistance(inputExpressions: Seq[Expression])
+ extends InferredTernaryExpression(Functions.hausdorffDistance) with FoldableExpression {
+ protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = {
+ copy(inputExpressions = newChildren)
+ }
+}

sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/st_functions.scala

@@ -333,4 +333,14 @@ object st_functions extends DataFrameAPI {
  def ST_BoundingDiagonal(geometry: String) =
  wrapExpression[ST_BoundingDiagonal](geometry)
 
+ def ST_HausdorffDistance(g1: Column, g2: Column) = wrapExpression[ST_HausdorffDistance](g1, g2, -1)
+
+ def ST_HausdorffDistance(g1: String, g2: String) = wrapExpression[ST_HausdorffDistance](g1, g2, -1);
+
+ def ST_HausdorffDistance(g1: Column, g2: Column, densityFrac: Column) = wrapExpression[ST_HausdorffDistance](g1, g2, densityFrac);
+
+ def ST_HausdorffDistance(g1: String, g2: String, densityFrac: Double) = wrapExpression[ST_HausdorffDistance](g1, g2, densityFrac);
+
+
+
 }

sql/common/src/test/scala/org/apache/sedona/sql/dataFrameAPITestScala.scala

@@ -1003,5 +1003,17 @@ class dataFrameAPITestScala extends TestBaseScala {
  val actual = wKTWriter.write(df.take(1)(0).get(0).asInstanceOf[Geometry])
  assertEquals(expected, actual)
  }
+
+ it("Passed ST_HausdorffDistance") {
+ val polyDf = sparkSession.sql("SELECT ST_GeomFromWKT('POLYGON ((1 2, 2 1, 2 0, 4 1, 1 2))') AS g1, " +
+ "ST_GeomFromWKT('MULTILINESTRING ((1 1, 2 1, 4 4, 5 5), (10 10, 11 11, 12 12, 14 14), (-11 -20, -11 -21, -15 -19))') AS g2")
+ val df = polyDf.select(ST_HausdorffDistance("g1", "g2", 0.05))
+ val dfDefaultValue = polyDf.select(ST_HausdorffDistance("g1", "g2"))
+ val expected = 25.495097567963924
+ val actual = df.take(1)(0).get(0).asInstanceOf[Double]
+ val actualDefaultValue = dfDefaultValue.take(1)(0).get(0).asInstanceOf[Double]
+ assert(expected == actual)
+ assert(expected == actualDefaultValue)
+ }
  }
 }

sql/common/src/test/scala/org/apache/sedona/sql/functionTestScala.scala

@@ -1987,4 +1987,25 @@ class functionTestScala extends TestBaseScala with Matchers with GeometrySample
  assertEquals(expected, actual)
  }
  }
+
+ it ("should pass ST_HausdorffDistance") {
+ val geomTestCases = Map (
+ ("'LINESTRING (1 2, 1 5, 2 6, 1 2)'", "'POINT (10 34)'", 0.34) -> (33.24154027718932, 33.24154027718932),
+ ("'LINESTRING (1 0, 1 1, 2 1, 2 0, 1 0)'", "'POINT EMPTY'", 0.23) -> (0.0, 0.0),
+ ("'POLYGON ((1 2, 2 1, 2 0, 4 1, 1 2))'", "'MULTIPOINT ((1 0), (40 10), (-10 -40))'", 0.0001) -> (41.7612260356422, 41.7612260356422)
+ )
+ for (((geom), expectedResult) <- geomTestCases) {
+ val geom1 = geom._1
+ val geom2 = geom._2
+ val densityFrac = geom._3
+ val df = sparkSession.sql(s"SELECT ST_HausdorffDistance(ST_GeomFromWKT($geom1), ST_GeomFromWKT($geom2), $densityFrac) AS dist")
+ val dfDefaultValue = sparkSession.sql(s"SELECT ST_HausdorffDistance(ST_GeomFromWKT($geom1), ST_GeomFromWKT($geom2)) as dist")
+ val actual = df.take(1)(0).get(0).asInstanceOf[Double]
+ val actualDefaultValue = dfDefaultValue.take(1)(0).get(0).asInstanceOf[Double]
+ val expected = expectedResult._1
+ val expectedDefaultValue = expectedResult._2
+ assert(expected == actual)
+ assert(expectedDefaultValue == actualDefaultValue)
+ }
+ }
 }