Repository: geopandas/geopandas
Pull request #2929 (2023-06-13 03:25:05)

ENH: add frechet_distance method
adding frechet distance method for geopandas from shapely


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 63689cb2..a3a450ef 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,25 +2,34 @@
 
 ## Development version
 
+GeoPandas will use Shapely 2.0 by default instead of PyGEOS when both Shapely >= 2.0 and
+PyGEOS are installed.  PyGEOS will continue to be used by default when PyGEOS is
+installed alongside Shapely < 2.0.  Support for PyGEOS and Shapely < 2.0 will be removed
+in GeoPandas 1.0.
+
 New methods:
 
+- Added ``frechet_distance()`` method from shapely to GeoSeries/GeoDataframe (#2929).
 - Added ``segmentize`` method from shapely to GeoSeries/GeoDataFrame (#2910).
 - 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).
+- Added ``shortest_line`` method from shapely to GeoSeries/GeoDataframe (#2960).
+- Added ``minimum_rotated_rectangle`` method from shapely to GeoSeries/GeoDataframe (#2541).
 
 New features and improvements:
 
 - Added ``exclusive`` parameter to ``sjoin_nearest`` method for Shapely >= 2.0 (#2877)
 
 Bug fixes:
+
 - Fix ambiguous error when GeoDataFrame is initialised with a column called "crs" (#2944)
 
 - 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)
 
@@ -889,4 +898,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).
\ No newline at end of file
+- Fix ``merge`` and ``concat`` to return correct GeoDataFrame (#247, #320, #322).
diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst
index 755b56a2..b76c60d6 100644
--- a/doc/source/docs/reference/geoseries.rst
+++ b/doc/source/docs/reference/geoseries.rst
@@ -25,6 +25,7 @@ General methods and attributes
    GeoSeries.offset_curve
    GeoSeries.distance
    GeoSeries.hausdorff_distance
+   GeoSeries.frechet_distance
    GeoSeries.representative_point
    GeoSeries.exterior
    GeoSeries.interiors
@@ -95,9 +96,11 @@ Constructive methods and attributes
    GeoSeries.extract_unique_points
    GeoSeries.make_valid
    GeoSeries.minimum_bounding_circle
+   GeoSeries.minimum_rotated_rectangle
    GeoSeries.normalize
    GeoSeries.sample_points
    GeoSeries.segmentize
+   GeoSeries.shortest_line
    GeoSeries.simplify
 
 Affine transformations
diff --git a/doc/source/docs/reference/sindex.rst b/doc/source/docs/reference/sindex.rst
index 5536d6d9..59b88ed6 100644
--- a/doc/source/docs/reference/sindex.rst
+++ b/doc/source/docs/reference/sindex.rst
@@ -3,10 +3,11 @@ Spatial index
 =============
 .. currentmodule:: geopandas
 
-GeoPandas offers built-in support for spatial indexing using an R-Tree algorithm.
-Depending on the ability to import PyGEOS, GeoPandas will either use
-``pygeos.STRtree`` or ``rtree.index.Index``. The main interface for both is the
-same and follows the PyGEOS model.
+GeoPandas will use the STRtree implementation provided by the default spatial engine
+(``shapely.STRtree`` for Shapely > 2.0, ``pygeos.STRtree`` for Shapely < 2.0 if PyGEOS
+is installed) if compatible, or ``rtree.index.Index`` for Shapely < 2.0 when PyGEOS is
+unavailable.
+
 
 ``GeoSeries.sindex`` creates a spatial index, which can use the methods and
 properties documented below.
@@ -36,7 +37,8 @@ methods:
     valid_query_predicates
 
 The concrete implementations currently available are
-``geopandas.sindex.PyGEOSSTRTreeIndex`` and ``geopandas.sindex.RTreeIndex``.
+``geopandas.sindex.PyGEOSSTRTreeIndex`` (when using Shapely 2.0 or PyGEOS) and
+``geopandas.sindex.RTreeIndex`` (when using Shapely <2).
 
 In addition to the methods listed above, the ``rtree``-based spatial index
 (``geopandas.sindex.RTreeIndex``) offers the full capability of
@@ -44,8 +46,9 @@ In addition to the methods listed above, the ``rtree``-based spatial index
 
 Similarly, the PyGEOS-based spatial index
 (``geopandas.sindex.PyGEOSSTRTreeIndex``) offers the full capability of
-``pygeos.STRtree``, including nearest-neighbor queries.
-See the full API in the `PyGEOS STRTree documentation`_.
+``pygeos.STRtree/shapely.STRtree``, including nearest-neighbor queries.
+See the full API in the `Shapely STRTree documentation`_ or `PyGEOS STRTree documentation`_.
 
 .. _rtree documentation: https://rtree.readthedocs.io/en/stable/class.html
 .. _PyGEOS STRTree documentation: https://pygeos.readthedocs.io/en/latest/strtree.html
+.. _Shapely STRTree documentation: https://shapely.readthedocs.io/en/stable/strtree.html
\ No newline at end of file
diff --git a/doc/source/docs/user_guide/pygeos_to_shapely.rst b/doc/source/docs/user_guide/pygeos_to_shapely.rst
index 94887abb..c63f22a6 100644
--- a/doc/source/docs/user_guide/pygeos_to_shapely.rst
+++ b/doc/source/docs/user_guide/pygeos_to_shapely.rst
@@ -7,7 +7,7 @@ Shapely project and released as part of Shapely 2.0. GeoPandas will therefore
 deprecate support of the PyGEOS backend and will go forward with Shapely 2.0 as the
 only geometry engine exposing GEOS functionality.
 
-Given that historically the PyGEOS engine is automatically used if the package is installed (this behaviour will change in GeoPandas 0.14 where Shapely 2.0 will be used by default if installed), some downstream code may depend on
+Given that historically the PyGEOS engine was automatically used if the package is installed (this behaviour will changed in GeoPandas 0.14 where Shapely 2.0 is used by default if installed), some downstream code may depend on
 PyGEOS geometries being available as underlying data of a ``GeometryArray``.
 
 This guide outlines the migration from the PyGEOS-based code to the Shapely-based code.
@@ -27,9 +27,9 @@ GeoPandas 0.13
 GeoPandas 0.14
 ^^^^^^^^^^^^^^
 
-- The default backend will become Shapely 2.0 and the PyGEOS will be used only
-  if Shapely 1.8 is installed instead of 2.0 or newer. The PyGEOS backend will still be
-  supported, but a user will need to opt in using the environment variable
+- The default backend is Shapely 2.0 and the PyGEOS is used only
+  if Shapely 1.8 is installed instead of 2.0 or newer. The PyGEOS backend is still
+  supported, but a user needs to opt in using the environment variable
   ``USE_PYGEOS`` as explained in the
   `installation instructions <../getting_started/install>`__.
 
diff --git a/doc/source/getting_started/install.rst b/doc/source/getting_started/install.rst
index 23883b41..94feeafe 100644
--- a/doc/source/getting_started/install.rst
+++ b/doc/source/getting_started/install.rst
@@ -164,17 +164,12 @@ Using the optional PyGEOS dependency
 
 .. attention::
 
-  The upcoming Shapely 2.0 release will absorb all improvements from PyGEOS.
-  If you are considering trying out those improvements, you can also test
-  the prerelease of Shapely instead.
-  See https://shapely.readthedocs.io/en/latest/release/2.x.html#version-2-0-0
-  for the release notes of Shapely 2.0, and https://github.com/shapely/shapely/discussions/1464
-  on how to install this and give feedback.
-
-Work is ongoing to improve the performance of GeoPandas. Currently, the
-fast implementations of basic spatial operations live in the `PyGEOS`_
-package (but work is under way to contribute those improvements to Shapely,
-coming to Shapely 2.0).
+  The Shapely 2.0 release absorbed all improvements from PyGEOS and the support of
+  PyGEOS in GeoPandas is deprecated. Please use Shapely 2 instead.
+
+Historically, the
+fast implementations of basic spatial operations used to live in the `PyGEOS`_
+package.
 Starting with GeoPandas 0.8, it is possible to optionally use those
 experimental speedups by installing PyGEOS. This can be done with conda
 (using the conda-forge channel) or pip::
@@ -184,11 +179,10 @@ experimental speedups by installing PyGEOS. This can be done with conda
     # pip
     pip install pygeos
 
-More specifically, whether the speedups are used or not is determined by:
+More specifically, whether the PyGEOS is used or not is determined by:
 
-- If PyGEOS >= 0.8 is installed, it will be used by default (but installing
-  GeoPandas will not yet automatically install PyGEOS as dependency, you need
-  to do this manually).
+- If Shapely >= 2 is installed, it will be used by default.
+- If PyGEOS >= 0.8 and Shapely < 2 are installed, PyGEOS will be used by default
 
 - You can still toggle the use of PyGEOS when it is available, by:
 
@@ -207,19 +201,18 @@ More specifically, whether the speedups are used or not is determined by:
     although this variable can be set during an interactive session, it will
     only work if the GeoDataFrames you use are created (e.g. reading a file
     with ``read_file``) after changing this value.
-    Attention: changing this option will no longer work in all cases when
-    having Shapely >=2.0 installed. In that case, use the environment variable
-    (see option above).
+
+.. attention::
+
+    Changing this option no longer works in all cases when
+    Shapely >=2.0 is installed. In that case, use the environment variable
+    (see option above) or simply use Shapely.
 
 .. warning::
 
-    The use of PyGEOS is experimental! Although it is passing all tests,
-    there might still be issues and not all functions of GeoPandas will
-    already benefit from speedups (one known issue: the `to_crs` coordinate
-    transformations lose the z coordinate). But trying this out is very welcome!
-    Any issues you encounter (but also reports of successful usage are
-    interesting!) can be reported at https://gitter.im/geopandas/geopandas
-    or https://github.com/geopandas/geopandas/issues
+    The use of PyGEOS is deprecated! If you depend on PyGEOS, see
+    the :doc:`migration guide <../docs/user_guide/pygeos_to_shapely>`
+    to Shapely 2.0. Otherwise, use Shapely 2 directly.
 
 
 .. _PyPI: https://pypi.python.org/pypi/geopandas
diff --git a/geopandas/_compat.py b/geopandas/_compat.py
index 3ccf3490..9035f010 100644
--- a/geopandas/_compat.py
+++ b/geopandas/_compat.py
@@ -21,7 +21,7 @@ PANDAS_GE_13 = Version(pd.__version__) >= Version("1.3.0")
 PANDAS_GE_14 = Version(pd.__version__) >= Version("1.4.0rc0")
 PANDAS_GE_15 = Version(pd.__version__) >= Version("1.5.0")
 PANDAS_GE_20 = Version(pd.__version__) >= Version("2.0.0")
-PANDAS_GE_21 = Version(pd.__version__) >= Version("2.1.0.dev0")
+PANDAS_GE_21 = Version(pd.__version__) >= Version("2.1.0")
 
 
 # -----------------------------------------------------------------------------
@@ -88,7 +88,10 @@ def set_use_pygeos(val=None):
         USE_PYGEOS = bool(val)
     else:
         if USE_PYGEOS is None:
-            USE_PYGEOS = HAS_PYGEOS
+            if SHAPELY_GE_20:
+                USE_PYGEOS = False
+            else:
+                USE_PYGEOS = HAS_PYGEOS
 
             if env_use_pygeos is not None:
                 USE_PYGEOS = bool(int(env_use_pygeos))
@@ -138,23 +141,12 @@ def set_use_pygeos(val=None):
         except ImportError:
             raise ImportError(INSTALL_PYGEOS_ERROR)
 
-    if USE_PYGEOS and env_use_pygeos is None and SHAPELY_GE_20:
+    if USE_PYGEOS:
         warnings.warn(
-            "Shapely 2.0 is installed, but because PyGEOS is also installed, "
-            "GeoPandas still uses PyGEOS by default. However, starting with "
-            "version 0.14, the default will switch to Shapely. To force to "
-            "use Shapely 2.0 now, you can either uninstall PyGEOS or set the "
-            "environment variable USE_PYGEOS=0. "
-            "You can do this before starting the Python process, or in your code "
-            "before importing geopandas:"
-            "\n\nimport os\nos.environ['USE_PYGEOS'] = '0'\nimport geopandas\n\n"
-            "In the next release, GeoPandas will switch to using Shapely by default, "
-            "even if PyGEOS is installed. If you only have PyGEOS installed to "
-            "get speed-ups, this switch should be smooth. However, if you are "
-            "using PyGEOS directly (calling PyGEOS functions on geometries "
-            "from GeoPandas), this will then stop working and you are encouraged to "
-            "migrate from PyGEOS to Shapely 2.0 "
-            "(https://shapely.readthedocs.io/en/latest/migration_pygeos.html).",
+            "GeoPandas is set to use PyGEOS over Shapely. PyGEOS support is deprecated"
+            "and will be removed in GeoPandas 1.0, released in the Q1 of 2024. "
+            "Please migrate to Shapely 2.0 "
+            "(https://geopandas.org/en/stable/docs/user_guide/pygeos_to_shapely.html).",
             DeprecationWarning,
             stacklevel=6,
         )
diff --git a/geopandas/_vectorized.py b/geopandas/_vectorized.py
index 596df147..77d7ab5e 100644
--- a/geopandas/_vectorized.py
+++ b/geopandas/_vectorized.py
@@ -658,6 +658,15 @@ def envelope(data):
         return _unary_geo("envelope", data)
 
 
+def minimum_rotated_rectangle(data):
+    if compat.USE_SHAPELY_20:
+        return shapely.oriented_envelope(data)
+    elif compat.USE_PYGEOS:
+        return pygeos.oriented_envelope(data)
+    else:
+        return _unary_geo("minimum_rotated_rectangle", data)
+
+
 def exterior(data):
     if compat.USE_SHAPELY_20:
         return shapely.get_exterior_ring(data)
@@ -941,6 +950,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
 #
@@ -969,6 +990,20 @@ def hausdorff_distance(data, other, densify=None, **kwargs):
         )
 
 
+def frechet_distance(data, other, densify=None, **kwargs):
+    if compat.USE_SHAPELY_20:
+        return shapely.frechet_distance(data, other, densify=densify, **kwargs)
+    elif compat.USE_PYGEOS:
+        return _binary_method(
+            "frechet_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 6869b27f..9ea601cb 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -529,6 +529,11 @@ class GeometryArray(ExtensionArray):
     def envelope(self):
         return GeometryArray(vectorized.envelope(self._data), crs=self.crs)
 
+    def minimum_rotated_rectangle(self):
+        return GeometryArray(
+            vectorized.minimum_rotated_rectangle(self.data), crs=self.crs
+        )
+
     @property
     def exterior(self):
         return GeometryArray(vectorized.exterior(self._data), crs=self.crs)
@@ -658,6 +663,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
     #
@@ -670,6 +680,10 @@ class GeometryArray(ExtensionArray):
         self.check_geographic_crs(stacklevel=6)
         return self._binary_method("hausdorff_distance", self, other, **kwargs)
 
+    def frechet_distance(self, other, **kwargs):
+        self.check_geographic_crs(stacklevel=6)
+        return self._binary_method("frechet_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)
@@ -1066,7 +1080,7 @@ class GeometryArray(ExtensionArray):
         self._data[idx] = value_arr
         return self
 
-    def fillna(self, value=None, method=None, limit=None):
+    def fillna(self, value=None, method=None, limit=None, copy=True):
         """
         Fill NA values with geometry (or geometries) or using the specified method.
 
@@ -1090,6 +1104,10 @@ class GeometryArray(ExtensionArray):
             maximum number of entries along the entire axis where NaNs will be
             filled.
 
+        copy : bool, default True
+            Whether to make a copy of the data before filling. If False, then
+            the original should be modified and no new memory should be allocated.
+
         Returns
         -------
         GeometryArray
@@ -1098,7 +1116,10 @@ class GeometryArray(ExtensionArray):
             raise NotImplementedError("fillna with a method is not yet supported")
 
         mask = self.isna()
-        new_values = self.copy()
+        if copy:
+            new_values = self.copy()
+        else:
+            new_values = self
         return new_values._fill(mask, value) if mask.any() else new_values
 
     def astype(self, dtype, copy=True):
diff --git a/geopandas/base.py b/geopandas/base.py
index 1fc7a038..af15c474 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -697,6 +697,46 @@ GeometryCollection
         """
         return _delegate_property("envelope", self)
 
+    def minimum_rotated_rectangle(self):
+        """Returns a ``GeoSeries`` of the general minimum bounding rectangle
+        that contains the object.
+
+        Unlike envelope this rectangle is not constrained to be parallel
+        to the coordinate axes. If the convex hull of the object is a
+        degenerate (line or point) this degenerate is returned.
+
+        Examples
+        --------
+
+        >>> from shapely.geometry import Polygon, LineString, Point, MultiPoint
+        >>> s = geopandas.GeoSeries(
+        ...     [
+        ...         Polygon([(0, 0), (1, 1), (0, 1)]),
+        ...         LineString([(0, 0), (1, 1), (1, 0)]),
+        ...         MultiPoint([(0, 0), (1, 1)]),
+        ...         Point(0, 0),
+        ...     ]
+        ... )
+        >>> s
+        0    POLYGON ((0.00000 0.00000, 1.00000 1.00000, 0....
+        1    LINESTRING (0.00000 0.00000, 1.00000 1.00000, ...
+        2        MULTIPOINT (0.00000 0.00000, 1.00000 1.00000)
+        3                              POINT (0.00000 0.00000)
+        dtype: geometry
+
+        >>> s.minimum_rotated_rectangle()
+        0    POLYGON ((1.00000 1.00000, 0.50000 1.50000, -0...
+        1    POLYGON ((0.00000 0.00000, 0.50000 -0.50000, 1...
+        2        LINESTRING (0.00000 0.00000, 1.00000 1.00000)
+        3                              POINT (0.00000 0.00000)
+        dtype: geometry
+
+        See also
+        --------
+        GeoSeries.envelope : bounding rectangle
+        """
+        return _delegate_geo_method("minimum_rotated_rectangle", self)
+
     @property
     def exterior(self):
         """Returns a ``GeoSeries`` of LinearRings representing the outer
@@ -2574,6 +2614,124 @@ GeometryCollection
         """
         return _binary_op("hausdorff_distance", self, other, align, densify=densify)
 
+    def frechet_distance(self, other, align=True, densify=None):
+        """Returns a ``Series`` containing the Frechet distance to aligned `other`.
+
+        The Fréchet distance is a measure of similarity: it is the greatest distance
+        between any point in A and the closest point in B. The discrete distance is an
+        approximation of this metric: only vertices are considered. The parameter
+        ``densify`` makes this approximation less coarse by splitting the line segments
+        between vertices before computing the distance.
+
+        Fréchet distance sweep continuously along their respective curves and the
+        direction of curves is significant. This makes it a better measure of similarity
+        than Hausdorff distance for curve or surface matching.
+
+        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 every 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 frechet distance of each geometry of GeoSeries
+        to a single geometry:
+
+        >>> point = Point(-1, 0)
+        >>> s.frechet_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.frechet_distance(s2, align=True)
+        0         NaN
+        1    2.121320
+        2    3.162278
+        3    2.000000
+        4         NaN
+        dtype: float64
+        >>> s.frechet_distance(s2, align=False)
+        0    0.707107
+        1    4.123106
+        2    2.000000
+        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([(0, 0), (10, 0), (0, 15)])])
+        >>> l2 = geopandas.GeoSeries([LineString([(0, 0), (20, 15), (9, 11)])])
+        >>> l1.frechet_distance(l2)
+        0    18.027756
+        dtype: float64
+        >>> l1.frechet_distance(l2, densify=0.25)
+        0    16.77051
+        dtype: float64
+        """
+        return _binary_op("frechet_distance", self, other, align, densify=densify)
+
     #
     # Binary operations that return a GeoSeries
     #
@@ -3099,6 +3257,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
     #
@@ -4136,7 +4397,7 @@ def _get_index_for_parts(orig_idx, outer_idx, ignore_index, index_parts):
             #    starting at 0 in each run
             run_start = np.r_[True, outer_idx[:-1] != outer_idx[1:]]
             counts = np.diff(np.r_[np.nonzero(run_start)[0], len(outer_idx)])
-            inner_index = (~run_start).cumsum()
+            inner_index = (~run_start).cumsum(dtype=outer_idx.dtype)
             inner_index -= np.repeat(inner_index[run_start], counts)
 
         else:
diff --git a/geopandas/tools/hilbert_curve.py b/geopandas/tools/hilbert_curve.py
index 4568dea5..4d42abac 100644
--- a/geopandas/tools/hilbert_curve.py
+++ b/geopandas/tools/hilbert_curve.py
@@ -103,6 +103,8 @@ def _continuous_to_discrete(vals, val_range, n):
 
     """
     width = val_range[1] - val_range[0]
+    if width == 0:
+        return np.zeros_like(vals, dtype=np.uint32)
     res = (vals - val_range[0]) * (n / width)
 
     np.clip(res, 0, n, out=res)


Test patch:
diff --git a/geopandas/tests/test_extension_array.py b/geopandas/tests/test_extension_array.py
index c4daa240..fe1cd106 100644
--- a/geopandas/tests/test_extension_array.py
+++ b/geopandas/tests/test_extension_array.py
@@ -131,7 +131,7 @@ def data_missing_for_sorting():
     This should be three items [B, NA, A] with
     A < B and NA missing.
     """
-    return from_shapely([Point(0, 1), None, Point(0, 0)])
+    return from_shapely([Point(1, 2), None, Point(0, 0)])
 
 
 @pytest.fixture
diff --git a/geopandas/tests/test_geom_methods.py b/geopandas/tests/test_geom_methods.py
index b7a753b7..f1f58dea 100644
--- a/geopandas/tests/test_geom_methods.py
+++ b/geopandas/tests/test_geom_methods.py
@@ -1,5 +1,4 @@
 import string
-import sys
 import warnings
 
 import numpy as np
@@ -178,6 +177,11 @@ class TestGeomMethods:
         self.g12 = GeoSeries([self.l5])
         self.g13 = GeoSeries([self.l6])
 
+        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
@@ -368,6 +372,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)])
@@ -619,6 +657,52 @@ class TestGeomMethods:
         ):
             self.g0.hausdorff_distance(self.g9)
 
+    @pytest.mark.skipif(
+        not (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="requires frechet_distance in shapely 2.0+",
+    )
+    def test_frechet_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.frechet_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.frechet_distance(self.g9, align=True)
+            )
+
+        # expected returns
+        val_1 = 1.0
+        val_2 = np.sqrt(2) / 4
+        val_3 = np.sqrt(2) / 2
+        val_4 = (np.sqrt(2) / 2) * 10
+        expected = Series(
+            np.array([val_1, val_1, val_2, val_3, val_4, np.nan, np.nan]), self.g0.index
+        )
+        assert_array_dtype_equal(
+            expected, self.g0.frechet_distance(self.g9, align=False)
+        )
+
+        expected = Series(
+            np.array([np.sqrt(100**2 + (100 - 5) ** 2)]), self.g12.index
+        )
+        assert_array_dtype_equal(
+            expected, self.g12.frechet_distance(self.g13, densify=0.25)
+        )
+
+    @pytest.mark.skipif(
+        (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
+        reason="frechet_distance not implemented for shapely<2",
+    )
+    def test_frechet_distance_not(self):
+        with pytest.raises(
+            NotImplementedError, match="shapely >= 2.0 or PyGEOS is required"
+        ):
+            self.g0.frechet_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))
@@ -1189,6 +1273,20 @@ class TestGeomMethods:
         assert isinstance(e, GeoSeries)
         assert self.g3.crs == e.crs
 
