diff --git a/dem_handler/dem/geoid.py b/dem_handler/dem/geoid.py index 75dfcf9..96964e2 100644 --- a/dem_handler/dem/geoid.py +++ b/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, diff --git a/dem_handler/dem/rema.py b/dem_handler/dem/rema.py index 3c9fe82..6cc5e8e 100644 --- a/dem_handler/dem/rema.py +++ b/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", diff --git a/dem_handler/utils/raster.py b/dem_handler/utils/raster.py index 3fe1b55..3f889c4 100644 --- a/dem_handler/utils/raster.py +++ b/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 diff --git a/dem_handler/utils/spatial.py b/dem_handler/utils/spatial.py index f29d7df..62a2bac 100644 --- a/dem_handler/utils/spatial.py +++ b/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 diff --git a/tests/rema/test_dem_rema.py b/tests/rema/test_dem_rema.py index cc1a56e..e79b36e 100644 --- a/tests/rema/test_dem_rema.py +++ b/tests/rema/test_dem_rema.py @@ -87,12 +87,44 @@ class TestDem: True, ) +# antimeridian and no tile intersections +antimeridian_no_intersection_bbox = BoundingBox(179.6, -71, -179.6, -70.6) +test_antimeridian_no_intersection_ellipsoid_h = TestDem( + antimeridian_no_intersection_bbox, + os.path.join( + TEST_DATA_PATH, "rema_32m_antimeridian_no_intersection_ellipsoid_h.tif" + ), + 32, + None, + os.path.join( + GEOID_DATA_PATH, + "egm_08_geoid_rema_32m_antimeridian_no_intersection.tif", + ), + True, +) + +# antimeridian two tile intersect +antimeridian_two_tile_bbox = BoundingBox(179.5, -78.5, -179.5, -78.1) +test_antimeridian_two_tiles_ellipsoid_h = TestDem( + antimeridian_two_tile_bbox, + os.path.join(TEST_DATA_PATH, "rema_32m_antimeridian_two_tile_ellipsoid_h.tif"), + 32, + None, + os.path.join( + GEOID_DATA_PATH, + "egm_08_geoid_rema_32m_antimeridian_two_tile.tif", + ), + True, +) + test_dems_ellipsoid = [ test_single_tile_ellipsoid_h, test_four_tiles_ellipsoid_h, test_one_tile_ocean_ellipsoid_h, test_no_intersection_ellipsoid_h, test_one_tile_and_no_tile_overlap_ellipsoid_h, + test_antimeridian_no_intersection_ellipsoid_h, + test_antimeridian_two_tiles_ellipsoid_h, ] # make the geoid test set @@ -124,6 +156,21 @@ class TestDem: ), ellipsoid_heights=False, ) +test_antimeridian_no_intersection_geoid_h = replace( + test_antimeridian_no_intersection_ellipsoid_h, + dem_file=test_antimeridian_no_intersection_ellipsoid_h.dem_file.replace( + "ellipsoid", "geoid" + ), + ellipsoid_heights=False, +) +test_antimeridian_two_tiles_geoid_h = replace( + test_antimeridian_two_tiles_ellipsoid_h, + dem_file=test_antimeridian_two_tiles_ellipsoid_h.dem_file.replace( + "ellipsoid", "geoid" + ), + ellipsoid_heights=False, +) + test_dems_geoid = [ test_single_tile_geoid_h, @@ -131,6 +178,8 @@ class TestDem: test_one_tile_ocean_geoid_h, test_no_intersection_geoid_h, test_one_tile_and_no_tile_overlap_geoid_h, + test_antimeridian_no_intersection_geoid_h, + test_antimeridian_two_tiles_geoid_h, ]