Handle anti-meridian issues during spatial indexing

Still doesn't handle possible long-line issues - assumes straight lines
are still straight lines even when reprojected, so assumes the corners
of bounding boxes are sufficient for reprojecting a bounding box.

kart/spatial_tree.py

@@ -462,41 +462,122 @@ def get_envelope_for_indexing(geom, transforms, feature_oid):
     It is always true that s <= n. Normally w <= e unless it crosses the anti-meridian, in which case e < w.
     """
 
-    # Result initialised to an infinite and negative area which disappears when unioned.
     result = None
 
     try:
-        raw_envelope = geom.envelope(only_2d=True, calculate_if_missing=True)
+        minmax_envelope = _transpose_gpkg_envelope(
+            geom.envelope(only_2d=True, calculate_if_missing=True)
+        )
 
         for transform in transforms:
-            envelope = _transform_envelope(raw_envelope, transform)
+            envelope = transform_minmax_envelope(minmax_envelope, transform)
             result = union_of_envelopes(result, envelope)
         return result
     except Exception as e:
         L.warning(f"Couldn't index feature {feature_oid}:\n{e}")
         return None
 
 
-def _transform_envelope(raw_envelope, transform):
+def _transpose_gpkg_envelope(envelope):
+    """
+    GPKG uses the envelope format (min-x, max-x, min-y, max-y). We use the envelope format (w, s, e, n).
+    We transpose GPKG envelope to (min-x, min-y, max-x, max-y), so that it least it has the same axis-order as our
+    format, and we handle anti-meridian issues seperately (see transform_minmax_envelope).
+    """
+    return envelope[0], envelope[2], envelope[1], envelope[3]
+
+
+def transform_minmax_envelope(minmax_envelope, transform):
+    """
+    Given an envelope in (min-x, min-y, max-x, max-y) format in any CRS, transforms it to EPSG:4326 using the given
+    transform, then returns an axis-aligned envelope in EPSG:4326 in (w, s, e, n) order that bounds the original
+    (but which may have a slightly larger area due to the axis-aligned edges not lining up with the original).
+    The returned envelope has w <= e unless it crosses the antimeridian, in which case e < w.
+    """
+
     # FIXME: This function is not 100% correct - it should segmentize each edge of the envelope and transform
     # all the segments, to account for the fact that these edges should sometimes curve when reprojected.
-    # But right now we just reproject the corners and draw lines between them - in rare cases the result
+    # But right now we just reproject the corners and assume straight lines between them - in rare cases the result
     # won't cover the entire geometry.
-    corners = ogr.Geometry(ogr.wkbLineString)
-    # At least we transform all 4 corners, which avoids some errors caused by only transforming 2 of them.
-
-    # The raw_envelope that we get from GPKG / OGR has the following format: min-x, max-x, min-y, max-y.
-    corners.AddPoint_2D(raw_envelope[0], raw_envelope[2])
-    corners.AddPoint_2D(raw_envelope[0], raw_envelope[3])
-    corners.AddPoint_2D(raw_envelope[1], raw_envelope[2])
-    corners.AddPoint_2D(raw_envelope[1], raw_envelope[3])
-    corners.Transform(transform)
-    values = [corners.GetPoint_2D(c) for c in range(4)]
-    x_values = [v[0] for v in values]
-    y_values = [v[1] for v in values]
-    # FIXME: This also assumes the geometry doesn't cross the anti-meridian.
-    # Return w, s, e, n:
-    return min(x_values), min(y_values), max(x_values), max(y_values)
+
+    ring = anticlockwise_ring_from_minmax_envelope(minmax_envelope)
+    ring.Transform(transform)
+    ring_points = [ring.GetPoint_2D(i) for i in range(5)]
+
+    x_values = [p[0] for p in ring_points]
+    y_values = [p[1] for p in ring_points]
+
+    wrapped_min_x = _wrap_lon(min(x_values))
+    min_y = min(y_values)
+    wrapped_max_x = _wrap_lon(max(x_values))
+    max_y = max(y_values)
+
+    if wrapped_min_x > wrapped_max_x:
+        # This happens if the transform was a no-op and leaves the points so that they obviously cross the antimeridian,
+        # eg has x-values 170 to 190. In this case we can just return the wrapped min and max as e and w, since e > w
+        # is how we indicate that the envelope crosses the meridian.
+        return (wrapped_min_x, min_y, wrapped_max_x, max_y)
+
+    if _is_anticlockwise(ring_points):
+        # The ring is still anticlockwise as we constructed it, so, it has not been inverted by reprojecting it into
+        # the longitude range [-180, 180]. So we can return it as is.
+        return (wrapped_min_x, min_y, wrapped_max_x, max_y)
+
+    # The ring was anticlockwise, but when projected and EPSG:4326 into the range [-180, 180] it became clockwise.
+    # We need to try different interprerations of the ring until we find one where it is anticlockwise (this will
+    # cross the meridian). Once we've found this interpretation, we can treat the min-x and max-x as w and e.
+    ring_points = _fix_ring_winding_order(ring_points)
+    x_values = [p[0] for p in ring_points]
+
+    wrapped_min_x = _wrap_lon(min(x_values))
+    wrapped_max_x = _wrap_lon(max(x_values))
+    return (wrapped_min_x, min_y, wrapped_max_x, max_y)
+
+
+def anticlockwise_ring_from_minmax_envelope(minmax_envelope):
+    """Given an envelope in (min-x, min-y, max-x, max-y) format, builds an anticlockwise ring around it."""
+    ring = ogr.Geometry(ogr.wkbLinearRing)
+    # The envelope has the following format: min-x, min-y, max-x, max-ys.
+    # We start at min-x, min-y and travel around it in an anti-clockwise direction:
+    ring.AddPoint_2D(minmax_envelope[0], minmax_envelope[1])
+    ring.AddPoint_2D(minmax_envelope[2], minmax_envelope[1])
+    ring.AddPoint_2D(minmax_envelope[2], minmax_envelope[3])
+    ring.AddPoint_2D(minmax_envelope[0], minmax_envelope[3])
+    ring.AddPoint_2D(minmax_envelope[0], minmax_envelope[1])
+    return ring
+
+
+def _is_anticlockwise(ring):
+    """
+    Does a polygon area calculation to determine whether the given ring or polygon is clockwise.
+    For explanation see https://en.wikipedia.org/wiki/Shoelace_formula
+    Ring should be a list of (x, y) tuples where the first and last point are the same.
+    """
+    result = 0
+    for i in range(len(ring) - 1):
+        result += ring[i][0] * ring[i + 1][1] - ring[i][1] * ring[i + 1][0]
+    return result >= 0
+
+
+def _fix_ring_winding_order(points):
+    """
+    Shifts each point in turn eastwards by 360 degrees around the globe until the winding order is fixed.
+    This works on any geometry, but has O(n^2) efficiency, so ideally should only be used on small geometries.
+    """
+    if _is_anticlockwise(points):
+        return points
+
+    # Make a copy of points, and make each point a list not a tuple so it can be modified easily.
+    points = [list(p) for p in points]
+
+    sorted_x_values = sorted(set(p[0] for p in points))
+    for x in sorted_x_values:
+        for p in points:
+            if p[0] <= x:
+                p[0] += 360
+        if _is_anticlockwise(points):
+            return points
+    raise AssertionError("This should never happen")
 
 
 def _unwrap_lon_envelope(w, e):

tests/test_spatial_tree.py

@@ -1,8 +1,16 @@
 from dataclasses import dataclass
 import pytest
 
+from osgeo import osr
+
+from kart.crs_util import make_crs
 from kart.sqlalchemy.sqlite import sqlite_engine
-from kart.spatial_tree import EnvelopeEncoder, union_of_envelopes
+from kart.spatial_tree import (
+    EnvelopeEncoder,
+    anticlockwise_ring_from_minmax_envelope,
+    transform_minmax_envelope,
+    union_of_envelopes,
+)
 from sqlalchemy.orm import sessionmaker
 
 H = pytest.helpers.helpers()
@@ -190,6 +198,65 @@ def _get_index_summary(repo_path):
         )
 
 
+# This transform leaves every point exactly where it is, even points past the antimeridian eg (185, 0)
+# We need to test this since its a special case - no other transform will result in longitudes outside
+# the range [-180, 180].
+EPSG_4326 = make_crs("EPSG:4326")
+IDENTITY_TRANSFORM = osr.CoordinateTransformation(EPSG_4326, EPSG_4326)
+
+# This transform leaves points more-or-less where they are, but it does ensure the output longitude
+# values are -180 <= x <= 180, which is what most transforms will do, so this is helps test that
+# we handle anti-meridians properly in the general case.
+NZGD2000 = make_crs("EPSG:2193")
+NZGD2000_TRANSFORM = osr.CoordinateTransformation(NZGD2000, EPSG_4326)
+
+
+@pytest.mark.parametrize(
+    "minmax_envelope,transform,expected_result",
+    [
+        ((1, 2, 3, 4), IDENTITY_TRANSFORM, (1, 2, 3, 4)),
+        ((177, -10, 184, 10), IDENTITY_TRANSFORM, (177, -10, -176, 10)),
+        ((185, 10, 190, 20), IDENTITY_TRANSFORM, (-175, 10, -170, 20)),
+        ((-190, -20, -185, -10), IDENTITY_TRANSFORM, (170, -20, 175, -10)),
+        (
+            (1347679, 5456907, 2021025, 6117225),
+            NZGD2000_TRANSFORM,
+            (170, -41, 178, -35),
+        ),
+        (
+            (1347679, 5456907, 2532792, 5740667),
+            NZGD2000_TRANSFORM,
+            (170, -41, -176, -38),
+        ),
+        (
+            (2367133, 5308558, 2513805, 5517072),
+            NZGD2000_TRANSFORM,
+            (-178, -42, -176, -40),
+        ),
+    ],
+)
+def test_transform_raw_envelope(transform, minmax_envelope, expected_result):
+    # This first part of the test just tests our assumptions about transforms:
+    # we need to be sure that the IDENTITY_TRANSFORM works differently to the other
+    # transforms, so that we can be sure the code handles both possibilities properly.
+    ring = anticlockwise_ring_from_minmax_envelope(minmax_envelope)
+    ring.Transform(transform)
+    x_values = [ring.GetPoint_2D(i)[0] for i in range(5)]
+
+    if transform == IDENTITY_TRANSFORM:
+        # IDENTITY_TRANSFORM leaves x-values as they are.
+        # This means output values can be > 180 or < -180.
+        assert set(x_values) == set([minmax_envelope[0], minmax_envelope[2]])
+    else:
+        # Any other transform will output x-values within this range,
+        # which means detecting anti-meridian crossings works a little differently.
+        assert all(-180 <= x <= 180 for x in x_values)
+
+    # This is the actual test that the transform works.
+    actual_result = transform_minmax_envelope(minmax_envelope, transform)
+    assert actual_result == pytest.approx(expected_result, abs=1e-5)
+
+
 @pytest.mark.parametrize(
     "env1,env2,expected_result",
     [