+    def test_minimum_rotated_rectangle(self):
+        s = GeoSeries([self.sq, self.t5], crs=3857)
+        r = s.minimum_rotated_rectangle()
+        exp = GeoSeries.from_wkt(
+            [
+                "POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))",
+                "POLYGON ((2 0, 2 3, 3 3, 3 0, 2 0))",
+            ]
+        )
+
+        assert np.all(r.normalize().geom_equals_exact(exp, 0.001))
+        assert isinstance(r, GeoSeries)
+        assert s.crs == r.crs
+
     @pytest.mark.skipif(
         not (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
         reason=(
@@ -1251,10 +1349,6 @@ class TestGeomMethods:
         with pytest.warns(FutureWarning, match="Currently, index_parts defaults"):
             assert_geoseries_equal(expected, s.explode())
 
-    @pytest.mark.skipif(
-        sys.platform == "win32" and sys.version_info < (3, 9),
-        reason="Inconsistent int dtype",
-    )
     @pytest.mark.parametrize("index_name", [None, "test"])
     def test_explode_geodataframe(self, index_name):
         s = GeoSeries([MultiPoint([Point(1, 2), Point(2, 3)]), Point(5, 5)])
@@ -1274,10 +1368,6 @@ class TestGeomMethods:
         expected_df = expected_df.set_index(expected_index)
         assert_frame_equal(test_df, expected_df)
 
-    @pytest.mark.skipif(
-        sys.platform == "win32" and sys.version_info < (3, 9),
-        reason="Inconsistent int dtype",
-    )
     @pytest.mark.parametrize("index_name", [None, "test"])
     def test_explode_geodataframe_level_1(self, index_name):
         # GH1393
@@ -1368,10 +1458,6 @@ class TestGeomMethods:
         exploded_df = gdf.explode(column="col1", ignore_index=True)
         assert_geodataframe_equal(exploded_df, expected_df)
 
-    @pytest.mark.skipif(
-        sys.platform == "win32" and sys.version_info < (3, 9),
-        reason="Inconsistent int dtype",
-    )
     @pytest.mark.parametrize("outer_index", [1, (1, 2), "1"])
     def test_explode_pandas_multi_index(self, outer_index):
         index = MultiIndex.from_arrays(
@@ -1479,10 +1565,6 @@ class TestGeomMethods:
         test_df = df.explode(ignore_index=True, index_parts=True)
         assert_frame_equal(test_df, expected_df)
 
-    @pytest.mark.skipif(
-        sys.platform == "win32" and sys.version_info < (3, 9),
-        reason="Inconsistent int dtype",
-    )
     def test_explode_order(self):
         df = GeoDataFrame(
             {"vals": [1, 2, 3]},
@@ -1512,10 +1594,6 @@ class TestGeomMethods:
         )
         assert_geodataframe_equal(test_df, expected_df)
 
-    @pytest.mark.skipif(
-        sys.platform == "win32" and sys.version_info < (3, 9),
-        reason="Inconsistent int dtype",
-    )
     def test_explode_order_no_multi(self):
         df = GeoDataFrame(
             {"vals": [1, 2, 3]},
@@ -1534,10 +1612,6 @@ class TestGeomMethods:
         )
         assert_geodataframe_equal(test_df, expected_df)
 
-    @pytest.mark.skipif(
-        sys.platform == "win32" and sys.version_info < (3, 9),
-        reason="Inconsistent int dtype",
-    )
     def test_explode_order_mixed(self):
         df = GeoDataFrame(
             {"vals": [1, 2, 3]},
@@ -1566,10 +1640,6 @@ class TestGeomMethods:
         )
         assert_geodataframe_equal(test_df, expected_df)
 
-    @pytest.mark.skipif(
-        sys.platform == "win32" and sys.version_info < (3, 9),
-        reason="Inconsistent int dtype",
-    )
     def test_explode_duplicated_index(self):
         df = GeoDataFrame(
             {"vals": [1, 2, 3]},
@@ -1692,10 +1762,6 @@ class TestGeomMethods:
         )
         assert_frame_equal(self.g11.get_coordinates(ignore_index=True), expected)
 
-    @pytest.mark.skipif(
-        sys.platform == "win32" and sys.version_info < (3, 9),
-        reason="Inconsistent int dtype",
-    )
     @pytest.mark.skipif(
         not (compat.USE_PYGEOS or compat.USE_SHAPELY_20),
         reason="get_coordinates not implemented for shapely<2",
diff --git a/geopandas/tests/test_pandas_methods.py b/geopandas/tests/test_pandas_methods.py
index 2f569b4d..d721f03e 100644
--- a/geopandas/tests/test_pandas_methods.py
+++ b/geopandas/tests/test_pandas_methods.py
@@ -1,4 +1,5 @@
 import os
+from packaging.version import Version
 
 import numpy as np
 from numpy.testing import assert_array_equal
@@ -249,8 +250,10 @@ def test_astype(s, df):
     assert s.astype(str)[0] == "POINT (0 0)"
 
     res = s.astype(object)
-    assert isinstance(res, pd.Series) and not isinstance(res, GeoSeries)
-    assert res.dtype == object
+    if not Version(pd.__version__) == Version("2.1.0"):
+        # https://github.com/geopandas/geopandas/issues/2948 - bug in pandas 2.1.0
+        assert isinstance(res, pd.Series) and not isinstance(res, GeoSeries)
+        assert res.dtype == object
 
     df = df.rename_geometry("geom_list")
 
@@ -433,6 +436,16 @@ def test_fillna_series(s):
     assert_geoseries_equal(res, s2)
 
 
+def test_fillna_inplace(s):
+    s2 = GeoSeries([Point(0, 0), None, Point(2, 2)])
+    arr = s2.array
+    s2.fillna(Point(1, 1), inplace=True)
+    assert_geoseries_equal(s2, s)
+    if compat.PANDAS_GE_21:
+        # starting from pandas 2.1, there is support to do this actually inplace
+        assert s2.array is arr
+
+
 def test_dropna():
     s2 = GeoSeries([Point(0, 0), None, Point(2, 2)])
     res = s2.dropna()
diff --git a/geopandas/tools/tests/test_hilbert_curve.py b/geopandas/tools/tests/test_hilbert_curve.py
index 143ee504..ca7c4fa6 100644
--- a/geopandas/tools/tests/test_hilbert_curve.py
+++ b/geopandas/tools/tests/test_hilbert_curve.py
@@ -1,3 +1,4 @@
+import numpy as np
 from shapely.geometry import Point
 from shapely.wkt import loads
 
@@ -63,3 +64,12 @@ def test_empty(geoseries_points, empty):
         ValueError, match="cannot be computed on a GeoSeries with empty"
     ):
         s.hilbert_distance()
+
+
+def test_zero_width():
+    # special case of all points on the same line -> avoid warnings because
+    # of division by 0 and introducing NaN
+    s = geopandas.GeoSeries([Point(0, 0), Point(0, 2), Point(0, 1)])
+    with np.errstate(all="raise"):
+        result = s.hilbert_distance()
+    assert np.array(result).argsort().tolist() == [0, 2, 1]