Merge pull request #29 from GeoscienceAustralia/fix/handle_am_for_dem_in_bounds

fix for checking dem intersection over the AM

Author: abradley60

dem_handler/utils/spatial.py

@@ -593,21 +593,36 @@ def check_dem_type_in_bounds(
     logger.info(f"Searching {dem_index_path}")
 
     gdf = gpd.read_file(dem_index_path, layer=layer)
-    bounding_box = shapely.geometry.box(*bounds)
+    intersecting_tiles = 0
 
-    if gdf.crs is not None:
-        # ensure same crs
-        bounding_box = (
-            gpd.GeoSeries([bounding_box], crs="EPSG:4326").to_crs(gdf.crs).iloc[0]
+    # handle antimeridian crossing
+    if bounds[0] > bounds[2]:
+        logger.warning(
+            "bounds cross the antimeridian, searching for intersections either side."
         )
+        bounds_east, bounds_west = split_bounds_at_antimeridian(bounds)
+        search_shapes = [
+            shapely.geometry.box(*bounds_east),
+            shapely.geometry.box(*bounds_west),
+        ]
     else:
-        logger.info('No crs found for index file. Assuming EPSG:4326"')
-        bounding_box = gpd.GeoSeries([bounding_box], crs="EPSG:4326")
-    # Find rows that intersect with the bounding box
-    intersecting_tiles = gdf[gdf.intersects(bounding_box)]
-    if len(intersecting_tiles) == 0:
+        search_shapes = [shapely.geometry.box(*bounds)]
+
+    for bounding_box in search_shapes:
+        if gdf.crs is not None:
+            # ensure same crs
+            bounding_box = (
+                gpd.GeoSeries([bounding_box], crs="EPSG:4326").to_crs(gdf.crs).iloc[0]
+            )
+        else:
+            logger.info('No crs found for index file. Assuming EPSG:4326"')
+            bounding_box = gpd.GeoSeries([bounding_box], crs="EPSG:4326")
+        # Count rows that intersect with the bounding box
+        intersecting_tiles += len(gdf[gdf.intersects(bounding_box)])
+
+    if intersecting_tiles == 0:
         logger.info(f"No intersecting tiles found")
         return False
     else:
-        logger.info(f"{len(intersecting_tiles)} intersecting tiles found")
+        logger.info(f"{intersecting_tiles} intersecting tiles found")
         return True

tests/utils/test_utils.py

@@ -214,6 +214,25 @@ class BoundsDEMCheckCase:
     is_error=False,
 )
 
+# crosses AM over the coast, no intersect
+test_antimeridian_no_intersect_with_rema = BoundsDEMCheckCase(
+    dem_type="rema",
+    resolution=32,
+    bounds=(173.430893, -71.618423, -178.032867, -68.765106),
+    in_bounds=False,
+    is_error=False,
+)
+
+# crosses AM over the coast, no intersect
+test_antimeridian_no_intersect_with_cop30 = BoundsDEMCheckCase(
+    dem_type="cop30",
+    resolution=30,
+    bounds=(173.430893, -71.618423, -178.032867, -68.765106),
+    in_bounds=False,
+    is_error=False,
+)
+
+
 test_cases = [
     test_invalid_dem_type,
     test_invalid_dem_resolution,
@@ -224,6 +243,8 @@ class BoundsDEMCheckCase:
     test_heard_island_bounds_with_cop30,
     test_antimeridian_with_cop30,
     test_antimeridian_with_rema,
+    test_antimeridian_no_intersect_with_rema,
+    test_antimeridian_no_intersect_with_cop30,
 ]