Repository: geopandas/geopandas

Issue:
ENH: implement `GeoSeries.sort_values` based on Hilbert and/or Morton distance
#### Is your feature request related to a problem?

We can implement `GeoSeries.sort_values` utilising Hilbert and/or Morton distance we have just merged in dask-geopandas. That would mirror [what PostGIS does](https://blog.crunchydata.com/blog/waiting-for-postgis-3-hilbert-geometry-sorting) in `ORDER BY geom`.

#### Describe the solution you'd like

I am not sure if we should just make a decision and implement sorting based on one of those algorithm or add a keyword to control which one is used. pandas already has a lot of [keywords](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.sort_values.html) so I am not sure we want to add more.

Hilbert curve is nicer than Morton as it preserves space a bit better but our implementation uses numba now. Morton is a bit faster, can be done with numpy but the result is not so ideal.

#### API breaking implications

None.




Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2950b2a0..3c288652 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,8 @@ New features and improvements:
   of the GeoParquet specification (geoparquet.org) (#2663).
 - New ``hilbert_distance()`` method that calculates the distance along a Hilbert curve
   for each geometry in a GeoSeries/GeoDataFrame (#2297).
+- Support for sorting geometries (for example, using ``sort_values()``) based on
+  the distance along the Hilbert curve (#2070).
 - Added ``minimum_bounding_circle()`` method from shapely to GeoSeries/GeoDataframe (#2621).
 - Support specifying ``min_zoom`` and ``max_zoom`` inside the ``map_kwds`` argument for ``.explore()`` (#2599).
 
diff --git a/geopandas/array.py b/geopandas/array.py
index 705caf04..41f75f7a 100644
--- a/geopandas/array.py
+++ b/geopandas/array.py
@@ -1264,7 +1264,55 @@ class GeometryArray(ExtensionArray):
         ExtensionArray.argsort
         """
         # Note: this is used in `ExtensionArray.argsort`.
-        raise TypeError("geometries are not orderable")
+        from geopandas.tools.hilbert_curve import _hilbert_distance
+
+        if self.size == 0:
+            # TODO _hilbert_distance fails for empty array
+            return np.array([], dtype="uint32")
+
+        mask_empty = self.is_empty
+        has_empty = mask_empty.any()
+        mask = self.isna() | mask_empty
+        if mask.any():
+            # if there are missing or empty geometries, we fill those with
+            # a dummy geometry so that the _hilbert_distance function can
+            # process those. The missing values are handled separately by
+            # pandas regardless of the values we return here (to sort
+            # first/last depending on 'na_position'), the distances for the
+            # empty geometries are substitued below with an appropriate value
+            geoms = self.copy()
+            indices = np.nonzero(~mask)[0]
+            if indices.size:
+                geom = self[indices[0]]
+            else:
+                # for all-empty/NA, just take random geometry
+                geom = shapely.geometry.Point(0, 0)
+
+            geoms[mask] = geom
+        else:
+            geoms = self
+        if has_empty:
+            # in case we have empty geometries, we need to expand the total
+            # bounds with a small percentage, so the empties can be
+            # deterministically sorted first
+            total_bounds = geoms.total_bounds
+            xoff = (total_bounds[2] - total_bounds[0]) * 0.01
+            yoff = (total_bounds[3] - total_bounds[1]) * 0.01
+            total_bounds += np.array([-xoff, -yoff, xoff, yoff])
+        else:
+            total_bounds = None
+        distances = _hilbert_distance(geoms, total_bounds=total_bounds)
+        if has_empty:
+            # empty geometries are sorted first ("smallest"), so fill in
+            # smallest possible value for uints
+            distances[mask_empty] = 0
+        return distances
+
+    def argmin(self, skipna: bool = True) -> int:
+        raise TypeError("geometries have no minimum or maximum")
+
+    def argmax(self, skipna: bool = True) -> int:
+        raise TypeError("geometries have no minimum or maximum")
 
     def _formatter(self, boxed=False):
         """Formatting function for scalar values.
diff --git a/geopandas/base.py b/geopandas/base.py
index f908f957..abe3c857 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -3618,9 +3618,11 @@ GeometryCollection
         """
         from geopandas.tools.hilbert_curve import _hilbert_distance
 
-        distances = _hilbert_distance(self, total_bounds=total_bounds, level=level)
+        distances = _hilbert_distance(
+            self.geometry.values, total_bounds=total_bounds, level=level
+        )
 
-        return distances
+        return pd.Series(distances, index=self.index, name="hilbert_distance")
 
 
 def _get_index_for_parts(orig_idx, outer_idx, ignore_index, index_parts):
diff --git a/geopandas/tools/hilbert_curve.py b/geopandas/tools/hilbert_curve.py
index 42189a41..69c6c423 100644
--- a/geopandas/tools/hilbert_curve.py
+++ b/geopandas/tools/hilbert_curve.py
@@ -1,8 +1,7 @@
 import numpy as np
-import pandas as pd
 
 
-def _hilbert_distance(gdf, total_bounds=None, level=16):
+def _hilbert_distance(geoms, total_bounds=None, level=16):
     """
     Calculate the distance along a Hilbert curve.
 
@@ -11,7 +10,7 @@ def _hilbert_distance(gdf, total_bounds=None, level=16):
 
     Parameters
     ----------
-    gdf : GeoDataFrame
+    geoms : GeometryArray
     total_bounds : 4-element array
         Total bounds of geometries - array
     level : int (1 - 16), default 16
@@ -20,23 +19,24 @@ def _hilbert_distance(gdf, total_bounds=None, level=16):
 
     Returns
     ---------
-    Pandas Series containing distances along the Hilbert curve
+    np.ndarray
+        Array containing distances along the Hilbert curve
 
     """
-    if gdf.is_empty.any() | gdf.geometry.isna().any():
+    if geoms.is_empty.any() | geoms.isna().any():
         raise ValueError(
             "Hilbert distance cannot be computed on a GeoSeries with empty or "
             "missing geometries.",
         )
     # Calculate bounds as numpy array
-    bounds = gdf.bounds.to_numpy()
+    bounds = geoms.bounds
 
     # Calculate discrete coords based on total bounds and bounds
     x, y = _continuous_to_discrete_coords(bounds, level, total_bounds)
     # Compute distance along hilbert curve
     distances = _encode(level, x, y)
 
-    return pd.Series(distances, index=gdf.index, name="hilbert_distance")
+    return distances
 
 
 def _continuous_to_discrete_coords(bounds, level, total_bounds):


Test patch:
diff --git a/geopandas/tests/test_extension_array.py b/geopandas/tests/test_extension_array.py
index 30a5a932..50914058 100644
--- a/geopandas/tests/test_extension_array.py
+++ b/geopandas/tests/test_extension_array.py
@@ -21,9 +21,10 @@ import pandas as pd
 from pandas.tests.extension import base as extension_tests
 
 import shapely.geometry
+from shapely.geometry import Point
 
 from geopandas.array import GeometryArray, GeometryDtype, from_shapely
-from geopandas._compat import ignore_shapely2_warnings
+from geopandas._compat import ignore_shapely2_warnings, SHAPELY_GE_20, PANDAS_GE_15
 
 import pytest
 
@@ -33,7 +34,10 @@ import pytest
 
 
 not_yet_implemented = pytest.mark.skip(reason="Not yet implemented")
-no_sorting = pytest.mark.skip(reason="Sorting not supported")
+no_minmax = pytest.mark.skip(reason="Min/max not supported")
+requires_shapely2 = pytest.mark.skipif(
+    not SHAPELY_GE_20, reason="Requires hashable geometries"
+)
 
 
 # -----------------------------------------------------------------------------
@@ -116,7 +120,7 @@ def data_for_sorting():
     This should be three items [B, C, A] with
     A < B < C
     """
-    raise NotImplementedError
+    return from_shapely([Point(0, 1), Point(1, 1), Point(0, 0)])
 
 
 @pytest.fixture
@@ -126,7 +130,7 @@ def data_missing_for_sorting():
     This should be three items [B, NA, A] with
     A < B and NA missing.
     """
-    raise NotImplementedError
+    return from_shapely([Point(0, 1), None, Point(0, 0)])
 
 
 @pytest.fixture
@@ -282,6 +286,15 @@ def all_compare_operators(request):
     return request.param
 
 
+@pytest.fixture(params=[None, lambda x: x])
+def sort_by_key(request):
+    """
+    Simple fixture for testing keys in sorting methods.
+    Tests None (no key) and the identity key.
+    """
+    return request.param
+
+
 # -----------------------------------------------------------------------------
 # Inherited tests
 # -----------------------------------------------------------------------------
@@ -458,60 +471,25 @@ class TestComparisonOps(extension_tests.BaseComparisonOpsTests):
 
 
 class TestMethods(extension_tests.BaseMethodsTests):
-    @no_sorting
+    @pytest.mark.skipif(
+        not PANDAS_GE_15, reason="sorting index not yet working with older pandas"
+    )
     @pytest.mark.parametrize("dropna", [True, False])
     def test_value_counts(self, all_data, dropna):
         pass
 
-    @no_sorting
+    @pytest.mark.skipif(
+        not PANDAS_GE_15, reason="sorting index not yet working with older pandas"
+    )
     def test_value_counts_with_normalize(self, data):
         pass
 
-    @no_sorting
-    def test_argsort(self, data_for_sorting):
-        result = pd.Series(data_for_sorting).argsort()
-        expected = pd.Series(np.array([2, 0, 1], dtype=np.int64))
-        self.assert_series_equal(result, expected)
-
-    @no_sorting
-    def test_argsort_missing(self, data_missing_for_sorting):
-        result = pd.Series(data_missing_for_sorting).argsort()
-        expected = pd.Series(np.array([1, -1, 0], dtype=np.int64))
-        self.assert_series_equal(result, expected)
-
-    @no_sorting
-    @pytest.mark.parametrize("ascending", [True, False])
-    def test_sort_values(self, data_for_sorting, ascending):
-        ser = pd.Series(data_for_sorting)
-        result = ser.sort_values(ascending=ascending)
-        expected = ser.iloc[[2, 0, 1]]
-        if not ascending:
-            expected = expected[::-1]
-
-        self.assert_series_equal(result, expected)
-
-    @no_sorting
-    @pytest.mark.parametrize("ascending", [True, False])
-    def test_sort_values_missing(self, data_missing_for_sorting, ascending):
-        ser = pd.Series(data_missing_for_sorting)
-        result = ser.sort_values(ascending=ascending)
-        if ascending:
-            expected = ser.iloc[[2, 0, 1]]
-        else:
-            expected = ser.iloc[[0, 2, 1]]
-        self.assert_series_equal(result, expected)
-
-    @no_sorting
+    @requires_shapely2
     @pytest.mark.parametrize("ascending", [True, False])
     def test_sort_values_frame(self, data_for_sorting, ascending):
-        df = pd.DataFrame({"A": [1, 2, 1], "B": data_for_sorting})
-        result = df.sort_values(["A", "B"])
-        expected = pd.DataFrame(
-            {"A": [1, 1, 2], "B": data_for_sorting.take([2, 0, 1])}, index=[2, 0, 1]
-        )
-        self.assert_frame_equal(result, expected)
-
-    @no_sorting
+        super().test_sort_values_frame(data_for_sorting, ascending)
+
+    @pytest.mark.skip(reason="searchsorted not supported")
     def test_searchsorted(self, data_for_sorting, as_series):
         pass
 
@@ -529,31 +507,23 @@ class TestMethods(extension_tests.BaseMethodsTests):
         with pytest.raises(ValueError, match=msg):
             data_missing.fillna(data_missing.take([1]))
 
-    @no_sorting
-    def test_nargsort(self):
-        pass
-
-    @no_sorting
-    def test_argsort_missing_array(self):
-        pass
-
-    @no_sorting
+    @no_minmax
     def test_argmin_argmax(self):
         pass
 
-    @no_sorting
+    @no_minmax
     def test_argmin_argmax_empty_array(self):
         pass
 
-    @no_sorting
+    @no_minmax
     def test_argmin_argmax_all_na(self):
         pass
 
-    @no_sorting
+    @no_minmax
     def test_argreduce_series(self):
         pass
 
-    @no_sorting
+    @no_minmax
     def test_argmax_argmin_no_skipna_notimplemented(self):
         pass
 
@@ -563,16 +533,16 @@ class TestCasting(extension_tests.BaseCastingTests):
 
 
 class TestGroupby(extension_tests.BaseGroupbyTests):
-    @no_sorting
+    @requires_shapely2
     @pytest.mark.parametrize("as_index", [True, False])
     def test_groupby_extension_agg(self, as_index, data_for_grouping):
-        pass
+        super().test_groupby_extension_agg(as_index, data_for_grouping)
 
-    @no_sorting
+    @requires_shapely2
     def test_groupby_extension_transform(self, data_for_grouping):
-        pass
+        super().test_groupby_extension_transform(data_for_grouping)
 
-    @no_sorting
+    @requires_shapely2
     @pytest.mark.parametrize(
         "op",
         [
@@ -584,7 +554,7 @@ class TestGroupby(extension_tests.BaseGroupbyTests):
         ids=["scalar", "list", "series", "object"],
     )
     def test_groupby_extension_apply(self, data_for_grouping, op):
-        pass
+        super().test_groupby_extension_apply(data_for_grouping, op)
 
 
 class TestPrinting(extension_tests.BasePrintingTests):
diff --git a/geopandas/tests/test_pandas_methods.py b/geopandas/tests/test_pandas_methods.py
index 5f0a7180..0373b551 100644
--- a/geopandas/tests/test_pandas_methods.py
+++ b/geopandas/tests/test_pandas_methods.py
@@ -453,6 +453,49 @@ def test_any_all():
 # Groupby / algos
 
 
+def test_sort_values():
+    s = GeoSeries([Point(0, 0), Point(2, 2), Point(0, 2)])
+    res = s.sort_values()
+    assert res.index.tolist() == [0, 2, 1]
+    res2 = s.sort_values(ascending=False)
+    assert res2.index.tolist() == [1, 2, 0]
+
+    # empty geoseries
+    assert_geoseries_equal(s.iloc[:0].sort_values(), s.iloc[:0])
+
+
+def test_sort_values_empty_missing():
+    s = GeoSeries([Point(0, 0), None, Point(), Point(1, 1)])
+    # default: NA sorts last, empty first
+    res = s.sort_values()
+    assert res.index.tolist() == [2, 0, 3, 1]
+
+    # descending: NA sorts last, empty last
+    res = s.sort_values(ascending=False)
+    assert res.index.tolist() == [3, 0, 2, 1]
+
+    # NAs first, empty first after NAs
+    res = s.sort_values(na_position="first")
+    assert res.index.tolist() == [1, 2, 0, 3]
+
+    # NAs first, descending with empyt last
+    res = s.sort_values(ascending=False, na_position="first")
+    assert res.index.tolist() == [1, 3, 0, 2]
+
+    # all missing / empty
+    s = GeoSeries([None, None, None])
+    res = s.sort_values()
+    assert res.index.tolist() == [0, 1, 2]
+
+    s = GeoSeries([Point(), Point(), Point()])
+    res = s.sort_values()
+    assert res.index.tolist() == [0, 1, 2]
+
+    s = GeoSeries([Point(), None, Point()])
+    res = s.sort_values()
+    assert res.index.tolist() == [0, 2, 1]
+
+
 def test_unique():
     s = GeoSeries([Point(0, 0), Point(0, 0), Point(2, 2)])
     exp = from_shapely([Point(0, 0), Point(2, 2)])
diff --git a/geopandas/tools/tests/test_hilbert_curve.py b/geopandas/tools/tests/test_hilbert_curve.py
index 14dfb072..143ee504 100644
--- a/geopandas/tools/tests/test_hilbert_curve.py
+++ b/geopandas/tools/tests/test_hilbert_curve.py
@@ -2,7 +2,6 @@ from shapely.geometry import Point
 from shapely.wkt import loads
 
 import geopandas
-from geopandas.tools.hilbert_curve import _hilbert_distance
 
 import pytest
 from pandas.testing import assert_series_equal
@@ -18,13 +17,13 @@ def test_hilbert_distance():
             "POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))",
         ]
     )
-    result = _hilbert_distance(geoms, total_bounds=(0, 0, 1, 1), level=2)
+    result = geoms.hilbert_distance(total_bounds=(0, 0, 1, 1), level=2)
     assert result.tolist() == [0, 10, 15, 2]
 
-    result = _hilbert_distance(geoms, total_bounds=(0, 0, 1, 1), level=3)
+    result = geoms.hilbert_distance(total_bounds=(0, 0, 1, 1), level=3)
     assert result.tolist() == [0, 42, 63, 10]
 
-    result = _hilbert_distance(geoms, total_bounds=(0, 0, 1, 1), level=16)
+    result = geoms.hilbert_distance(total_bounds=(0, 0, 1, 1), level=16)
     assert result.tolist() == [0, 2863311530, 4294967295, 715827882]