Repository: geopandas/geopandas
Pull request #2915 (2023-06-02 14:37:39)

ENH: add extract_unique_points as a GeoSeries property
Implements the shapely extract_unique_points method as part of https://github.com/geopandas/geopandas/issues/2010


Patch:
diff --git a/.gitignore b/.gitignore
index be513b86..10c62aa9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -70,3 +70,5 @@ geopandas/datasets/ne_110m_admin_0_countries.zip
 
 .asv
 doc/source/getting_started/my_file.geojson
+
+.env
\ No newline at end of file
diff --git a/CHANGELOG.md b/CHANGELOG.md
index d631f664..7862f92d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,17 +4,18 @@
 
 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)
+- Added ``extract_unique_points`` method from shapely to GeoSeries/GeoDataframe (#2915).
+- 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).
 
 New features and improvements:
 
-Bug fixes:
+- Added ``exclusive`` parameter to ``sjoin_nearest`` method for Shapely >= 2.0 (#2877)
 
+Bug fixes:
 - ``assert_geodataframe_equal`` now handles GeoDataFrames with no active geometry (#2498)
-
 ## Version 0.13.2 (Jun 6, 2023)
 
 Bug fix:
@@ -880,4 +881,4 @@ Bug fixes:
 - Properly handle rows with missing geometries (#139, #193).
 - Fix ``GeoSeries.to_json`` (#263).
 - Correctly serialize metadata when pickling (#199, #206).
-- Fix ``merge`` and ``concat`` to return correct GeoDataFrame (#247, #320, #322).
+- Fix ``merge`` and ``concat`` to return correct GeoDataFrame (#247, #320, #322).
\ No newline at end of file
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 134a9b86..b6726b6b 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -92,6 +92,7 @@ Constructive methods and attributes
    GeoSeries.convex_hull
    GeoSeries.delaunay_triangles
    GeoSeries.envelope
+   GeoSeries.extract_unique_points
    GeoSeries.make_valid
    GeoSeries.minimum_bounding_circle
    GeoSeries.normalize
diff --git a/doc/source/docs/user_guide/geometric_manipulations.rst b/doc/source/docs/user_guide/geometric_manipulations.rst
index 1df5e5f2..8aa07a7d 100644
--- a/doc/source/docs/user_guide/geometric_manipulations.rst
+++ b/doc/source/docs/user_guide/geometric_manipulations.rst
@@ -50,6 +50,11 @@ Constructive methods
   smallest rectangular polygon (with sides parallel to the coordinate
   axes) that contains each object.
 
+.. method:: GeoSeries.extract_unique_points
+
+  Returns a :class:`~geopandas.GeoSeries` of geometries containing all distinct
+  vertices of each input geometry as a multipoint.
+  
 .. method:: GeoSeries.offset_curve(distance, quad_segs=8, join_style="round", mitre_limit=5.0)
 
   Returns a :class:`~geopandas.GeoSeries` containing a `Linestring` or `MultiLineString`
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index be4c1d54..e667584f 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -667,6 +667,18 @@ def exterior(data):
         return _unary_geo("exterior", data)
 
 
+def extract_unique_points(data):
+    if compat.USE_SHAPELY_20:
+        return shapely.extract_unique_points(data)
+    elif compat.USE_PYGEOS:
+        return pygeos.extract_unique_points(data)
+    else:
+        raise NotImplementedError(
+            f"shapely >= 2.0 or PyGEOS is required, "
+            f"version {shapely.__version__} is installed"
+        )
+
+
 def offset_curve(data, distance, quad_segs=8, join_style="round", mitre_limit=5.0):
     if compat.USE_SHAPELY_20:
         return shapely.offset_curve(
diff --git a/geopandas/array.py b/geopandas/array.py
index b9d2ca1f..fa860e62 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -533,6 +533,9 @@ class GeometryArray(ExtensionArray):
     def exterior(self):
         return GeometryArray(vectorized.exterior(self._data), crs=self.crs)
 
+    def extract_unique_points(self):
+        return GeometryArray(vectorized.extract_unique_points(self._data), crs=self.crs)
+
     def offset_curve(self, distance, quad_segs=8, join_style="round", mitre_limit=5.0):
         return GeometryArray(
             vectorized.offset_curve(
diff --git a/geopandas/base.py b/geopandas/base.py
index d15cecaf..8f0eb4fc 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -733,10 +733,41 @@ GeometryCollection
         # TODO: return empty geometry for non-polygons
         return _delegate_property("exterior", self)
 
+    def extract_unique_points(self):
+        """Returns a ``GeoSeries`` of MultiPoints representing all
+        distinct vertices of an input geometry.
+
+        Examples
+        --------
+
+        >>> from shapely import LineString, Polygon
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         LineString([(0, 0), (0, 0), (1, 1), (1, 1)]),
+        ...         Polygon([(0, 0), (0, 0), (1, 1), (1, 1)])
+        ...     ],
+        ...     crs=3857
+        ... )
+        >>> s
+        0    LINESTRING (0.000 0.000, 0.000 0.000, 1.000 1....
+        1    POLYGON ((0.000 0.000, 0.000 0.000, 1.000 1.00...
+        dtype: geometry
+
+        >>> s.extract_unique_points()
+        0    MULTIPOINT (0.000 0.000, 1.000 1.000)
+        1    MULTIPOINT (0.000 0.000, 1.000 1.000)
+        dtype: geometry
+
+        See also
+        --------
+
+        GeoSeries.get_coordinates : extract coordinates as a :class:`~pandas.DataFrame`
+        """
+        return _delegate_geo_method("extract_unique_points", self)
+
     def offset_curve(self, distance, quad_segs=8, join_style="round", mitre_limit=5.0):
         """Returns a ``LineString`` or ``MultiLineString`` geometry at a
         distance from the object on its right or its left side.
-
         Parameters
         ----------
         distance : float | array-like
diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py
index 038656ed..a5ac58a8 100644
--- a/geopandas/geodataframe.py
+++ b/geopandas/geodataframe.py
@@ -2184,6 +2184,7 @@ individually so that features may have different properties
         lsuffix="left",
         rsuffix="right",
         distance_col=None,
+        exclusive=False,
     ):
         """
         Spatial join of two GeoDataFrames based on the distance between their
@@ -2221,6 +2222,10 @@ individually so that features may have different properties
         distance_col : string, default None
             If set, save the distances computed between matching geometries under a
             column of this name in the joined GeoDataFrame.
+        exclusive : bool, optional, default False
+            If True, the nearest geometries that are equal to the input geometry
+            will not be returned, default False.
+            Requires Shapely >= 2.0
 
         Examples
         --------
@@ -2303,6 +2308,7 @@ chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
             lsuffix=lsuffix,
             rsuffix=rsuffix,
             distance_col=distance_col,
+            exclusive=exclusive,
         )
 
     def clip(self, mask, keep_geom_type=False):
diff --git a/geopandas/sindex.py b/geopandas/sindex.py
index 7334f4c3..0ef0eeae 100644
--- a/geopandas/sindex.py
+++ b/geopandas/sindex.py
@@ -232,7 +232,12 @@ class BaseSpatialIndex:
         raise NotImplementedError
 
     def nearest(
-        self, geometry, return_all=True, max_distance=None, return_distance=False
+        self,
+        geometry,
+        return_all=True,
+        max_distance=None,
+        return_distance=False,
+        exclusive=False,
     ):
         """
         Return the nearest geometry in the tree for each input geometry in
@@ -281,6 +286,9 @@ geometries}
             Must be greater than 0. By default None, indicating no distance limit.
         return_distance : bool, optional
             If True, will return distances in addition to indexes. By default False
+        exclusive : bool, optional
+            if True, the nearest geometries that are equal to the input geometry
+            will not be returned. By default False.  Requires Shapely >= 2.0.
 
         Returns
         -------
@@ -839,13 +847,23 @@ if compat.SHAPELY_GE_20 or compat.HAS_PYGEOS:
 
         @doc(BaseSpatialIndex.nearest)
         def nearest(
-            self, geometry, return_all=True, max_distance=None, return_distance=False
+            self,
+            geometry,
+            return_all=True,
+            max_distance=None,
+            return_distance=False,
+            exclusive=False,
         ):
             if not (compat.USE_SHAPELY_20 or compat.PYGEOS_GE_010):
                 raise NotImplementedError(
                     "sindex.nearest requires shapely >= 2.0 or pygeos >= 0.10"
                 )
 
+            if exclusive and not compat.USE_SHAPELY_20:
+                raise NotImplementedError(
+                    "sindex.nearest exclusive parameter requires shapely >= 2.0"
+                )
+
             geometry = self._as_geometry_array(geometry)
             if isinstance(geometry, BaseGeometry) or geometry is None:
                 geometry = [geometry]
@@ -856,6 +874,7 @@ if compat.SHAPELY_GE_20 or compat.HAS_PYGEOS:
                     max_distance=max_distance,
                     return_distance=return_distance,
                     all_matches=return_all,
+                    exclusive=exclusive,
                 )
             else:
                 if not return_all and max_distance is None and not return_distance:
diff --git a/geopandas/tools/sjoin.py b/geopandas/tools/sjoin.py
index 3b7b01eb..8bffdd16 100644
--- a/geopandas/tools/sjoin.py
+++ b/geopandas/tools/sjoin.py
@@ -356,12 +356,14 @@ def _nearest_query(
     max_distance: float,
     how: str,
     return_distance: bool,
+    exclusive: bool,
 ):
     if not (compat.USE_SHAPELY_20 or (compat.USE_PYGEOS and compat.PYGEOS_GE_010)):
         raise NotImplementedError(
             "Currently, only PyGEOS >= 0.10.0 or Shapely >= 2.0 supports "
             "`nearest_all`. " + compat.INSTALL_PYGEOS_ERROR
         )
+
     # use the opposite of the join direction for the index
     use_left_as_sindex = how == "right"
     if use_left_as_sindex:
@@ -376,6 +378,7 @@ def _nearest_query(
             return_all=True,
             max_distance=max_distance,
             return_distance=return_distance,
+            exclusive=exclusive,
         )
         if return_distance:
             (input_idx, tree_idx), distances = res
@@ -409,6 +412,7 @@ def sjoin_nearest(
     lsuffix: str = "left",
     rsuffix: str = "right",
     distance_col: Optional[str] = None,
+    exclusive: bool = False,
 ) -> GeoDataFrame:
     """Spatial join of two GeoDataFrames based on the distance between their geometries.
 
@@ -446,6 +450,10 @@ def sjoin_nearest(
     distance_col : string, default None
         If set, save the distances computed between matching geometries under a
         column of this name in the joined GeoDataFrame.
+    exclusive : bool, default False
+        If True, the nearest geometries that are equal to the input geometry
+        will not be returned, default False.
+        Requires Shapely >= 2.0.
 
     Examples
     --------
@@ -527,7 +535,9 @@ chicago_w_groceries[chicago_w_groceries["community"] == "UPTOWN"]
 
     return_distance = distance_col is not None
 
-    join_df = _nearest_query(left_df, right_df, max_distance, how, return_distance)
+    join_df = _nearest_query(
+        left_df, right_df, max_distance, how, return_distance, exclusive
+    )
 
     if return_distance:
         join_df = join_df.rename(columns={"distances": distance_col})


Test patch:
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index 83879ff5..0e907b0e 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -49,6 +49,7 @@ class TestGeomMethods:
         self.sqz = Polygon([(1, 1, 1), (2, 2, 2), (3, 3, 3), (4, 4, 4)])
         self.t4 = Polygon([(0, 0), (3, 0), (3, 3), (0, 2)])
         self.t5 = Polygon([(2, 0), (3, 0), (3, 3), (2, 3)])
+        self.t6 = Polygon([(2, 0), (2, 0), (3, 0), (3, 0)])
         self.inner_sq = Polygon(
             [(0.25, 0.25), (0.75, 0.25), (0.75, 0.75), (0.25, 0.75)]
         )
@@ -1124,6 +1125,27 @@ class TestGeomMethods:
         assert isinstance(e, GeoSeries)
         assert self.g3.crs == e.crs
 
+    @pytest.mark.skipif(
+        not (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason=(
+            "extract_unique_points is only implemented for pygeos and shapely >= 2.0"
+        ),
+    )
+    def test_extract_unique_points(self):
+        eup = GeoSeries([self.t6]).extract_unique_points()
+        expected = GeoSeries([MultiPoint([(2, 0), (3, 0)])])
+        assert_series_equal(eup, expected)
+
+    @pytest.mark.skipif(
+        (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="extract_unique_points not implemented for shapely<2",
+    )
+    def test_extract_unique_points_not_implemented(self):
+        with pytest.raises(
+            NotImplementedError, match="shapely >= 2.0 or PyGEOS is required"
+        ):
+            self.g1.extract_unique_points()
+
     @pytest.mark.skipif(
         not (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
         reason="minimum_bounding_circle is only implemented for pygeos, not shapely",
diff --git a/geopandas/tests/test_sindex.py b/geopandas/tests/test_sindex.py
index e60127cb..bc909a33 100644
--- a/geopandas/tests/test_sindex.py
+++ b/geopandas/tests/test_sindex.py
@@ -834,6 +834,66 @@ class TestPygeosInterface:
         else:
             assert_array_equal(res, expected[0])
 
+    @pytest.mark.skipif(
+        not (compat.USE_SHAPELY_20),
+        reason=(
+            "shapely >= 2.0 is required to test sindex.nearest with parameter exclusive"
+        ),
+    )
+    @pytest.mark.parametrize("return_distance", [True, False])
+    @pytest.mark.parametrize(
+        "return_all,max_distance,exclusive,expected",
+        [
+            (False, None, False, ([[0, 1, 2, 3, 4], [0, 1, 2, 3, 4]], 5 * [0])),
+            (False, None, True, ([[0, 1, 2, 3, 4], [1, 0, 1, 2, 3]], 5 * [sqrt(2)])),
+            (True, None, False, ([[0, 1, 2, 3, 4], [0, 1, 2, 3, 4]], 5 * [0])),
+            (
+                True,
+                None,
+                True,
+                ([[0, 1, 1, 2, 2, 3, 3, 4], [1, 0, 2, 1, 3, 2, 4, 3]], 8 * [sqrt(2)]),
+            ),
+            (False, 1.1, True, ([[1, 2, 5], [5, 5, 1]], 3 * [1])),
+            (True, 1.1, True, ([[1, 2, 5, 5], [5, 5, 1, 2]], 4 * [1])),
+        ],
+    )
+    def test_nearest_exclusive(
+        self, expected, max_distance, return_all, return_distance, exclusive
+    ):
+        geoms = mod.points(np.arange(5), np.arange(5))
+        if max_distance:
+            # add a non grid point
+            geoms = np.append(geoms, [Point(1, 2)])
+
+        df = geopandas.GeoDataFrame({"geometry": geoms})
+
+        ps = geoms
+        res = df.sindex.nearest(
+            ps,
+            return_all=return_all,
+            max_distance=max_distance,
+            return_distance=return_distance,
+            exclusive=exclusive,
+        )
+        if return_distance:
+            assert_array_equal(res[0], expected[0])
+            assert_array_equal(res[1], expected[1])
+        else:
+            assert_array_equal(res, expected[0])
+
+    @pytest.mark.skipif(
+        compat.USE_SHAPELY_20 or not (compat.USE_PYGEOS and not compat.PYGEOS_GE_010),
+        reason="sindex.nearest exclusive parameter requires shapely >= 2.0",
+    )
+    def test_nearest_exclusive_unavailable(self):
+        from shapely.geometry import Point
+
+        geoms = [Point((x, y)) for (x, y) in zip(np.arange(5), np.arange(5))]
+        df = geopandas.GeoDataFrame(geometry=geoms)
+
+        with pytest.raises(NotImplementedError, match="requires shapely >= 2.0"):
+            df.sindex.nearest(geoms, exclusive=True)
+
     # --------------------------- misc tests ---------------------------- #
 
     def test_empty_tree_geometries(self):
diff --git a/geopandas/tools/tests/test_sjoin.py b/geopandas/tools/tests/test_sjoin.py
index 1b15d21e..a2cb087c 100644
--- a/geopandas/tools/tests/test_sjoin.py
+++ b/geopandas/tools/tests/test_sjoin.py
@@ -1,18 +1,18 @@
 import math
 from typing import Sequence
-from geopandas.testing import assert_geodataframe_equal
 
 import numpy as np
 import pandas as pd
+import shapely
 
 from shapely.geometry import Point, Polygon, GeometryCollection
 
 import geopandas
 import geopandas._compat as compat
 from geopandas import GeoDataFrame, GeoSeries, read_file, sjoin, sjoin_nearest
-from geopandas.testing import assert_geoseries_equal
+from geopandas.testing import assert_geodataframe_equal, assert_geoseries_equal
 
-from pandas.testing import assert_frame_equal
+from pandas.testing import assert_frame_equal, assert_series_equal
 import pytest
 
 
@@ -926,3 +926,31 @@ class TestNearest:
         result5 = result5.dropna()
         result5["index_right"] = result5["index_right"].astype("int64")
         assert_geodataframe_equal(result5, result4, check_like=True)
+
+    @pytest.mark.skipif(
+        not (compat.USE_SHAPELY_20),
+        reason=(
+            "shapely >= 2.0 is required to run sjoin_nearest"
+            "with parameter `exclusive` set"
+        ),
+    )
+    @pytest.mark.parametrize(
+        "max_distance,expected", [(None, [1, 1, 3, 3, 2]), (1.1, [3, 3, 1, 2])]
+    )
+    def test_sjoin_nearest_exclusive(self, max_distance, expected):
+        geoms = shapely.points(np.arange(3), np.arange(3))
+        geoms = np.append(geoms, [Point(1, 2)])
+
+        df = geopandas.GeoDataFrame({"geometry": geoms})
+        result = df.sjoin_nearest(
+            df, max_distance=max_distance, distance_col="dist", exclusive=True
+        )
+
+        assert_series_equal(
+            result["index_right"].reset_index(drop=True),
+            pd.Series(expected),
+            check_names=False,
+        )
+
+        if max_distance:
+            assert result["dist"].max() <= max_distance