Merge pull request #24 from GeoscienceAustralia/feature/check_dem_bounds_utility

adding utility function for checking bounds in a dem

Author: abradley60

dem_handler/__init__.py

@@ -1 +1,14 @@
 from dem_handler._version import __version__
+from pathlib import Path
+import typing
+
+DATA_DIR = Path(__file__).parent / Path("data")
+REMA_GPKG_PATH = DATA_DIR / Path("REMA_Mosaic_Index_v2.gpkg")
+COP30_GPKG_PATH = DATA_DIR / Path("copdem_tindex_filename.gpkg")
+
+REMAResolutions = typing.Literal[2, 10, 32]
+COPResolutions = typing.Literal[30]
+ValidDEMResolutions = typing.Literal[REMAResolutions, COPResolutions]
+
+REMA_VALID_RESOLUTIONS = typing.get_args(REMAResolutions)
+COP_VALID_RESOLUTIONS = typing.get_args(COPResolutions)

dem_handler/dem/cop_glo30.py

@@ -35,8 +35,7 @@
 # Create a custom type that allows use of BoundingBox or tuple(xmin, ymin, xmax, ymax)
 BBox = BoundingBox | tuple[float | int, float | int, float | int, float | int]
 
-DATA_DIR = Path(__file__).parents[1] / Path("data")
-COP30_GPKG_PATH = DATA_DIR / Path("copdem_tindex_filename.gpkg")
+from dem_handler import COP30_GPKG_PATH
 
 
 @log_timing

dem_handler/dem/rema.py

@@ -24,13 +24,7 @@
 # Create a custom type that allows use of BoundingBox or tuple(xmin, ymin, xmax, ymax)
 BBox = BoundingBox | tuple[float | int, float | int, float | int, float | int]
 
-DATA_DIR = Path(__file__).parents[1] / Path("data")
-REMA_GPKG_PATH = DATA_DIR / Path("REMA_Mosaic_Index_v2.gpkg")
-REMA_VALID_RESOLUTIONS = [
-    2,
-    10,
-    32,
-]  # [2, 10, 32, 100, 500, 1000] It seems there are no higher resolutions in the new index
+from dem_handler import REMA_GPKG_PATH, REMA_VALID_RESOLUTIONS, REMAResolutions
 
 
 @log_timing
