Repository: geopandas/geopandas
Pull request #2297 (2022-01-13 20:38:04)

ENH: add distance along the Hilbert Curve
This is a numpy-based implementation of the distance along the Hilbert Curve, similar as we have in [dask-geopandas](https://github.com/geopandas/dask-geopandas/blob/main/dask_geopandas/hilbert_distance.py) but thus not using numba but a vectorized implementation based on the version in GEOS (which is basically https://github.com/rawrunprotected/hilbert_curves)

With the version of GEOS / https://github.com/rawrunprotected/hilbert_curves, it was actually trivially to have it run vectorized on a numpy array of discrete coordinates as well (literally the same code).

I am directly creating a PR in geopandas here, and not in dask-geopandas, since this is a feature we were thinking to move to base geopandas anyway (we can use it here as well to sort geometries).

This still needs tests, docs, and hooking it up in the public API, etc.



jorisvandenbossche: Based on a first test, this new implementation gives the same result as the dask-geopandas one, while actually being faster. On 1 million points, I get ca 500ms for the numba-based implementation of dask-geopandas vs ~80ms for this numpy-based implementation. 
(this is only for the lower-level encoding of discrete coordinates to distance, so the actual difference starting from a GeoSeries will be smaller)

Notebook: https://notebooksharing.space/view/fb83908941339cc7654dbe271c02b2df49eaf8cc4d434077da58066f9cd1ff86#displayOptions=
martinfleis: I suggest using this implementation in `dask-geopandas` as well (before we do 0.1 release) as it allows us dropping `numba` as a dependency. 
jorisvandenbossche: > I suggest using this implementation in dask-geopandas as well (before we do 0.1 release) as it allows us dropping numba as a dependency.

Ah, yes, that's a good point. I think we can duplicate it on the short term, and later remove it again from dask-geopandas to rely on the geopandas version.
brendan-ward: What about adding this in Shapely instead of GeoPandas?  This seems more tightly coupled to working with geometries than to GeoPandas specifically, and may be easier to leverage anything from GEOS C API if this becomes available there and is more performant than the numpy solution.
martinfleis: @brendan-ward I think that long-term, that is the best solution. I think that the main benefit of having it here right now is that we don't have to wait for shapely to release it and then expose it here about 3 releases from now.
brendan-ward: Simplification and likely speedup: instead of calculating the `total_bounds` of the input geometries, calculate the extent (in numpy) of their midpoints calculated here.  This should also lead to better distribution of geometries along the curve for geometries that have large extents (since it would be scaling their midpoints relative to each other rather than midpoints to the total bounds of the original geometries).
jorisvandenbossche: > (side note: the codecov messages make it quite a bit harder to review code, maybe those should be turned off or handled differently during review)

Yeah, those are indeed annoying (it's because I didn't add any tests ..). Small hint that you can easily hide them all at once with:

![Screenshot from 2022-01-28 09-03-43](https://user-images.githubusercontent.com/1020496/151511646-a954f7bf-12d2-4699-a2ec-e24235287a7c.png)


jorisvandenbossche: Thanks for the review @brendan-ward !

So I am first going to update https://github.com/geopandas/dask-geopandas/pull/134 with those comments, have that merged, and then just copy it over here (as mentioned by Martin, short term I think it's good to include this code in those several places, to make it easier to start using it, agreed that long term it would fit in shapely).

> Simplification and likely speedup: instead of calculating the `total_bounds` of the input geometries, calculate the extent (in numpy) of their midpoints calculated here. This should also lead to better distribution of geometries along the curve for geometries that have large extents (since it would be scaling their midpoints relative to each other rather than midpoints to the total bounds of the original geometries).

Yes, I actually already had a commit locally (pushed it now) that calculate the "total bounds" from the bounds, which we access anyway (and the underlying implementation of `total_bounds` in pygeos/shapely is exactly that, getting the `bounds` and calculating the min/max of that, so we were effectively extracting the bounds twice). 
But indeed I can take that a step further and only do this "total bounds" calculation after we calculated the midpoints.
jorisvandenbossche: I think this was actually ready.
martinfleis: Can you add a changelog note? Otherwise this seems to be ready.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index d8facce8..1fbfb46f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,8 @@ New features and improvements:
 
 - The Parquet and Feather IO functions now support the latest 1.0.0-beta.1 version
   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).
 - 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/base.py b/geopandas/base.py
index b4c6a908..d72a7edc 100644
--- a/geopandas/base.py
+++ b/geopandas/base.py
@@ -3484,6 +3484,38 @@ GeometryCollection
             return False
         return self._data.equals(other._data)
 
+    def hilbert_distance(self, total_bounds=None, level=16):
+        """
+        Calculate the distance along a Hilbert curve.
+
+        The distances are calculated for the midpoints of the geometries in the
+        GeoDataFrame, and using the total bounds of the GeoDataFrame.
+
+        The Hilbert distance can be used to spatially sort GeoPandas
+        objects, by mapping two dimensional geometries along the Hilbert curve.
+
+        Parameters
+        ----------
+        total_bounds : 4-element array, optional
+            The spatial extent in which the curve is constructed (used to
+            rescale the geometry midpoints). By default, the total bounds
+            of the full GeoDataFrame or GeoSeries will be computed. If known,
+            you can pass the total bounds to avoid this extra computation.
+        level : int (1 - 16), default 16
+            Determines the precision of the curve (points on the curve will
+            have coordinates in the range [0, 2^level - 1]).
+
+        Returns
+        -------
+        Series
+            Series containing distance along the curve for geometry
+        """
+        from geopandas.tools.hilbert_curve import _hilbert_distance
+
+        distances = _hilbert_distance(self, total_bounds=total_bounds, level=level)
+
+        return distances
+
 
 class _CoordinateIndexer(object):
     # see docstring GeoPandasBase.cx property above
diff --git a/geopandas/tools/hilbert_curve.py b/geopandas/tools/hilbert_curve.py
new file mode 100644
index 00000000..364b0ef4
--- /dev/null
+++ b/geopandas/tools/hilbert_curve.py
@@ -0,0 +1,187 @@
+import numpy as np
+import pandas as pd
+
+
+def _hilbert_distance(gdf, total_bounds=None, level=16):
+    """
+    Calculate the distance along a Hilbert curve.
+
+    The distances are calculated for the midpoints of the geometries in the
+    GeoDataFrame.
+
+    Parameters
+    ----------
+    gdf : GeoDataFrame
+    total_bounds : 4-element array
+        Total bounds of geometries - array
+    level : int (1 - 16), default 16
+        Determines the precision of the curve (points on the curve will
+        have coordinates in the range [0, 2^level - 1]).
+
+    Returns
+    ---------
+    Pandas Series containing distances along the Hilbert curve
+
+    """
+    if gdf.is_empty.any() | gdf.geometry.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()
+
+    # 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")
+
+
+def _continuous_to_discrete_coords(bounds, level, total_bounds):
+    """
+    Calculates mid points & ranges of geoms and returns
+    as discrete coords
+
+    Parameters
+    ----------
+
+    bounds : Bounds of each geometry - array
+
+    p : The number of iterations used in constructing the Hilbert curve
+
+    total_bounds : Total bounds of geometries - array
+
+    Returns
+    ---------
+    Discrete two-dimensional numpy array
+    Two-dimensional array Array of hilbert distances for each geom
+
+    """
+    # Hilbert Side length
+    side_length = (2**level) - 1
+
+    # Calculate mid points for x and y bound coords - returns array
+    x_mids = (bounds[:, 0] + bounds[:, 2]) / 2.0
+    y_mids = (bounds[:, 1] + bounds[:, 3]) / 2.0
+
+    # Calculate x and y range of total bound coords - returns array
+    if total_bounds is None:
+        total_bounds = (
+            np.nanmin(x_mids),
+            np.nanmin(y_mids),
+            np.nanmax(x_mids),
+            np.nanmax(y_mids),
+        )
+
+    xmin, ymin, xmax, ymax = total_bounds
+
+    # Transform continuous value to discrete integer for each dimension
+    x_int = _continuous_to_discrete(x_mids, (xmin, xmax), side_length)
+    y_int = _continuous_to_discrete(y_mids, (ymin, ymax), side_length)
+
+    return x_int, y_int
+
+
+def _continuous_to_discrete(vals, val_range, n):
+    """
+    Convert a continuous one-dimensional array to discrete integer values
+    based their ranges
+
+    Parameters
+    ----------
+    vals : Array of continuous values
+
+    val_range : Tuple containing range of continuous values
+
+    n : Number of discrete values
+
+    Returns
+    ---------
+    One-dimensional array of discrete ints
+
+    """
+    width = val_range[1] - val_range[0]
+    res = (vals - val_range[0]) * (n / width)
+
+    np.clip(res, 0, n, out=res)
+    return res.astype(np.uint32)
+
+
+# Fast Hilbert curve algorithm by http://threadlocalmutex.com/
+# From C++ https://github.com/rawrunprotected/hilbert_curves
+# (public domain)
+
+
+MAX_LEVEL = 16
+
+
+def _interleave(x):
+    x = (x | (x << 8)) & 0x00FF00FF
+    x = (x | (x << 4)) & 0x0F0F0F0F
+    x = (x | (x << 2)) & 0x33333333
+    x = (x | (x << 1)) & 0x55555555
+    return x
+
+
+def _encode(level, x, y):
+
+    x = np.asarray(x, dtype="uint32")
+    y = np.asarray(y, dtype="uint32")
+
+    if level > MAX_LEVEL:
+        raise ValueError("Level out of range")
+
+    x = x << (16 - level)
+    y = y << (16 - level)
+
+    # Initial prefix scan round, prime with x and y
+    a = x ^ y
+    b = 0xFFFF ^ a
+    c = 0xFFFF ^ (x | y)
+    d = x & (y ^ 0xFFFF)
+
+    A = a | (b >> 1)
+    B = (a >> 1) ^ a
+    C = ((c >> 1) ^ (b & (d >> 1))) ^ c
+    D = ((a & (c >> 1)) ^ (d >> 1)) ^ d
+
+    a = A.copy()
+    b = B.copy()
+    c = C.copy()
+    d = D.copy()
+
+    A = (a & (a >> 2)) ^ (b & (b >> 2))
+    B = (a & (b >> 2)) ^ (b & ((a ^ b) >> 2))
+    C ^= (a & (c >> 2)) ^ (b & (d >> 2))
+    D ^= (b & (c >> 2)) ^ ((a ^ b) & (d >> 2))
+
+    a = A.copy()
+    b = B.copy()
+    c = C.copy()
+    d = D.copy()
+
+    A = (a & (a >> 4)) ^ (b & (b >> 4))
+    B = (a & (b >> 4)) ^ (b & ((a ^ b) >> 4))
+    C ^= (a & (c >> 4)) ^ (b & (d >> 4))
+    D ^= (b & (c >> 4)) ^ ((a ^ b) & (d >> 4))
+
+    # Final round and projection
+    a = A.copy()
+    b = B.copy()
+    c = C.copy()
+    d = D.copy()
+
+    C ^= (a & (c >> 8)) ^ (b & (d >> 8))
+    D ^= (b & (c >> 8)) ^ ((a ^ b) & (d >> 8))
+
+    # Undo transformation prefix scan
+    a = C ^ (C >> 1)
+    b = D ^ (D >> 1)
+
+    # Recover index bits
+    i0 = x ^ y
+    i1 = b | (0xFFFF ^ (i0 | a))
+
+    return ((_interleave(i1) << 1) | _interleave(i0)) >> (32 - 2 * level)


Test patch:
diff --git a/geopandas/tools/tests/test_hilbert_curve.py b/geopandas/tools/tests/test_hilbert_curve.py
new file mode 100644
index 00000000..14dfb072
--- /dev/null
+++ b/geopandas/tools/tests/test_hilbert_curve.py
@@ -0,0 +1,66 @@
+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
+
+
+def test_hilbert_distance():
+    # test the actual Hilbert Code algorithm against some hardcoded values
+    geoms = geopandas.GeoSeries.from_wkt(
+        [
+            "POINT (0 0)",
+            "POINT (1 1)",
+            "POINT (1 0)",
+            "POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))",
+        ]
+    )
+    result = _hilbert_distance(geoms, 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)
+    assert result.tolist() == [0, 42, 63, 10]
+
+    result = _hilbert_distance(geoms, total_bounds=(0, 0, 1, 1), level=16)
+    assert result.tolist() == [0, 2863311530, 4294967295, 715827882]
+
+
+@pytest.fixture
+def geoseries_points():
+    p1 = Point(1, 2)
+    p2 = Point(2, 3)
+    p3 = Point(3, 4)
+    p4 = Point(4, 1)
+    return geopandas.GeoSeries([p1, p2, p3, p4])
+
+
+def test_hilbert_distance_level(geoseries_points):
+    with pytest.raises(ValueError):
+        geoseries_points.hilbert_distance(level=20)
+
+
+def test_specified_total_bounds(geoseries_points):
+    result = geoseries_points.hilbert_distance(
+        total_bounds=geoseries_points.total_bounds
+    )
+    expected = geoseries_points.hilbert_distance()
+    assert_series_equal(result, expected)
+
+
+@pytest.mark.parametrize(
+    "empty",
+    [
+        None,
+        loads("POLYGON EMPTY"),
+    ],
+)
+def test_empty(geoseries_points, empty):
+    s = geoseries_points
+    s.iloc[-1] = empty
+    with pytest.raises(
+        ValueError, match="cannot be computed on a GeoSeries with empty"
+    ):
+        s.hilbert_distance()