Repository: geopandas/geopandas
Pull request #2909 (2023-05-29 16:46:45)

ENH: add hausdorff distance method
implements hausdorff distance, from #2010 

kaushiksrini: it looks like a test is timing out and failing despite only changelog change. is there a way to rerun the checks?

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 876a9ea3..d631f664 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 New methods:
 
+- Added ``hausdorff_distance()`` method from shapely to GeoSeries/GeoDataframe (#2909).
 - Added ``delaunay_triangles`` method from shapely to GeoSeries/GeoDataframe (#2907)
 - Added ``concave_hull`` method from shapely to GeoSeries/GeoDataframe. (#2903)
 - Added ``offset_curve`` method from shapely to GeoSeries/GeoDataframe. (#2902)
@@ -12,6 +13,7 @@ New features and improvements:
 
 Bug fixes:
 
+- ``assert_geodataframe_equal`` now handles GeoDataFrames with no active geometry (#2498)
 
 ## Version 0.13.2 (Jun 6, 2023)
 
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index f5844795..134a9b86 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -24,6 +24,7 @@ General methods and attributes
    GeoSeries.geom_type
    GeoSeries.offset_curve
    GeoSeries.distance
+   GeoSeries.hausdorff_distance
    GeoSeries.representative_point
    GeoSeries.exterior
    GeoSeries.interiors
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index 92d3666e..be4c1d54 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -930,6 +930,20 @@ def distance(data, other):
         return _binary_op_float("distance", data, other)
 
 
+def hausdorff_distance(data, other, densify=None, **kwargs):
+    if compat.USE_SHAPELY_20:
+        return shapely.hausdorff_distance(data, other, densify=densify, **kwargs)
+    elif compat.USE_PYGEOS:
+        return _binary_method(
+            "hausdorff_distance", data, other, densify=densify, **kwargs
+        )
+    else:
+        raise NotImplementedError(
+            f"shapely >= 2.0 or PyGEOS is required, "
+            f"version {shapely.__version__} is installed"
+        )
+
+
 def buffer(data, distance, resolution=16, **kwargs):
     if compat.USE_SHAPELY_20:
         if compat.SHAPELY_G_20a1:
diff --git a/geopandas/array.py b/geopandas/array.py
index 37dfcfd6..b9d2ca1f 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -657,6 +657,10 @@ class GeometryArray(ExtensionArray):
         self.check_geographic_crs(stacklevel=6)
         return self._binary_method("distance", self, other)
 
+    def hausdorff_distance(self, other, **kwargs):
+        self.check_geographic_crs(stacklevel=6)
+        return self._binary_method("hausdorff_distance", self, other, **kwargs)
+
     def buffer(self, distance, resolution=16, **kwargs):
         if not (isinstance(distance, (int, float)) and distance == 0):
             self.check_geographic_crs(stacklevel=5)
diff --git a/geopandas/base.py b/geopandas/base.py
index 09d641fe..d15cecaf 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -2384,6 +2384,121 @@ GeometryCollection
         """
         return _binary_op("distance", self, other, align)
 
+    def hausdorff_distance(self, other, align=True, densify=None):
+        """Returns a ``Series`` containing the Hausdorff distance to aligned `other`.
+
+        The Hausdorff distance is the largest distance consisting of any point in `self`
+        with the nearest point in `other`.
+
+        The operation works on a 1-to-1 row-wise manner:
+
+        .. image:: ../../../_static/binary_op-01.svg
+           :align: center
+
+        Parameters
+        ----------
+        other : GeoSeries or geometric object
+            The Geoseries (elementwise) or geometric object to find the
+            distance to.
+        align : bool (default True)
+            If True, automatically aligns GeoSeries based on their indices.
+            If False, the order of elements is preserved.
+        densify : float (default None)
+            A value between 0 and 1, that splits each subsegment of a line string
+            into equal length segments, making the approximation less coarse.
+            A densify value of 0.5 will add a point halfway between each pair of
+            points. A densify value of 0.25 will add a point a quarter of the way
+            between each pair of points.
+
+
+        Returns
+        -------
+        Series (float)
+
+        Examples
+        --------
+        >>> from shapely.geometry import Polygon, LineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (1, 0), (1, 1)]),
+        ...         Polygon([(0, 0), (-1, 0), (-1, 1)]),
+        ...         LineString([(1, 1), (0, 0)]),
+        ...         Point(0, 0),
+        ...     ],
+        ... )
+        >>> s2 = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0.5, 0.5), (1.5, 0.5), (1.5, 1.5), (0.5, 1.5)]),
+        ...         Point(3, 1),
+        ...         LineString([(1, 0), (2, 0)]),
+        ...         Point(0, 1),
+        ...     ],
+        ...     index=range(1, 5),
+        ... )
+
+        >>> s
+        0    POLYGON ((0.00000 0.00000, 1.00000 0.00000, 1....
+        1    POLYGON ((0.00000 0.00000, -1.00000 0.00000, -...
+        2        LINESTRING (1.00000 1.00000, 0.00000 0.00000)
+        3                              POINT (0.00000 0.00000)
+        dtype: geometry
+
+        >>> s2
+        1    POLYGON ((0.50000 0.50000, 1.50000 0.50000, 1....
+        2                              POINT (3.00000 1.00000)
+        3        LINESTRING (1.00000 0.00000, 2.00000 0.00000)
+        4                              POINT (0.00000 1.00000)
+        dtype: geometry
+
+        We can check the hausdorff distance of each geometry of GeoSeries
+        to a single geometry:
+
+        >>> point = Point(-1, 0)
+        >>> s.hausdorff_distance(point)
+        0    2.236068
+        1    1.000000
+        2    2.236068
+        3    1.000000
+        dtype: float64
+
+        We can also check two GeoSeries against each other, row by row.
+        The GeoSeries above have different indices. We can either align both GeoSeries
+        based on index values and use elements with the same index using
+        ``align=True`` or ignore index and use elements based on their matching
+        order using ``align=False``:
+
+        .. image:: ../../../_static/binary_op-02.svg
+
+        >>> s.hausdorff_distance(s2, align=True)
+        0         NaN
+        1    2.121320
+        2    3.162278
+        3    2.000000
+        4         NaN
+        dtype: float64
+
+        >>> s.hausdorff_distance(s2, align=False)
+        0    0.707107
+        1    4.123106
+        2    1.414214
+        3    1.000000
+        dtype: float64
+
+        We can also set a densify value, which is a float between 0 and 1 and
+        signifies the fraction of the distance between each pair of points that will
+        be used as the distance between the points when densifying.
+
+        >>> l1 = geopandas.GeoSeries([LineString([(130, 0), (0, 0), (0, 150)])])
+        >>> l2 = geopandas.GeoSeries([LineString([(10, 10), (10, 150), (130, 10)])])
+        >>> l1.hausdorff_distance(l2)
+        0    14.142136
+        dtype: float64
+        >>> l1.hausdorff_distance(l2, densify=0.25)
+        0    70.0
+        dtype: float64
+        """
+        return _binary_op("hausdorff_distance", self, other, align, densify=densify)
+
     #
     # Binary operations that return a GeoSeries
     #


Test patch:
diff --git a/geopandas/testing.py b/geopandas/testing.py
index ae89cd0b..54b68686 100644
--- a/geopandas/testing.py
+++ b/geopandas/testing.py
@@ -288,8 +288,19 @@ def assert_geodataframe_equal(
         assert isinstance(left, type(right))
 
         if check_crs:
+            # allow if neither left and right has an active geometry column
+            if (
+                left._geometry_column_name is None
+                and right._geometry_column_name is None
+            ):
+                pass
+            elif (
+                left._geometry_column_name not in left.columns
+                and right._geometry_column_name not in right.columns
+            ):
+                pass
             # no crs can be either None or {}
-            if not left.crs and not right.crs:
+            elif not left.crs and not right.crs:
                 pass
             else:
                 assert left.crs == right.crs
@@ -328,9 +339,12 @@ def assert_geodataframe_equal(
                 check_crs=check_crs,
             )
 
+    # ensure the active geometry column is the same
+    assert left._geometry_column_name == right._geometry_column_name
+
     # drop geometries and check remaining columns
-    left2 = left.drop([left._geometry_column_name], axis=1)
-    right2 = right.drop([right._geometry_column_name], axis=1)
+    left2 = left.select_dtypes(exclude="geometry")
+    right2 = right.select_dtypes(exclude="geometry")
     assert_frame_equal(
         left2,
         right2,
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 7e1cad22..83879ff5 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -172,6 +172,11 @@ class TestGeomMethods:
         )
         self.squares = GeoSeries([self.sq for _ in range(3)])
 
+        self.l5 = LineString([(100, 0), (0, 0), (0, 100)])
+        self.l6 = LineString([(5, 5), (5, 100), (100, 5)])
+        self.g12 = GeoSeries([self.l5])
+        self.g13 = GeoSeries([self.l6])
+
     def _test_unary_real(self, op, expected, a):
         """Tests for 'area', 'length', 'is_valid', etc."""
         fcmp = assert_series_equal
@@ -566,6 +571,53 @@ class TestGeomMethods:
         with pytest.warns(UserWarning, match="Geometry is in a geographic CRS"):
             self.g4.distance(self.p0)
 
+    @pytest.mark.skipif(
+        not (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="requires hausdorff_distance in shapely 2.0+",
+    )
+    def test_hausdorff_distance(self):
+        # closest point is (0, 0) in self.p1
+        expected = Series(
+            np.array([np.sqrt(5**2 + 5**2), np.nan]), self.na_none.index
+        )
+        assert_array_dtype_equal(expected, self.na_none.hausdorff_distance(self.p0))
+
+        expected = Series(
+            np.array([np.sqrt(5**2 + 5**2), np.nan]), self.na_none.index
+        )
+        assert_array_dtype_equal(expected, self.na_none.hausdorff_distance(self.p0))
+
+        expected = Series(np.array([np.nan, 0, 0, 0, 0, 0, np.nan, np.nan]), range(8))
+        with pytest.warns(UserWarning, match="The indices .+ different"):
+            assert_array_dtype_equal(
+                expected, self.g0.hausdorff_distance(self.g9, align=True)
+            )
+
+        val_1 = self.g0.iloc[0].hausdorff_distance(self.g9.iloc[0])
+        val_2 = self.g0.iloc[2].hausdorff_distance(self.g9.iloc[2])
+        val_3 = self.g0.iloc[4].hausdorff_distance(self.g9.iloc[4])
+        expected = Series(
+            np.array([val_1, val_1, val_2, val_2, val_3, np.nan, np.nan]), self.g0.index
+        )
+        assert_array_dtype_equal(
+            expected, self.g0.hausdorff_distance(self.g9, align=False)
+        )
+
+        expected = Series(np.array([52.5]), self.g12.index)
+        assert_array_dtype_equal(
+            expected, self.g12.hausdorff_distance(self.g13, densify=0.25)
+        )
+
+    @pytest.mark.skipif(
+        (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="hausdorff_distance not implemented for shapely<2",
+    )
+    def test_hausdorff_distance_not(self):
+        with pytest.raises(
+            NotImplementedError, match="shapely >= 2.0 or PyGEOS is required"
+        ):
+            self.g0.hausdorff_distance(self.g9)
+
     def test_intersects(self):
         expected = [True, True, True, True, True, False, False]
         assert_array_dtype_equal(expected, self.g0.intersects(self.t1))
diff --git a/geopandas/tests/test_pandas_methods.py b/geopandas/tests/test_pandas_methods.py
index 99da2185..c10a4183 100644
--- a/geopandas/tests/test_pandas_methods.py
+++ b/geopandas/tests/test_pandas_methods.py
@@ -776,9 +776,7 @@ def test_pivot(df):
     # pivot failing due to creating a MultiIndex
     result = df.pivot(columns="value1")
     expected = GeoDataFrame(pd.DataFrame(df).pivot(columns="value1"))
-    # TODO assert_geodataframe_equal crashes
-    assert isinstance(result, GeoDataFrame)
-    assert_frame_equal(result, expected)
+    assert_geodataframe_equal(result, expected)
 
 
 def test_preserve_attrs(df):
diff --git a/geopandas/tests/test_testing.py b/geopandas/tests/test_testing.py
index 1783dd23..b119f72b 100644
--- a/geopandas/tests/test_testing.py
+++ b/geopandas/tests/test_testing.py
@@ -3,7 +3,8 @@ import warnings
 import numpy as np
 
 from shapely.geometry import Point, Polygon
-from pandas import Series
+import pandas as pd
+from pandas import DataFrame, Series
 
 from geopandas import GeoDataFrame, GeoSeries
 from geopandas.array import from_shapely
@@ -139,3 +140,44 @@ def test_almost_equal_but_not_equal():
     assert_geoseries_equal(s_origin, s_almost_origin, check_less_precise=True)
     with pytest.raises(AssertionError):
         assert_geoseries_equal(s_origin, s_almost_origin)
+
+
+def test_geodataframe_no_active_geometry_column():
+    def create_dataframe():
+        gdf = GeoDataFrame({"value": [1, 2], "geometry": [Point(1, 1), Point(2, 2)]})
+        gdf["geom2"] = GeoSeries([Point(3, 3), Point(4, 4)])
+        return gdf
+
+    # no active geometry column (None)
+    df1 = create_dataframe()
+    df1._geometry_column_name = None
+    df2 = create_dataframe()
+    df2._geometry_column_name = None
+    assert_geodataframe_equal(df1, df2)
+
+    # active geometry column ("geometry") not present
+    df1 = create_dataframe()[["value", "geom2"]]
+    df2 = create_dataframe()[["value", "geom2"]]
+    assert_geodataframe_equal(df1, df2)
+
+    df1 = GeoDataFrame(create_dataframe()[["value"]])
+    df2 = GeoDataFrame(create_dataframe()[["value"]])
+    assert_geodataframe_equal(df1, df2)
+
+
+def test_geodataframe_multiindex():
+    def create_dataframe():
+        gdf = DataFrame([[Point(0, 0), Point(1, 1)], [Point(2, 2), Point(3, 3)]])
+        gdf = GeoDataFrame(gdf.astype("geometry"))
+        gdf.columns = pd.MultiIndex.from_product([["geometry"], [0, 1]])
+        return gdf
+
+    df1 = create_dataframe()
+    df2 = create_dataframe()
+    assert_geodataframe_equal(df1, df2)
+
+    df1 = create_dataframe()
+    df1._geometry_column_name = None
+    df2 = create_dataframe()
+    df2._geometry_column_name = None
+    assert_geodataframe_equal(df1, df2)