[Backport to 3.3.x][Fixes #9070] Filter by extent issue with dateline… (

#9090) * [Backport to 3.3.x][Fixes #9070] Filter by extent issue with dateline (#9071) * [Fixes #9070] Filter by extent issue with dateline * - Remove wrong code accordingly to the review (cherry picked from commit cda0d06) # Conflicts: # geonode/base/tests.py * [CircleCI] Fix test cases
GeoNode · Apr 14, 2022 · 35c2168 · 35c2168
1 parent acb0d5d
commit 35c2168
Show file tree

Hide file tree

Showing 2 changed files with 178 additions and 25 deletions.
diff --git a/geonode/base/bbox_utils.py b/geonode/base/bbox_utils.py
@@ -17,10 +17,21 @@
 #
 #########################################################################
 
+import math
+import copy
+import json
+
 from decimal import Decimal
+from typing import Union, List, Generator
+from shapely import affinity
+from shapely.ops import split
+from shapely.geometry import (
+    mapping,
+    Polygon,
+    LineString,
+    GeometryCollection)
 
-from django.contrib.gis.geos import Polygon
-from django.db.models import Q
+from django.contrib.gis.geos import Polygon as DjangoPolygon
 
 
 class BBOXHelper:
@@ -43,7 +54,7 @@ def from_xy(cls, xy):
         return cls(xy)
 
     def as_polygon(self):
-        return Polygon.from_bbox((self.xmin, self.ymin, self.xmax, self.ymax))
+        return DjangoPolygon.from_bbox((self.xmin, self.ymin, self.xmax, self.ymax))
 
 
 def normalize_x_value(value):
@@ -57,7 +68,7 @@ def polygon_from_bbox(bbox, srid=4326):
     """
     Constructs a Polygon object with srid from a provided bbox.
     """
-    poly = Polygon.from_bbox(bbox)
+    poly = DjangoPolygon.from_bbox(bbox)
     poly.srid = srid
     return poly
 
@@ -78,27 +89,144 @@ def filter_bbox(queryset, bbox):
             _bbox_index += 1
         bboxes[_bbox_index].append(_y)
 
+    search_queryset = None
     for _bbox in bboxes:
         _bbox = list(map(Decimal, _bbox))
+        search_polygon = polygon_from_bbox((_bbox[0], _bbox[1], _bbox[2], _bbox[3]))
+        for search_polygon_dl in [DjangoPolygon.from_ewkt(_p.wkt) for _p in split_polygon(json.loads(search_polygon.json), output_format="polygons")]:
+            _qs = queryset.filter(ll_bbox_polygon__intersects=search_polygon_dl)
+            search_queryset = _qs if search_queryset is None else search_queryset | _qs
+
+    return search_queryset
+
+
+def check_crossing(lon1: float, lon2: float, validate: bool = False, dlon_threshold: float = 180.0):
+    """
+    ref.: https://towardsdatascience.com/around-the-world-in-80-lines-crossing-the-antimeridian-with-python-and-shapely-c87c9b6e1513
+    Assuming a minimum travel distance between two provided longitude coordinates,
+    checks if the 180th meridian (antimeridian) is crossed.
+    """
+
+    if validate and any([abs(x) > 180.0 for x in [lon1, lon2]]):
+        raise ValueError("longitudes must be in degrees [-180.0, 180.0]")
+    return abs(lon2 - lon1) > dlon_threshold
+
 
-        # Return all layers when the search extent exceeds 360deg
-        if abs(_bbox[0] - _bbox[2]) >= 360:
-            return queryset.all()
-
-        x_min = normalize_x_value(_bbox[0])
-        x_max = normalize_x_value(_bbox[2])
-
-        # When the search extent crosses the 180th meridian, we'll need to search
-        # on two conditions
-        if x_min > x_max:
-            left_polygon = polygon_from_bbox((-180, _bbox[1], x_max, _bbox[3]))
-            right_polygon = polygon_from_bbox((x_min, _bbox[1], 180, _bbox[3]))
-            queryset = queryset.filter(
-                Q(ll_bbox_polygon__intersects=left_polygon) |
-                Q(ll_bbox_polygon__intersects=right_polygon)
-            )
-        # Otherwise, we do a simple polygon-based search
+def translate_polygons(geometry_collection: GeometryCollection,
+                       output_format: str = "geojson") -> Generator[
+    Union[List[dict], List[Polygon]], None, None
+]:
+    """
+    ref.: https://towardsdatascience.com/around-the-world-in-80-lines-crossing-the-antimeridian-with-python-and-shapely-c87c9b6e1513
+    """
+    for polygon in geometry_collection:
+        (minx, _, maxx, _) = polygon.bounds
+        if minx < -180:
+            geo_polygon = affinity.translate(polygon, xoff=360)
+        elif maxx > 180:
+            geo_polygon = affinity.translate(polygon, xoff=-360)
+        else:
+            geo_polygon = polygon
+
+        yield_geojson = output_format == "geojson"
+        yield json.dumps(mapping(geo_polygon)) if (yield_geojson) else geo_polygon
+
+
+def split_polygon(geojson: dict, output_format: str = "geojson", validate: bool = False) -> Union[
+    List[dict], List[Polygon], GeometryCollection
+]:
+    """
+    ref.: https://towardsdatascience.com/around-the-world-in-80-lines-crossing-the-antimeridian-with-python-and-shapely-c87c9b6e1513
+    Given a GeoJSON representation of a Polygon, returns a collection of
+    'antimeridian-safe' constituent polygons split at the 180th meridian,
+    ensuring compliance with GeoJSON standards (https://tools.ietf.org/html/rfc7946#section-3.1.9)
+    Assumptions:
+      - Any two consecutive points with over 180 degrees difference in
+        longitude are assumed to cross the antimeridian
+      - The polygon spans less than 360 degrees in longitude (i.e. does not wrap around the globe)
+      - However, the polygon may cross the antimeridian on multiple occasions
+    Parameters:
+        geojson (dict): GeoJSON of input polygon to be split. For example:
+                        {
+                        "type": "Polygon",
+                        "coordinates": [
+                          [
+                            [179.0, 0.0], [-179.0, 0.0], [-179.0, 1.0],
+                            [179.0, 1.0], [179.0, 0.0]
+                          ]
+                        ]
+                        }
+        output_format (str): Available options: "geojson", "polygons", "geometrycollection"
+                             If "geometrycollection" returns a Shapely GeometryCollection.
+                             Otherwise, returns a list of either GeoJSONs or Shapely Polygons
+        validate (bool): Checks if all longitudes are within [-180.0, 180.0]
+
+    Returns:
+        List[dict]/List[Polygon]/GeometryCollection: antimeridian-safe polygon(s)
+    """
+    output_format = output_format.replace("-", "").strip().lower()
+    coords_shift = copy.deepcopy(geojson["coordinates"])
+    shell_minx = shell_maxx = None
+    split_meridians = set()
+    splitter = None
+
+    for ring_index, ring in enumerate(coords_shift):
+        if len(ring) < 1:
+            continue
         else:
-            search_polygon = polygon_from_bbox((x_min, _bbox[1], x_max, _bbox[3]))
-            queryset = queryset.filter(ll_bbox_polygon__intersects=search_polygon)
-    return queryset
+            ring_minx = ring_maxx = ring[0][0]
+            crossings = 0
+
+        for coord_index, (lon, _) in enumerate(ring[1:], start=1):
+            lon_prev = ring[coord_index - 1][0]  # [0] corresponds to longitude coordinate
+            if check_crossing(lon, lon_prev, validate=validate):
+                direction = math.copysign(1, lon - lon_prev)
+                coords_shift[ring_index][coord_index][0] = lon - (direction * 360.0)
+                crossings += 1
+
+            x_shift = coords_shift[ring_index][coord_index][0]
+            if x_shift < ring_minx:
+                ring_minx = x_shift
+            if x_shift > ring_maxx:
+                ring_maxx = x_shift
+
+        # Ensure that any holes remain contained within the (translated) outer shell
+        if (ring_index == 0):  # by GeoJSON definition, first ring is the outer shell
+            shell_minx, shell_maxx = (ring_minx, ring_maxx)
+        elif (ring_minx < shell_minx):
+            ring_shift = [[x + 360, y] for (x, y) in coords_shift[ring_index]]
+            coords_shift[ring_index] = ring_shift
+            ring_minx, ring_maxx = (x + 360 for x in (ring_minx, ring_maxx))
+        elif (ring_maxx > shell_maxx):
+            ring_shift = [[x - 360, y] for (x, y) in coords_shift[ring_index]]
+            coords_shift[ring_index] = ring_shift
+            ring_minx, ring_maxx = (x - 360 for x in (ring_minx, ring_maxx))
+
+        if crossings:  # keep track of meridians to split on
+            if ring_minx < -180:
+                split_meridians.add(-180)
+            if ring_maxx > 180:
+                split_meridians.add(180)
+
+    n_splits = len(split_meridians)
+    if n_splits > 1:
+        raise NotImplementedError(
+            """Splitting a Polygon by multiple meridians (MultiLineString)
+               not supported by Shapely"""
+        )
+    elif n_splits == 1:
+        split_lon = next(iter(split_meridians))
+        meridian = [[split_lon, -90.0], [split_lon, 90.0]]
+        splitter = LineString(meridian)
+
+    shell, *holes = coords_shift if splitter else geojson["coordinates"]
+    if splitter:
+        split_polygons = split(Polygon(shell, holes), splitter)
+    else:
+        split_polygons = GeometryCollection([Polygon(shell, holes)])
+
+    geo_polygons = list(translate_polygons(split_polygons, output_format))
+    if output_format == "geometrycollection":
+        return GeometryCollection(geo_polygons)
+    else:
+        return geo_polygons
diff --git a/geonode/base/tests.py b/geonode/base/tests.py
@@ -50,12 +50,13 @@
 )
 
 from django.conf import settings