@@ -39,7 +33,7 @@ def get_rema_dem_for_bounds(
     save_path: Path | str = "",
     rema_index_path: Path | str = REMA_GPKG_PATH,
     local_dem_dir: Path | str | None = None,
-    resolution: int = 2,
+    resolution: REMAResolutions = 2,
     bounds_src_crs: int = 3031,
     buffer_metres: int = 0,
     buffer_pixels: int = 0,

dem_handler/utils/spatial.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 from dataclasses import dataclass
+import geopandas as gpd
 import pyproj
 from shapely import segmentize
 from shapely.geometry import Polygon, box
@@ -20,6 +21,14 @@
 
 logger = logging.getLogger(__name__)
 
+from dem_handler import (
+    REMA_GPKG_PATH,
+    COP30_GPKG_PATH,
+    REMA_VALID_RESOLUTIONS,
+    ValidDEMResolutions,
+    COP_VALID_RESOLUTIONS,
+)
+
 
 # Construct a dataclass for bounding boxes
 @dataclass
@@ -509,3 +518,72 @@ def crop_datasets_to_bounds(
             # )
 
     return dem_array, dem_profile
+
+
+def check_dem_type_in_bounds(
+    dem_type: str, resolution: ValidDEMResolutions, bounds: BBox
+) -> bool:
+    """Check if the specified dem has data within the provided bounds. The provided dem_type is matched to either the
+    Copernicus Global 30m DEM or REMA DEM (currently implemented options). True is returned if the provided bounds
+    intersect with any tiles of the specified DEM.
+
+    Parameters
+    ----------
+    dem_type : str
+        dem type to check. Can be variations of REMA and COP. e.g. REMA_32, REMA_10,
+        cop30m, cop_30m, cop_glo30.
+    resolution: ValidDEMResolutions
+        resolution of the dem. Required to read the correct layer of the GPKG.
+    bounds : BBox
+        Bounds to check if data exists
+
+    Returns
+    -------
+    bool
+        True if the bounds intersects a tile, False otherwise.
+
+    Raises
+    -------
+    ValueError
+        If the provided dem_type cannot be matched to either the Copernicus 30m global DEM or the REMA dem.
+    """
+
+    if isinstance(bounds, BoundingBox):
+        bounds = bounds.bounds
+
+    dem_type_match = dem_type.upper()
+    if "COP" in dem_type_match and resolution in COP_VALID_RESOLUTIONS:
+        dem_index_path = COP30_GPKG_PATH
+        dem_type_formal = "Copernicus 30m global DEM"
+        layer = None  # only one layer
+    elif "REMA" in dem_type_match and resolution in REMA_VALID_RESOLUTIONS:
+        dem_index_path = REMA_GPKG_PATH
+        dem_type_formal = "REMA DEM"
+        layer = f"REMA_Mosaic_Index_v2_{resolution}m"
+    else:
+        raise ValueError(
+            f"DEM type `{dem_type}` and resolution {resolution} could not be matched to either the Copernicus 30m global DEM or a valid REMA DEM with resolutions of {REMA_VALID_RESOLUTIONS}"
+        )
+
+    logger.info(f"Checking if bounds intersect with tiles of the {dem_type_formal}")
+    logger.info(f"Searching {COP30_GPKG_PATH}")
+
+    gdf = gpd.read_file(dem_index_path, layer=layer)
+    bounding_box = shapely.geometry.box(*bounds)
+
+    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")
+    # Find rows that intersect with the bounding box
+    intersecting_tiles = gdf[gdf.intersects(bounding_box)]
+    if len(intersecting_tiles) == 0:
+        logger.info(f"No intersecting tiles found")
+        return False
+    else:
+        logger.info(f"{len(intersecting_tiles)} intersecting tiles found")
+        return True

tests/test_utils.py

@@ -0,0 +1,172 @@
+import pytest
+from dataclasses import dataclass
+from shapely.geometry import Polygon, MultiPolygon
+from dem_handler.utils.spatial import (
+    BoundingBox,
+    get_correct_bounds_from_shape_at_antimeridian,
+    check_s1_bounds_cross_antimeridian,
+    check_dem_type_in_bounds,
+)
+
+# Separate lists for clarity
+poly1 = Polygon(
+    [
+        (178.576126, -71.618423),
+        (-178.032867, -70.167343),
+        (176.938004, -68.765106),
+        (173.430893, -70.119957),
+        (178.576126, -71.618423),
+    ],
+)
+poly2 = Polygon(
+    [
+        (179.5, -10.0),  # Just west of the antimeridian
+        (-179.8, 5.0),  # Just east of the antimeridian
+        (-178.5, 15.0),  # Further east
+        (178.0, 12.0),  # West again
+        (179.5, -10.0),  # Close the ring
+    ]
+)
+poly3 = Polygon(
+    [
+        (170.0, -55.0),  # Western point
+        (-175.0, -50.0),  # Just east of antimeridian
+        (-160.0, -45.0),  # Farther east
+        (175.0, -40.0),  # Back west
+        (170.0, -55.0),  # Close the loop
+    ]
+)
+
+poly4 = MultiPolygon([poly1, poly2])
+
+shapes = [poly1, poly2, poly3, poly4]
+
+expected_bounds = [
+    BoundingBox(-178.032867, -71.618423, 173.430893, -68.765106),
+    BoundingBox(-178.5, -10, 178.0, 15),
+    BoundingBox(-160, -55, 170, -40),
+    BoundingBox(-178.032867, -71.618423, 173.430893, 15),
+]
+
+# Combine for parameterization
+test_cases = list(zip(shapes, expected_bounds))
+
+
+@pytest.mark.parametrize("shape, expected_bound", test_cases)
+def test_get_correct_bounds_from_shape_at_antimeridian(shape, expected_bound):
+    # assert the incorrect bounds still cross the AM
+    assert check_s1_bounds_cross_antimeridian(shape.bounds)
+    result = get_correct_bounds_from_shape_at_antimeridian(shape)
+    assert result == expected_bound
+
+
+@dataclass
+class BoundsDEMCheckCase:
+    dem_type: str
+    resolution: int
+    bounds: tuple[float, float, float, float]
+    in_bounds: bool
+    is_error: bool
+
+
+test_invalid_dem_type = BoundsDEMCheckCase(
+    dem_type="spaghetti",
+    resolution=90,
+    bounds=(161.96252, -70.75924, 162.10388, -70.72293),
+    in_bounds=False,
+    is_error=True,
+)
+
+test_invalid_dem_resolution = BoundsDEMCheckCase(
+    dem_type="REMA",
+    resolution=15,
+    bounds=(161.96252, -70.75924, 162.10388, -70.72293),
+    in_bounds=False,
+    is_error=True,
+)
+
+test_antarctic_bounds_with_rema = BoundsDEMCheckCase(
+    dem_type="REMA",
+    resolution=32,
+    bounds=(161.96252, -70.75924, 162.10388, -70.72293),
+    in_bounds=True,
+    is_error=False,
+)
+
+test_antarctic_bounds_with_cop30 = BoundsDEMCheckCase(
+    dem_type="copernicus",
+    resolution=30,
+    bounds=(161.96252, -70.75924, 162.10388, -70.72293),
+    in_bounds=True,
+    is_error=False,
+)
+
+test_australia_bounds_with_rema = BoundsDEMCheckCase(
+    dem_type="REMA",
+    resolution=10,
+    bounds=(133.0, -25.0, 134.0, -24.0),
+    in_bounds=False,
+    is_error=False,
+)
+
+test_heard_island_bounds_with_rema = BoundsDEMCheckCase(
+    dem_type="rema_2",
+    resolution=32,
+    bounds=(73.2144, -53.224, 73.87, -52.962),
+    in_bounds=False,
+    is_error=False,
+)
+
+test_heard_island_bounds_with_cop30 = BoundsDEMCheckCase(
+    dem_type="cop_glo30",
+    resolution=30,
+    bounds=(73.2144, -53.224, 73.87, -52.962),
+    in_bounds=True,
+    is_error=False,
+)
+
+test_antimeridian_with_cop30 = BoundsDEMCheckCase(
+    dem_type="cop_glo30",
+    resolution=30,
+    bounds=(-178.032867, -80.618423, 173.430893, -68.765106),
+    in_bounds=True,
+    is_error=False,
+)
+
+test_antimeridian_with_rema = BoundsDEMCheckCase(
+    dem_type="rema",
+    resolution=2,
+    bounds=(-178.032867, -80.618423, 173.430893, -68.765106),
+    in_bounds=True,
+    is_error=False,
+)
+
+test_cases = [
+    test_invalid_dem_type,
+    test_invalid_dem_resolution,
+    test_antarctic_bounds_with_rema,
+    test_antarctic_bounds_with_cop30,
+    test_australia_bounds_with_rema,
+    test_heard_island_bounds_with_rema,
+    test_heard_island_bounds_with_cop30,
+    test_antimeridian_with_cop30,
+    test_antimeridian_with_rema,
+]
+
+
+@pytest.mark.parametrize("test_case", test_cases)
+def test_check_dem_type_in_bounds(test_case: BoundsDEMCheckCase):
+
+    if test_case.is_error:
+        with pytest.raises(ValueError):
+            check_dem_type_in_bounds(
+                test_case.dem_type, test_case.resolution, test_case.bounds
+            )
+
+    else:
+        assert (
+            check_dem_type_in_bounds(
+                test_case.dem_type, test_case.resolution, test_case.bounds
+            )
+            == test_case.in_bounds
+        )