Repository: geopandas/geopandas
Pull request #2960 (2023-07-14 23:46:36)

ENH: add shortest_line from shapely


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 63689cb2..271c0df6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@ New methods:
 - 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 ``shortest_line`` method from shapely to GeoSeries/GeoDataframe (#2960).
 
 New features and improvements:
 
@@ -20,7 +21,7 @@ Bug fixes:
 
 - Fix a color assignment in ``explore`` when using ``UserDefined`` bins (#2923)
 - ``assert_geodataframe_equal`` now handles GeoDataFrames with no active geometry (#2498)
-- Fix bug in `apply` with `axis=1` where the given user defined function returns nested 
+- Fix bug in `apply` with `axis=1` where the given user defined function returns nested
   data in the geometry column (#2959)
 - Properly infer schema for np.int32 and pd.Int32Dtype columns (#2950)
 
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 755b56a2..0b9872ca 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -98,6 +98,7 @@ Constructive methods and attributes
    GeoSeries.normalize
    GeoSeries.sample_points
    GeoSeries.segmentize
+   GeoSeries.shortest_line
    GeoSeries.simplify
 
 Affine transformations
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index 596df147..9924dc60 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -941,6 +941,18 @@ def union(data, other):
         return _binary_geo("union", data, other)
 
 
+def shortest_line(data, other):
+    if compat.USE_SHAPELY_20:
+        return shapely.shortest_line(data, other)
+    elif compat.USE_PYGEOS:
+        return _binary_method("shortest_line", data, other)
+    else:
+        raise NotImplementedError(
+            f"shapely >= 2.0 or PyGEOS is required, "
+            f"version {shapely.__version__} is installed"
+        )
+
+
 #
 # Other operations
 #
diff --git a/geopandas/array.py b/geopandas/array.py
index 6869b27f..e1892bab 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -658,6 +658,11 @@ class GeometryArray(ExtensionArray):
     def union(self, other):
         return GeometryArray(self._binary_method("union", self, other), crs=self.crs)
 
+    def shortest_line(self, other):
+        return GeometryArray(
+            self._binary_method("shortest_line", self, other), crs=self.crs
+        )
+
     #
     # Other operations
     #
diff --git a/geopandas/base.py b/geopandas/base.py
index 1fc7a038..fd8cc1b0 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -3099,6 +3099,109 @@ GeometryCollection
         clipped_geometry = geometry_array.clip_by_rect(xmin, ymin, xmax, ymax)
         return GeoSeries(clipped_geometry, index=self.index, crs=self.crs)
 
+    def shortest_line(self, other, align=True):
+        """
+        Returns the shortest two-point line between two geometries.
+
+        The resulting line consists of two points, representing the nearest points
+        between the geometry pair. The line always starts in the first geometry a
+        and ends in he second geometry b. The endpoints of the line will not
+        necessarily be existing vertices of the input geometries a and b, but
+        can also be a point along a line segment.
+
+
+        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
+            shortest line with.
+        align : bool (default True)
+            If True, automatically aligns GeoSeries based on their indices.
+            If False, the order of elements is preserved.
+
+        Returns
+        -------
+        GeoSeries
+
+        Examples
+        --------
+        >>> from shapely.geometry import Polygon, LineString, Point
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (2, 2), (0, 2)]),
+        ...         Polygon([(0, 0), (2, 2), (0, 2)]),
+        ...         LineString([(0, 0), (2, 2)]),
+        ...         LineString([(2, 0), (0, 2)]),
+        ...         Point(0, 1),
+        ...     ],
+        ...     crs=5514
+        ... )
+        >>> s
+        0    POLYGON ((0.000 0.000, 2.000 2.000, 0.000 2.00...
+        1    POLYGON ((0.000 0.000, 2.000 2.000, 0.000 2.00...
+        2                LINESTRING (0.000 0.000, 2.000 2.000)
+        3                LINESTRING (2.000 0.000, 0.000 2.000)
+        4                                  POINT (0.000 1.000)
+        dtype: geometry
+
+        We can also do intersection of each geometry and a single
+        shapely geometry:
+
+        .. image:: ../../../_static/binary_op-03.svg
+           :align: center
+
+        >>> p = Point(3, 3)
+        >>> s.shortest_line(p)
+        0    LINESTRING (2.000 2.000, 3.000 3.000)
+        1    LINESTRING (2.000 2.000, 3.000 3.000)
+        2    LINESTRING (2.000 2.000, 3.000 3.000)
+        3    LINESTRING (1.000 1.000, 3.000 3.000)
+        4    LINESTRING (0.000 1.000, 3.000 3.000)
+        dtype: geometry
+
+        We can also check two GeoSeries against each other, row by row.
+        The GeoSeries above have different indices than the one below. We can either
+        align both GeoSeries based on index values and compare elements with the same
+        index using ``align=True`` or ignore index and compare elements based on their
+        matching order using ``align=False``:
+
+        .. image:: ../../../_static/binary_op-02.svg
+
+        >>> 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(10, 15),
+        ...         Point(0, 1),
+        ...     ],
+        ...     index=range(1, 6),
+        ...     crs=5514,
+        ... )
+        >>> s.shortest_line(s2, align=True)
+        0                                       None
+        1      LINESTRING (0.500 0.500, 0.500 0.500)
+        2      LINESTRING (2.000 2.000, 3.000 1.000)
+        3      LINESTRING (2.000 0.000, 2.000 0.000)
+        4    LINESTRING (0.000 1.000, 10.000 15.000)
+        5                                       None
+        dtype: geometry
+
+        >>> s.shortest_line(s2, align=False)
+        0      LINESTRING (0.500 0.500, 0.500 0.500)
+        1      LINESTRING (2.000 2.000, 3.000 1.000)
+        2      LINESTRING (0.500 0.500, 1.000 0.000)
+        3    LINESTRING (0.000 2.000, 10.000 15.000)
+        4      LINESTRING (0.000 1.000, 0.000 1.000)
+        dtype: geometry
+        """
+        return _binary_geo("shortest_line", self, other, align)
+
     #
     # Other operations
     #


Test patch:
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index b7a753b7..b66ce57b 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -368,6 +368,40 @@ class TestGeomMethods:
         expected = GeoSeries([self.t1, self.t1])
         self._test_binary_topological("difference", expected, self.g1, self.t2)
 
+    @pytest.mark.skipif(
+        not (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="get_coordinates not implemented for shapely<2",
+    )
+    def test_shortest_line(self):
+        expected = GeoSeries([LineString([(1, 1), (5, 5)]), None])
+        assert_array_dtype_equal(expected, self.na_none.shortest_line(self.p0))
+
+        expected = GeoSeries(
+            [
+                LineString([(5, 5), (1, 1)]),
+                LineString([(2, 0), (2, 0)]),
+            ]
+        )
+        assert_array_dtype_equal(expected, self.g6.shortest_line(self.g7))
+
+        expected = GeoSeries(
+            [LineString([(0.5, 0.5), (0.5, 0.5)]), LineString([(0.5, 0.5), (0.5, 0.5)])]
+        )
+        crossed_lines_inv = self.crossed_lines[::-1]
+        assert_array_dtype_equal(
+            expected, self.crossed_lines.shortest_line(crossed_lines_inv, align=False)
+        )
+
+    @pytest.mark.skipif(
+        (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="get_coordinates not implemented for shapely<2",
+    )
+    def test_shortest_line_not(self):
+        with pytest.raises(
+            NotImplementedError, match="shapely >= 2.0 or PyGEOS is required"
+        ):
+            self.na_none.shortest_line(self.p0)
+
     def test_geo_op_empty_result(self):
         l1 = LineString([(0, 0), (1, 1)])
         l2 = LineString([(2, 2), (3, 3)])