-from django.shortcuts import reverse
+from django.contrib.gis.geos import Polygon
 from django.template import Template, Context
 from django.contrib.auth import get_user_model
 from django.core.exceptions import ObjectDoesNotExist
 from django.core.files.storage import default_storage as storage
 from django.test import Client, TestCase, override_settings, SimpleTestCase
+from django.shortcuts import reverse
 
 from geonode.base.models import CuratedThumbnail
 from geonode.base.middleware import ReadOnlyMiddleware, MaintenanceMiddleware
@@ -748,6 +749,30 @@ def test_visible_notifications(self):
         self.assertFalse(show_notification('monitoring_alert', self.user))
         self.assertTrue(show_notification('request_download_resourcebase', self.user))
 
+    def test_extent_filter_crossing_dateline(self):
+        from .bbox_utils import filter_bbox
+
+        _ll = None
+        try:
+            bbox = [166.06619, -22.40043, 172.09202, -13.03425]
+            _ll = Layer.objects.create(
+                uuid=str(uuid4()),
+                owner=self.user,
+                name='test_extent_filter_crossing_dateline',
+                title='test_extent_filter_crossing_dateline',
+                alternate='geonode:test_extent_filter_crossing_dateline',
+                is_approved=True,
+                is_published=True,
+                ll_bbox_polygon=Polygon.from_bbox(bbox)
+            )
+            self.assertListEqual(list(_ll.ll_bbox_polygon.extent), bbox, _ll.ll_bbox_polygon.extent)
+            self.assertTrue(Layer.objects.filter(title=_ll.title).exists(), Layer.objects.all())
+            _qs = filter_bbox(Layer.objects.all(), '-180.0000,-39.7790,-164.2456,9.2702,134.0552,-39.7790,180.0000,9.2702')
+            self.assertTrue(_qs.filter(title=_ll.title), Layer.objects.all() | _qs.all())
+        finally:
+            if _ll:
+                _ll.delete()
+
 
 class TestHtmlTagRemoval(SimpleTestCase):