Merge pull request #32 from GeoscienceAustralia/30-pipeline-rema-antimeridian-handling

30 pipeline rema antimeridian handling

Author: abradley60

dem_handler/dem/geoid.py

@@ -61,6 +61,7 @@ def apply_geoid(
     dem_profile: dict,
     geoid_path=str | Path,
     buffer_pixels: int = 2,
+    dem_mask_buffer: int = 0,
     save_path: str | Path = "",
     mask_array: np.ndarray | None = None,
     method: Literal["add", "subtract"] = "add",
@@ -77,6 +78,10 @@ def apply_geoid(
         Path to the geoid file
     buffer_pixels : int, optional
         Additional pixel buffer for geoid, by default 2.
+    dem_mask_buffer : int, optional
+        An additional buffer for the mask that gets applied when
+        reading the geoid. The mask is based on the dem bounds. Value
+        is in geographic units. e.g. degrees for 4326 and metres for 3031.
     save_path : str | Path, optional
         Location to save dem, by default "".
     mask_array : np.ndarray | None, optional
@@ -107,11 +112,13 @@ def apply_geoid(
         dem_crs = dem_profile["crs"].to_epsg()
         if geoid_crs != dem_crs:
             mask_dem_bounds = transform_polygon(
-                shapely.geometry.box(*dem_bounds), dem_crs, geoid_crs
+                shapely.geometry.box(*dem_bounds).buffer(dem_mask_buffer),
+                dem_crs,
+                geoid_crs,
             ).bounds
         geoid_array, geoid_transform = rasterio.mask.mask(
             src,
-            [shapely.geometry.box(*mask_dem_bounds)],
+            [shapely.geometry.box(*mask_dem_bounds).buffer(dem_mask_buffer)],
             all_touched=True,
             crop=True,
             pad=True,

dem_handler/dem/rema.py

@@ -21,7 +21,11 @@
 
 from dem_handler.dem.geoid import apply_geoid
 from dem_handler.download.aws import download_egm_08_geoid
-
+from dem_handler.utils.spatial import (
+    check_bounds_likely_cross_antimeridian,
+    split_antimeridian_shape_into_east_west_bounds,
+)
+from dem_handler.utils.raster import reproject_and_merge_rasters
 
 # Create a custom type that allows use of BoundingBox or tuple(left, bottom, right, top)
 BBox = BoundingBox | tuple[float | int, float | int, float | int, float | int]
@@ -43,6 +47,7 @@ def get_rema_dem_for_bounds(
     return_over_ocean: bool = True,
     geoid_tif_path: Path | str = "egm_08_geoid.tif",
     download_geoid: bool = False,
+    check_geoid_crosses_antimeridian: bool = True,
     num_cpus: int = 1,
     num_tasks: int | None = None,
     return_paths: bool = False,
@@ -77,6 +82,9 @@ def get_rema_dem_for_bounds(
         Path to the existing ellipsoid file, by default "egm_08_geoid.tif"
     download_geoid : bool, optional
         Flag to download the ellipsoid file, by default False
+    check_geoid_crosses_antimeridian : bool, optional
+        Check if the geoid crosses the antimeridian. Set to False if it is known
+        The requested bounds do not cross the antimeridian to stop false positives.
     num_cpus : int, optional
         Number of cpus to be used for parallel download, by default 1.
         Setting to -1 will use all available cpus
@@ -125,6 +133,26 @@ def get_rema_dem_for_bounds(
     # Log the requested bounds
     logging.info(f"Getting REMA DEM for bounds: {bounds.bounds}")
 
+    # Check if bounds cross the antimeridian
+    if bounds.left > bounds.right:
+        logging.warning(
+            f"left longitude value ({bounds[0]}) is greater than the right longitude value {({bounds[2]})} "
+            f"for the bounds provided. Assuming the bounds cross the antimeridian : {bounds}"
+        )
+        dem_crosses_antimeridian = True
+    else:
+        dem_crosses_antimeridian = False
+        # run a basic to check if the bounds likely cross the antimeridian but
+        # are just formatted wrong. If so, warn the user.
+        if bounds_src_crs == 4326:
+            if check_bounds_likely_cross_antimeridian(bounds):
+                logging.warning(
+                    "Provided bounds have very large longitude extent. If the shape crosses the "
+                    f"antimeridian, reformat the bounds as : ({bounds[2]}, {bounds[1]}, {bounds[0]}, {bounds[3]}). "
+                    "For large areas, provide the inputs bounds in 3031 to avoid transform errors between "
+                    "coordinate systems."
+                )
+
     if bounds_src_crs != REMA_CRS:
         logging.warning(
             f"Transforming bounds from {bounds_src_crs} to {REMA_CRS}. This may return data beyond the requested bounds. If this is not desired, provide the bounds in EPSG:{REMA_CRS}."
@@ -228,47 +256,114 @@ def get_rema_dem_for_bounds(
 
     if dem_novalues_count == 0 and ellipsoid_heights:
         # we have data everywhere and the values are already ellipsoid referenced
+        logging.info(f"All DEM values already referenced to ellipsoid")
+        logging.info(f"Returning DEM referenced to ellipsoidal heights")
         if save_path:
             logging.info(f"DEM saved to : {save_path}")
         logging.info(f"Dem array shape = {dem_array.shape}")
         return dem_array, dem_profile, raster_paths
     else:
-        geoid_bounds = bounds
-        if bounds_crs != GEOID_CRS:
-            geoid_bounds = transform_polygon(
-                box(*bounds.bounds), bounds_crs, GEOID_CRS
-            ).bounds
-        if not download_geoid and not geoid_tif_path.exists():
+        # we need to apply the geoid to the DEM
+        if bounds_crs == GEOID_CRS:
+            # bounds already in correct CRS (4326 for egm_08)
+            geoid_bounds = bounds
+            geoid_poly = box(*bounds.bounds)
+        else:
+            geoid_poly = transform_polygon(box(*bounds.bounds), bounds_crs, GEOID_CRS)
+            geoid_bounds = geoid_poly.bounds
+
+        # check if the bounds likely the antimeridian
+        if check_geoid_crosses_antimeridian:
+            logging.info(
+                "Checking if geoid likely crosses the antimeridian. "
+                "If this is not desired, set check_geoid_crosses_antimeridian = False"
+            )
+            geoid_crosses_antimeridian = check_bounds_likely_cross_antimeridian(
+                geoid_bounds
+            )
+        else:
+            geoid_crosses_antimeridian = False
+        if geoid_crosses_antimeridian or dem_crosses_antimeridian:
+            logging.info(
+                "Geoid crosses antimeridian, splitting into east and west hemispheres"
+            )
+            # separate the geoid into east and west sections
+            east_geoid_tif_path = geoid_tif_path.with_stem(
+                geoid_tif_path.stem + "_eastern"
+            )
+            west_geoid_tif_path = geoid_tif_path.with_stem(
+                geoid_tif_path.stem + "_western"
+            )
+            # construct the bounds for east and west hemisphere
+            east_geoid_bounds, west_geoid_bounds = (
+                split_antimeridian_shape_into_east_west_bounds(
+                    geoid_poly, buffer_degrees=0.1
+                )
+            )
+            logging.info(f"East geoid bounds : {east_geoid_bounds.bounds}")
+            logging.info(f"West geoid bounds : {west_geoid_bounds.bounds}")
+            geoid_tifs_to_apply = [east_geoid_tif_path, west_geoid_tif_path]
+            geoid_bounds_to_apply = [
+                east_geoid_bounds.bounds,
+                west_geoid_bounds.bounds,
+            ]
+        else:
+            geoid_tifs_to_apply = [geoid_tif_path]
+            geoid_bounds_to_apply = [geoid_bounds]  # bounds already tuple
+
+        if not download_geoid and not all(p.exists() for p in geoid_tifs_to_apply):
             raise FileExistsError(
-                f"Geoid file does not exist: {geoid_tif_path}. "
+                f"Required geoid files do not exist: {geoid_tifs_to_apply}. "
                 "correct path or set download_geoid = True"
             )
-        elif download_geoid and not geoid_tif_path.exists():
+        elif download_geoid and not all(p.exists() for p in geoid_tifs_to_apply):
             logging.info(f"Downloading the egm_08 geoid")
-            download_egm_08_geoid(geoid_tif_path, bounds=geoid_bounds)
-        elif download_geoid and geoid_tif_path.exists():
-            # Check that the existing geiod covers the dem
-            with rasterio.open(geoid_tif_path) as src:
-                existing_geoid_bounds = shapely.geometry.box(*src.bounds)
-            if existing_geoid_bounds.covers(shapely.geometry.box(*bounds.bounds)):
-                logging.info(
-                    f"Skipping geoid download. The existing geoid file covers the DEM bounds. Existing geoid file: {geoid_tif_path}."
-                )
-            else:
-                logging.info(
-                    f"The existing geoid file does not cover the DEM bounds. A new geoid file covering the bounds will be downloaded, overwriting the existing geiod file: {geoid_tif_path}."
-                )
-                download_egm_08_geoid(geoid_tif_path, bounds=geoid_bounds)
-
-        logging.info(f"Using geoid file: {geoid_tif_path}")
+            for geoid_path, geoid_bounds in zip(
+                geoid_tifs_to_apply, geoid_bounds_to_apply
+            ):
+                download_egm_08_geoid(geoid_path, bounds=geoid_bounds)
+        elif download_geoid and all(p.exists() for p in geoid_tifs_to_apply):
+            # Check that the existing geoid covers the dem
+            for geoid_path, geoid_bounds in zip(
+                geoid_tifs_to_apply, geoid_bounds_to_apply
+            ):
+                with rasterio.open(geoid_path) as src:
+                    existing_geoid_bounds = shapely.geometry.box(*src.bounds)
+                if existing_geoid_bounds.covers(shapely.geometry.box(*bounds.bounds)):
+                    logging.info(
+                        f"Skipping geoid download. The existing geoid file covers the DEM bounds. Existing geoid file: {geoid_path}."
+                    )
+                else:
+                    logging.info(
+                        f"The existing geoid file does not cover the DEM bounds. A new geoid file covering the bounds will be downloaded, overwriting the existing geiod file: {geoid_tif_path}."
+                    )
+                    download_egm_08_geoid(geoid_path, bounds=geoid_bounds)
+
+    if geoid_crosses_antimeridian:
+        logging.info(
+            f"Reproject and merge east and west hemisphere geoid rasters to EPSG:{REMA_CRS}"
+        )
+        reproject_and_merge_rasters(
+            geoid_tifs_to_apply, REMA_CRS, save_path=geoid_tif_path
+        )
+        # add a larger buffer to ensure geoid is applied correctly to all of dem
+        dem_mask_buffer = 5_000
+    else:
+        # no buffer is required
+        dem_mask_buffer = 0
 
     if ellipsoid_heights:
         logging.info(f"Returning DEM referenced to ellipsoidal heights")
+        logging.info(f"Applying geoid to DEM : {geoid_tif_path}")
+        # Apply the geoid only on areas where no rema data was found
+        # i.e. these values were set to zero and will be replaced with
+        # ellipsoid heights
         dem_array = apply_geoid(
             dem_array=dem_array,
             dem_profile=dem_profile,
             geoid_path=geoid_tif_path,
             buffer_pixels=2,
+            dem_mask_buffer=dem_mask_buffer,
             save_path=save_path,
             mask_array=dem_novalues_mask,
             method="add",
@@ -279,11 +374,13 @@ def get_rema_dem_for_bounds(
         # rema dem is by default referenced to the ellipsoid. We do this only in
         # areas with data, leaving the nodata areas at zero values
         logging.info(f"Returning DEM referenced to geoid heights")
+        logging.info(f"Applying geoid to DEM : {geoid_tif_path}")
         dem_array = apply_geoid(
             dem_array=dem_array,
             dem_profile=dem_profile,
             geoid_path=geoid_tif_path,
             buffer_pixels=2,
+            dem_mask_buffer=dem_mask_buffer,
             save_path=save_path,
             mask_array=dem_values_mask,
             method="subtract",

dem_handler/utils/raster.py

@@ -483,3 +483,51 @@ def read_raster_with_bounds(file_path, bounds, buffer_pixels=0):
         )
 
     return data, profile
+
+
+def reproject_and_merge_rasters(
+    file_paths: list[str] | list[Path], target_crs, save_path=None
+) -> tuple[np.ndarray, dict]:
+    """
+    Reproject multiple raster files to a common CRS and merge them into a single array.
+
+    Parameters
+    ----------
+    file_paths : list of str or Path
+        List of file paths to the input raster files.
+    target_crs : rasterio CRS or compatible
+        The target coordinate reference system to which all rasters will be reprojected.
+    save_path : str or Path, optional
+        Path to save the merged raster to disk. If None, the merged raster is not saved.
+
+    Returns
+    -------
+    merged_array : np.ndarray
+        The merged raster data as a NumPy array.
+    merged_profile : dict
+        The rasterio profile (metadata) associated with the merged raster.
+
+    Notes
+    -----
+    - Each raster in `file_paths` is first reprojected individually using `reproject_raster`.
+    - The resulting arrays are then merged using `merge_arrays_with_geometadata` with a
+      pixel-wise maximum as the default method.
+    - The output CRS of the merged array matches `target_crs`.
+    """
+
+    arrays = []
+    profiles = []
+
+    for raster_path in file_paths:
+        array, profile = reproject_raster(raster_path, target_crs)
+        arrays.append(array)
+        profiles.append(profile)
+
+    array, profile = merge_arrays_with_geometadata(
+        arrays=arrays,
+        profiles=profiles,
+        method="max",
+        output_path=save_path,
+    )
+
+    return array, profile

dem_handler/utils/spatial.py

@@ -3,7 +3,7 @@
 import geopandas as gpd
 import pyproj
 from shapely import segmentize
-from shapely.geometry import Polygon, box
+from shapely.geometry import Polygon, box, MultiPolygon
 from pyproj.database import query_utm_crs_info
 from pyproj.aoi import AreaOfInterest
 from pyproj import CRS
@@ -349,7 +349,7 @@ def split_bounds_at_antimeridian(
     return (bounds_eastern_hemisphere, bounds_western_hemisphere)
 
 
-def get_all_lat_lon_coords(
+def get_all_lon_lat_coords(
     geom: shapely.geometry.Polygon | shapely.geometry.MultiPolygon,
 ) -> tuple[list, list]:
     """
@@ -626,3 +626,67 @@ def check_dem_type_in_bounds(
     else:
         logger.info(f"{intersecting_tiles} intersecting tiles found")
         return True
+
+
+def split_antimeridian_shape_into_east_west_bounds(
+    antimeridian_shape: Polygon | MultiPolygon, buffer_degrees: float = 0
+) -> tuple[BoundingBox, BoundingBox]:
+    """
+    Split a geometry that crosses the antimeridian into eastern and western bounding boxes.
+
+    This function derives longitude/latitude bounds from a polygon or multipolygon
+    that spans the antimeridian (±180°) and returns two bounding boxes:
+    one covering the eastern hemisphere (-180 to <0) and one covering the
+    western hemisphere (>0 to 180).
+
+    An optional buffer can be applied to slightly expand the bounds, which is
+    useful when downloading or resampling rasters to avoid edge effects.
+
+    Parameters
+    ----------
+    antimeridian_shape : Polygon or MultiPolygon
+        Geometry that crosses the antimeridian, expressed in a geographic
+        coordinate system (longitude, latitude).
+    buffer_degrees : float, optional
+        Buffer (in degrees) to apply to the resulting bounds in all directions.
+        Longitude buffers expand away from the antimeridian, and latitude buffers
+        are clipped to valid ranges [-90, 90]. Default is 0 (no buffer).
+
+    Returns
+    -------
+    east_bounds : BoundingBox
+        Bounding box covering the eastern hemisphere portion
+        (-180, min_lat, east_lon_min, max_lat).
+    west_bounds : BoundingBox
+        Bounding box covering the western hemisphere portion
+        (west_lon_min, min_lat, 180, max_lat).
+
+    Notes
+    -----
+    - The split is determined by the maximum negative longitude (east side)
+      and the minimum positive longitude (west side) present in the geometry.
+    - This is commonly used when downloading or applying global rasters
+      (e.g. geoids) that cannot be queried across the antimeridian in a single request.
+    """
+    # get the lat and lon points to construct the bounds either side of the AM
+    lons, lats = get_all_lon_lat_coords(antimeridian_shape)
+    east_lon_min = max([x for x in lons if x < 0])
+    west_lon_min = min([x for x in lons if x > 0])
+    min_lat, max_lat = min(lats), max(lats)
+
+    # add slight buffer to the bounds
+    if buffer_degrees:
+        logging.info(
+            f"Adding small buffer of {buffer_degrees} degrees to bounds either side of the antimeridian"
+        )
+        east_lon_min = (
+            east_lon_min + buffer_degrees
+        )  # push further east (-177 -> -176.9)
+        west_lon_min = west_lon_min - buffer_degrees  # push further west (177 -> 176.9)
+        min_lat = max(min_lat - buffer_degrees, -90)  # ensure valid
+        max_lat = min(max_lat + buffer_degrees, 90)
+
+    # construct the bounds
+    east_bounds = BoundingBox(-180, min_lat, east_lon_min, max_lat)
+    west_bounds = BoundingBox(west_lon_min, min_lat, 180, max_lat)
+    return east_bounds, west_bounds