diff --git a/README.md b/README.md index 0949465..dea92e4 100644 --- a/README.md +++ b/README.md @@ -97,6 +97,40 @@ get_cop30_dem_for_bounds( ``` +### Antimeridian + +Requesting data across the antimeridian requires correctly formatted bounds. The antimeridian is a special case where the `left` value of the bounds is greater than the `right`, e.g left = 178 (eastern hemisphere) and right = -179 (western hemisphere). If bounds such as these are passed in, the dem-handler logic will return a mosaiced DEM that crossess the antimeridian projected in the most appropriate crs. For example, 3031 for high latitudes. + +```python +# specify bounds over the antimeridian +# bounds tuple must be (left, bottom, right, top) +bounds = (178, -66, -179, -60) + +get_cop30_dem_for_bounds( + bounds = bounds, + ..., +) +``` + +If a set of bounds look like they *may* cross the antimeridian, but are incorrectly formatted, a warning will be raised, but the process will run as normal. + +```python +# A valid set of bounds that almost entirely wraps the earth, and may therefore +# be a misformatted set of bounds over the antimeridian +bounds = (-179, -66, 178, -60) + +# the functions will run as intended but a warning will be raised +get_cop30_dem_for_bounds( + bounds = bounds, + ..., +) + +>>>WARNING:dem_handler.dem.cop_glo30:Provided bounds have very large longitude extent. If the shape crosses the antimeridian, reformat the bounds as : (178, -66, -179, -60) +``` + + + + ## Install Currently, we only support installing `dem-handler` from source. diff --git a/dem_handler/dem/cop_glo30.py b/dem_handler/dem/cop_glo30.py index cce6358..3faa391 100644 --- a/dem_handler/dem/cop_glo30.py +++ b/dem_handler/dem/cop_glo30.py @@ -14,9 +14,9 @@ from dem_handler.utils.spatial import ( BoundingBox, - check_s1_bounds_cross_antimeridian, + check_bounds_likely_cross_antimeridian, get_target_antimeridian_projection, - split_s1_bounds_at_am_crossing, + split_bounds_at_antimeridian, adjust_bounds_at_high_lat, crop_datasets_to_bounds, ) @@ -32,7 +32,7 @@ logger = logging.getLogger(__name__) -# Create a custom type that allows use of BoundingBox or tuple(xmin, ymin, xmax, ymax) +# 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] from dem_handler import COP30_GPKG_PATH @@ -56,6 +56,93 @@ def get_cop30_dem_for_bounds( return_paths: bool = False, download_dir: Path | None = None, ): + """ + Retrieve and mosaic COPDEM GLO-30 (Copernicus 30 m DEM) tiles covering a specified geographic bounding box. + + This function locates and optionally downloads the Copernicus GLO-30 Digital Elevation Model (DEM) + tiles that intersect the requested bounding box. It handles high-latitude adjustments, buffering, + ellipsoidal height conversion using a geoid model, and special cases where the bounds cross + the antimeridian (±180° longitude). + + The function can return either a merged DEM array with metadata or a list of intersecting DEM tile paths. + + Parameters + ---------- + bounds : BBox or tuple + The geographic bounding box of interest, either as a `BBox` object or a tuple + `(left, bottom, right, top)` in degrees. + save_path : Path + File path where the output DEM (GeoTIFF) will be saved. + ellipsoid_heights : bool, optional + If True, converts DEM heights from the geoid reference to ellipsoidal heights + using the EGM08 geoid model. Default is True. + adjust_at_high_lat : bool, optional + If True, expands the bounds near the poles to ensure adequate DEM coverage. + Default is False. + buffer_pixels : int or None, optional + Optional pixel buffer applied around the requested bounds to include a margin of DEM data. + Default is None. + buffer_degrees : int or float or None, optional + Optional geographic buffer in degrees around the bounds. Default is None. + cop30_index_path : Path, optional + Path to the COPDEM GLO-30 index GeoPackage (`.gpkg`) used to locate intersecting DEM tiles. + Default is `COP30_GPKG_PATH`. + cop30_folder_path : Path, optional + Directory containing the COPDEM GLO-30 tiles. Default is the current directory. + geoid_tif_path : Path, optional + Path to the local geoid model GeoTIFF (e.g., `egm_08_geoid.tif`) used for height conversion. + Default is `"egm_08_geoid.tif"`. + download_dem_tiles : bool, optional + If True, automatically downloads any missing DEM tiles required to cover the requested bounds. + Default is False. + download_geoid : bool, optional + If True, downloads the EGM08 geoid model if it does not exist locally. + Default is False. + num_cpus : int, optional + Number of CPU cores to use for parallel tasks such as downloading or merging tiles. + Default is 1. + num_tasks : int or None, optional + Number of parallel tasks to execute when searching or downloading tiles. + Default is None (serial execution). + return_paths : bool, optional + If True, returns only a list of file paths to intersecting DEM tiles rather than reading or merging them. + Default is False. + download_dir : Path or None, optional + Directory where downloaded DEM tiles or geoid files should be saved. Default is None (current directory). + + Returns + ------- + tuple or list + If `return_paths=True`, returns: + list of Path + File paths to the intersecting COPDEM GLO-30 tiles. + Otherwise, returns: + tuple (dem_array, dem_profile, dem_paths) + - dem_array : numpy.ndarray + The merged DEM raster data covering the requested bounds. + - dem_profile : dict + Raster metadata/profile dictionary compatible with `rasterio`. + - dem_paths : list of Path + The DEM tile paths used to produce the merged output. + + Raises + ------ + FileExistsError + If the geoid model file does not exist locally and `download_geoid=False`. + ValueError + If the input bounds are invalid or cannot be processed. + RuntimeError + If DEM merging or reprojection fails during processing. + + Notes + ----- + - If the bounds cross the antimeridian the function recursively processes each + side of the antimeridian, reprojects them into a common coordinate reference system, + and merges them into a continuous DEM. + - The DEM heights are geoid-referenced by default (EGM08 model). Set `ellipsoid_heights=True` + to obtain ellipsoidal heights (WGS84). + - At high latitudes or near the poles, `adjust_at_high_lat=True` can help include complete DEM coverage. + """ # Convert bounding box to built-in bounding box type if isinstance(bounds, tuple): @@ -65,9 +152,21 @@ def get_cop30_dem_for_bounds( logger.info(f"Getting cop30m dem for bounds: {bounds.bounds}") # Check if bounds cross the antimeridian - antimeridian_crossing = check_s1_bounds_cross_antimeridian( - bounds, max_scene_width=15 - ) + if bounds.left > bounds.right: + logger.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}" + ) + antimeridian_crossing = True + else: + antimeridian_crossing = False + # run a basic to check if the bounds likely cross the antimeridian but + # are just formatted wrong. If so, warn the user. + if check_bounds_likely_cross_antimeridian(bounds): + logger.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]})" + ) if antimeridian_crossing: logger.warning( @@ -77,7 +176,7 @@ def get_cop30_dem_for_bounds( target_crs = get_target_antimeridian_projection(bounds) logger.info(f"Splitting bounds into left and right side of antimeridian") - bounds_eastern, bounds_western = split_s1_bounds_at_am_crossing(bounds) + bounds_eastern, bounds_western = split_bounds_at_antimeridian(bounds) # Use recursion to process each side of the AM. The function is rerun # This time, antimeridian_crossing will be False enabling each side to be @@ -168,7 +267,7 @@ def get_cop30_dem_for_bounds( output_path=save_path, ) - return dem_array, dem_profile, eastern_output[2] + western_output[2] + return dem_array[0], dem_profile, eastern_output[2] + western_output[2] else: # Adjust bounds at high latitude if requested @@ -202,7 +301,7 @@ def get_cop30_dem_for_bounds( "The Cop30 DEM bounds do not fully cover the requested bounds. " "Try increasing the 'buffer_pixels' value. Note at the antimeridian " "This is expected, with bounds being slightly smaller on +ve side. " - "e.g. max_lon is 179.9999 < 180." + "e.g. right is 179.9999 < 180." ) logging.warning(warn_msg) @@ -391,7 +490,7 @@ def buffer_bounds_cop_glo30( Parameters ---------- bounds : BoundingBox | tuple[float | int, float | int, float | int, float | int] - The set of bounds (min_lon, min_lat, max_lon, max_lat) + The set of bounds (left, bottom, right, top) pixel_buffer : int | None, optional Number of pixels to buffer, by default None degree_buffer : float | int | None, optional @@ -423,12 +522,12 @@ def buffer_bounds_cop_glo30( if degree_buffer: buffer = (degree_buffer, degree_buffer) - new_xmin = max(bounds.xmin - buffer[0], -180) - new_ymin = max(bounds.ymin - buffer[1], -90) - new_xmax = min(bounds.xmax + buffer[0], 180) - new_ymax = min(bounds.ymax + buffer[1], 90) + new_left = max(bounds.left - buffer[0], -180) + new_bottom = max(bounds.bottom - buffer[1], -90) + new_right = min(bounds.right + buffer[0], 180) + new_top = min(bounds.top + buffer[1], 90) - return BoundingBox(new_xmin, new_ymin, new_xmax, new_ymax) + return BoundingBox(new_left, new_bottom, new_right, new_top) def get_cop_glo30_spacing( @@ -439,7 +538,7 @@ def get_cop_glo30_spacing( Parameters ---------- bounds : BoundingBox | tuple[float | int, float | int, float | int, float | int] - The set of bounds (min_lon, min_lat, max_lon, max_lat) + The set of bounds (left, bottom, right, top) Returns ------- @@ -455,7 +554,7 @@ def get_cop_glo30_spacing( if isinstance(bounds, tuple): bounds = BoundingBox(*bounds) - mean_latitude = abs((bounds.ymin + bounds.ymax) / 2) + mean_latitude = abs((bounds.bottom + bounds.top) / 2) minimum_pixel_spacing = 0.0002777777777777778 @@ -531,7 +630,7 @@ def make_empty_cop_glo30_profile_for_bounds( Parameters ---------- bounds : BoundingBox | tuple[float | int, float | int, float | int, float | int] - The set of bounds (min_lon, min_lat, max_lon, max_lat) + The set of bounds (left, bottom, right, top) pixel_buffer | int The number of pixels to add as a buffer to the profile @@ -552,7 +651,7 @@ def make_empty_cop_glo30_profile_for_bounds( spacing_lon, spacing_lat = get_cop_glo30_spacing(bounds) glo30_transform = get_cop_glo30_tile_transform( - bounds.xmin, bounds.ymax, spacing_lon, spacing_lat + bounds.left, bounds.top, spacing_lon, spacing_lat ) # Expand the bounds to the edges of pixels diff --git a/dem_handler/dem/geoid.py b/dem_handler/dem/geoid.py index ff670e2..75dfcf9 100644 --- a/dem_handler/dem/geoid.py +++ b/dem_handler/dem/geoid.py @@ -28,7 +28,7 @@ def read_geoid( geoid_path : str | Path Path to the GEOID file bounds : tuple - the set of bounds (min_lon, min_lat, max_lon, max_lat) + the set of bounds (left, bottom, right, top) buffer_pixels : int, optional additional pixels to buffern around bounds, by default 0 diff --git a/dem_handler/dem/rema.py b/dem_handler/dem/rema.py index 49910f0..8aeedcb 100644 --- a/dem_handler/dem/rema.py +++ b/dem_handler/dem/rema.py @@ -21,7 +21,7 @@ from dem_handler.download.aws import download_egm_08_geoid -# Create a custom type that allows use of BoundingBox or tuple(xmin, ymin, xmax, ymax) +# 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] from dem_handler import REMA_GPKG_PATH, REMA_VALID_RESOLUTIONS, REMAResolutions diff --git a/dem_handler/utils/spatial.py b/dem_handler/utils/spatial.py index 0bbc11a..49153ae 100644 --- a/dem_handler/utils/spatial.py +++ b/dem_handler/utils/spatial.py @@ -33,36 +33,72 @@ # Construct a dataclass for bounding boxes @dataclass class BoundingBox: - xmin: float | int - ymin: float | int - xmax: float | int - ymax: float | int + left: float | int + bottom: float | int + right: float | int + top: float | int @property def bounds(self) -> tuple[float | int, float | int, float | int, float | int]: - return (self.xmin, self.ymin, self.xmax, self.ymax) + return (self.left, self.bottom, self.right, self.top) @property def top_left(self) -> tuple[float | int, float | int]: - return (self.xmin, self.ymax) + return (self.left, self.top) @property def bottom_right(self) -> tuple[float | int, float | int]: - return (self.xmax, self.ymin) + return (self.right, self.bottom) + + def __getitem__(self, index: int) -> float | int: + """Allow list-style indexing for bounding box values.""" + try: + return self.bounds[index] + except IndexError: + raise IndexError("BoundingBox index out of range (valid indices: 0–3)") # Run checks on the bounding box values def __post_init__(self): - if self.ymin >= self.ymax: + if self.bottom >= self.top: raise ValueError( - "The bounding box's ymin value is greater than or equal to ymax value. Check ordering" + "The bounding box's bottom value is greater than or equal to the top value. Check ordering" ) - if self.xmin >= self.xmax: - raise ValueError( - "The bounding box's x_min value is greater than or equal to x_max value. Check ordering" + if self.left >= self.right: + logger.warning( + "The bounding box's left value is greater than or equal to the right value. " + "Assuming the bounds cross the antimeridian. Refactor the bounds if this is not correct." ) + def _check_valid(self, lon_tol: float = 0, lat_tol: float = 0): + """Check that the bounds are valid. + + Ensures: + - left and right are within [-180, 180] + - bottom and top are within [-90, 90] + + Parameters: + ------- + lon_tol: float + a longitude tolerance if the bounds can be slightly outside of [-180,180]. + lon_tol: float + a latitude tolerance if the bounds can be slightly outside of [-90,90]. + + Returns + ------- + bool + True if all bounds are within valid ranges, False otherwise. + """ + return all( + [ + -180 - lon_tol <= self.left <= 180 + lon_tol, + -180 - lon_tol <= self.right <= 180 + lon_tol, + -90 - lat_tol <= self.bottom <= 90 + lat_tol, + -90 - lat_tol <= self.top <= 90 + lat_tol, + ] + ) + -# Create a custom type that allows use of BoundingBox or tuple(xmin, ymin, xmax, ymax) +# 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] @@ -96,7 +132,8 @@ def adjust_bounds( ref_crs: int, segment_length: float = 0.1, ) -> tuple: - """_summary_ + """Adjust the bounds by considering the area covered in the ref_crs. The bounds + will generally be larger due to the rotation that occurs in reprojections. Parameters ---------- @@ -136,7 +173,7 @@ def get_local_utm( Parameters ---------- bounds : BoundingBox | tuple[float | int, float | int, float | int, float | int] - The set of bounds (min_lon, min_lat, max_lon, max_lat) + The set of bounds (left, bottom, top, right) antimeridian : bool, optional Whether the bounds cross the antimeridian, by default False @@ -149,8 +186,8 @@ def get_local_utm( bounds = BoundingBox(*bounds) logger.info("Finding best crs for area") - centre_lat = (bounds.ymin + bounds.ymax) / 2 - centre_lon = (bounds.xmin + bounds.xmax) / 2 + centre_lat = (bounds.bottom + bounds.top) / 2 + centre_lon = (bounds.left + bounds.right) / 2 if antimeridian: # force the lon to be next to antimeridian on the side with the scene centre. # e.g. (-177 + 178)/2 = 1, this is > 0 more data on -'ve side @@ -169,27 +206,27 @@ def get_local_utm( return int(crs) -def bounds_to_geojson(bounds, save_path=""): +def bounds_to_polygon(bounds, save_path=""): """ Convert a bounding box to a GeoJSON Polygon. Parameters: - bounds (tuple): A tuple of (min_lon, min_lat, max_lon, max_lat). + bounds (tuple): A tuple of (left, bottom, right, top). save_path (str): path to save the geojson Returns: dict: A GeoJSON FeatureCollection with a single Polygon feature. """ - min_lon, min_lat, max_lon, max_lat = bounds + left, bottom, right, top = bounds # Define the polygon coordinates coordinates = [ [ - [min_lon, min_lat], - [min_lon, max_lat], - [max_lon, max_lat], - [max_lon, min_lat], - [min_lon, min_lat], # Closing the polygon + [left, bottom], + [left, top], + [right, top], + [right, bottom], + [left, bottom], # Closing the polygon ] ] @@ -212,36 +249,41 @@ def bounds_to_geojson(bounds, save_path=""): return geojson -def check_s1_bounds_cross_antimeridian(bounds: BBox, max_scene_width: int = 20) -> bool: - """Check if the s1 scene bounds cross the antimeridian. The bounds of a sentinel-1 - are valid at the antimeridian, just very large. By setting a max scene width, we - can determine if the antimeridian is crossed. Alternate scenario is a bounds - with a very large width (i.e. close to the width of the earth). +def check_bounds_likely_cross_antimeridian( + bounds: BBox, max_crossing_width: int = 20 +) -> bool: + """Check if the bounds likely cross the antimeridian and may not be correctly + formatted. Can occur if the bounds are provided as (xmin, ymin, xmax, ymax) + rather than (left, bottom, right, top). For example, the bounds (-178, -60, 179, -57) + are valid, but stretch 357 degrees across the earth. It is therefore likely + the scene crosses the antimeridian with correct bounds (179, -60, -178, -57) representing + the left, bottom, right and top values respectively. The function works based on + allowable antimeridian crossing width (max_crossing_width) Parameters ---------- bounds : BoundingBox - the set of bounds (xmin, ymin, xmax, ymax) - max_scene_width : int, optional - maximum allowable width of the scene bounds in degrees, by default 20 + the set of bounds (xmin, ymin, xmax, ymax) / (left, bottom, right, top) + max_crossing_width : int, optional + maximum allowable width of bounds degrees, by default 20 Returns ------- bool - if the bounds cross the antimeridian + if the bounds likely cross the antimeridian """ if isinstance(bounds, tuple): bounds = BoundingBox(*bounds) antimeridian_xmin = -180 - bounding_xmin = antimeridian_xmin + max_scene_width # -160 by default + bounding_xmin = antimeridian_xmin + max_crossing_width # -160 by default antimeridian_xmax = 180 - bounding_xmax = antimeridian_xmax - max_scene_width # 160 by default + bounding_xmax = antimeridian_xmax - max_crossing_width # 160 by default - if (bounds.xmin < bounding_xmin) and (bounds.xmin > antimeridian_xmin): - if bounds.xmax > bounding_xmax and bounds.xmax < antimeridian_xmax: + if (bounds.left < bounding_xmin) and (bounds.left > antimeridian_xmin): + if bounds.right > bounding_xmax and bounds.right < antimeridian_xmax: return True return False @@ -255,34 +297,34 @@ def get_target_antimeridian_projection(bounds: BoundingBox) -> int: Parameters ---------- bounds : BoundingBox - The set of bounds (min_lon, min_lat, max_lon, max_lat) + The set of bounds (left, bottom, right, top) Returns ------- int The CRS in integer form (e.g. 3031 for Polar Stereographic) """ - min_lat = min(bounds.ymin, bounds.ymax) + bottom = min(bounds.bottom, bounds.top) target_crs = ( 3031 - if min_lat < -50 - else 3995 if min_lat > 50 else get_local_utm(bounds.bounds, antimeridian=True) + if bottom < -50 + else 3995 if bottom > 50 else get_local_utm(bounds.bounds, antimeridian=True) ) logger.warning(f"Data will be returned in EPSG:{target_crs} projection") return target_crs -def split_s1_bounds_at_am_crossing( +def split_bounds_at_antimeridian( bounds: BBox, lat_buff: float = 0 ) -> tuple[BoundingBox]: - """Split the s1 bounds at the antimeridian, producing one set of bounds for the + """Split the bounds at the antimeridian, producing one set of bounds for the Eastern Hemisphere (left of the antimeridian) and one set for the Western Hemisphere (right of the antimeridian) Parameters ---------- bounds : BBox (BoundingBox | tuple[float | int, float | int, float | int, float | int]) - The set of bounds (xmin, ymin, xmax, ymax) + The set of bounds (left, bottom, right, top) lat_buff : float, optional An additional buffer to subtract from lat, by default 0. @@ -295,23 +337,11 @@ def split_s1_bounds_at_am_crossing( if isinstance(bounds, tuple): bounds = BoundingBox(*bounds) - eastern_hemisphere_x = min([x for x in [bounds.xmin, bounds.xmax] if x > 0]) - if eastern_hemisphere_x > 180: - raise ValueError( - f"Eastern Hemisphere coordinate of {eastern_hemisphere_x} is more than 180 degrees, but should be less." - ) - - western_hemisphere_x = max([x for x in [bounds.xmin, bounds.xmax] if x < 0]) - if western_hemisphere_x < -180: - raise ValueError( - f"Western Hemisphere coordinate of {western_hemisphere_x} is less than -180 degrees, but should be greater." - ) - - min_y = max(-90, bounds.ymin - lat_buff) - max_y = min(90, bounds.ymax + lat_buff) + min_y = max(-90, bounds.bottom - lat_buff) + max_y = min(90, bounds.top + lat_buff) - bounds_western_hemisphere = BoundingBox(-180, min_y, western_hemisphere_x, max_y) - bounds_eastern_hemisphere = BoundingBox(eastern_hemisphere_x, min_y, 180, max_y) + bounds_western_hemisphere = BoundingBox(-180, min_y, bounds.right, max_y) + bounds_eastern_hemisphere = BoundingBox(bounds.left, min_y, 180, max_y) logger.info(f"Eastern Hemisphere bounds: {bounds_eastern_hemisphere.bounds}") logger.info(f"Western Hemisphere bounds: {bounds_western_hemisphere.bounds}") @@ -361,39 +391,6 @@ def coords_from_polygon(polygon): return list(longitudes), list(latitudes) -def get_correct_bounds_from_shape_at_antimeridian( - shape: shapely.geometry.shape, -) -> BoundingBox: - """Get the correct set of bounds for a polygon that crosses the antimeridian. For example - - shape = POLYGON ((178.576126 -71.618423, -178.032867 -70.167343, 176.938004 -68.765106, 173.430893 -70.119957, 178.576126 -71.618423)) - This is a valid shape that nearly stretches the width of the earth, thus shapely is not aware the polygon actually crosses the AM. - By taking bounds with shapely: - - bounds -> (-178.032867, -71.618423, 178.576126, -68.765106) - However, the desired bounds as it crosses the AM are: - - bounds -> (-178.032867, -71.618423, 173.430893, -68.765106) - - Parameters - ---------- - shape : shapely.Polygon.shape - Shape, for example of a sentinel-1 scene footprint. - - Returns - ------- - BoundingBox - The corrected set of bounds with the full span across the antimeridian - """ - - longitudes, latitudes = get_all_lat_lon_coords(shape) - # get furthest point west on the eastern side - min_x = max([x for x in longitudes if x < 0]) - # get furthest point east on the western side - max_x = min([x for x in longitudes if x > 0]) - min_y = min(latitudes) - max_y = max(latitudes) - - return BoundingBox(min_x, min_y, max_x, max_y) - - def adjust_bounds_at_high_lat(bounds: BBox) -> tuple: """Expand the bounds for high latitudes. The provided bounds sometimes do not contain the full scene due to @@ -404,20 +401,20 @@ def adjust_bounds_at_high_lat(bounds: BBox) -> tuple: Parameters ---------- bounds : BBox (BoundingBox | tuple[float | int, float | int, float | int, float | int]) - The set of bounds (min_lon, min_lat, max_lon, max_lat) + The set of bounds (left, bottom, right, top) Returns ------- BoundingBox - The expanded bounds (min_lon, min_lat, max_lon, max_lat) + The expanded bounds (left, bottom, right, top) """ if isinstance(bounds, tuple): bounds = BoundingBox(*bounds) - if bounds.ymin < -50: + if bounds.bottom < -50: logging.info(f"Adjusting bounds at high southern latitudes") bounds = adjust_bounds(bounds, src_crs=4326, ref_crs=3031) - if bounds.ymin > 50: + if bounds.bottom > 50: logging.info(f"Adjusting bounds at high northern latitudes") bounds = adjust_bounds(bounds, src_crs=4326, ref_crs=3995) @@ -439,15 +436,42 @@ def resize_bounds(bounds: BoundingBox, scale_factor: float = 1.0) -> BoundingBox BoundingBox Resized bounding box. """ - x_dim = bounds.xmax - bounds.xmin - y_dim = bounds.ymax - bounds.ymin + + if bounds.left >= bounds.right: + logger.warning( + "The bounding box's left value is greater than or equal to the right value. " + "Assuming the bounds cross the antimeridian. Refactor the bounds if this is not correct." + ) + x_dim = (180 - bounds.left) + (bounds.right - (-180)) # width crossing the am + y_dim = bounds.top - bounds.bottom + else: + x_dim = bounds.right - bounds.left + y_dim = bounds.top - bounds.bottom dx = ((scale_factor - 1) * x_dim) / 2 dy = ((scale_factor - 1) * y_dim) / 2 - return BoundingBox( - bounds.xmin - dx, bounds.ymin - dy, bounds.xmax + dx, bounds.ymax + dy - ) + # adjust the bounds by requested factor + new_left = bounds.left - dx + new_bottom = bounds.bottom - dy + new_right = bounds.right + dx + new_top = bounds.top + dy + + resized_bounds = BoundingBox(new_left, new_bottom, new_right, new_top) + + if not resized_bounds._check_valid(): + logger.warning( + "Resized bounds are not valid. " + f"Bounds : {bounds.bounds}, scale_factor : {scale_factor} -> Resized Bounds : {resized_bounds.bounds}. " + ) + # Clamp each coordinate to valid geographic limits + resized_bounds.left = max(-180, min(180, resized_bounds.left)) + resized_bounds.right = max(-180, min(180, resized_bounds.right)) + resized_bounds.bottom = max(-90, min(90, resized_bounds.bottom)) + resized_bounds.top = max(-90, min(90, resized_bounds.top)) + f"Setting to maximum allowable extents : {resized_bounds.bounds}" + + return resized_bounds def crop_datasets_to_bounds( diff --git a/tests/cop_glo30/test_dem_cop.py b/tests/cop_glo30/test_dem_cop.py index 1c9eb25..554af90 100644 --- a/tests/cop_glo30/test_dem_cop.py +++ b/tests/cop_glo30/test_dem_cop.py @@ -46,17 +46,22 @@ class TestDem: / "Copernicus_DSM_COG_10_two_tiles_same_latitude_S71_1155_2603_171_131.tif", ) +test_antimeridian_dem = TestDem( + (178.5, -78.1, -179.1, -77.6), + TEST_DATA_PATH / "Copernicus_DSM_COG_10_antimeridian.tif", +) test_dems = [ test_one_land_tile_with_ocean, test_one_land_tile_no_ocean, test_dem_three_land_tiles_and_one_ocean_tile, test_dem_two_tiles_same_latitude, + test_antimeridian_dem, ] @pytest.mark.parametrize("test_input", test_dems) -def test_local_get_cop30_dem_for_bounds_ocean_and_land(test_input: TestDem): +def test_local_get_cop30_dem_for_bounds(test_input: TestDem): bounds = test_input.requested_bounds bounds_array_file = test_input.bounds_array_file @@ -68,7 +73,7 @@ def test_local_get_cop30_dem_for_bounds_ocean_and_land(test_input: TestDem): INDEX_PATH = TMP_PATH / Path("TMP.gpkg") # Find relevant test tiles and build tile index - TEST_TILES = find_tiles(TEST_DATA_PATH, "Copernicus_DSM_COG_10_S??_00_E16?_00_DEM") + TEST_TILES = find_tiles(TEST_DATA_PATH, "Copernicus_DSM_COG_10_???_00_?1??_00_DEM") build_tileindex( TEST_TILES, INDEX_PATH, @@ -102,7 +107,7 @@ def test_local_get_cop30_dem_for_bounds_ocean_and_land(test_input: TestDem): @pytest.mark.parametrize("test_input", test_dems) -def test_download_get_cop30_dem_for_bounds_ocean_and_land(test_input: TestDem): +def test_download_get_cop30_dem_for_bounds(test_input: TestDem): bounds = test_input.requested_bounds bounds_array_file = test_input.bounds_array_file @@ -117,7 +122,7 @@ def test_download_get_cop30_dem_for_bounds_ocean_and_land(test_input: TestDem): TEST_DOWNLOAD_PATH.mkdir(parents=True, exist_ok=True) # Find relevant test tiles and build tile index - TEST_TILES = find_tiles(TEST_DATA_PATH, "Copernicus_DSM_COG_10_S??_00_E16?_00_DEM") + TEST_TILES = find_tiles(TEST_DATA_PATH, "Copernicus_DSM_COG_10_???_00_?1??_00_DEM") build_tileindex( TEST_TILES, INDEX_PATH, diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index b3a0793..b81a85c 100644 --- a/tests/utils/test_utils.py +++ b/tests/utils/test_utils.py @@ -3,11 +3,96 @@ 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_bounds_likely_cross_antimeridian, check_dem_type_in_bounds, + resize_bounds, ) + +@dataclass +class TestValidBounds: + bounds: BoundingBox + lat_tol: float + lon_tol: float + is_valid: int + + +# Test cases for bounding box validity +TEST_VALID_BOUNDS = [ + # Fully valid within standard ranges + TestValidBounds( + BoundingBox(-170, 0, -165, 5), lat_tol=0.0, lon_tol=0.0, is_valid=True + ), + TestValidBounds( + BoundingBox(0, 10, 160, 50), lat_tol=0.0, lon_tol=0.0, is_valid=True + ), + # Crosses antimeridian but still valid logically + TestValidBounds( + BoundingBox(170, -70, -170, -65), lat_tol=0.0, lon_tol=0.0, is_valid=True + ), + # Longitude slightly out of range, invalid without tolerance + TestValidBounds( + BoundingBox(-181, 0, -170, 15), lat_tol=0.0, lon_tol=0.0, is_valid=False + ), + # Latitude slightly out of range, invalid without tolerance + TestValidBounds( + BoundingBox(0, 0, 30, 92), lat_tol=0.0, lon_tol=0.0, is_valid=False + ), + # Becomes valid with longitude tolerance + TestValidBounds( + BoundingBox(-181, 0, -170, 15), lat_tol=0.0, lon_tol=2.0, is_valid=True + ), + # Becomes valid with latitude tolerance + TestValidBounds(BoundingBox(0, 0, 30, 92), lat_tol=3.0, lon_tol=0.0, is_valid=True), +] + + +@pytest.mark.parametrize("case", TEST_VALID_BOUNDS) +def test_check_valid_bounding_box(case: TestValidBounds): + """Test BoundingBox validity logic with and without tolerance.""" + result = case.bounds._check_valid(lat_tol=case.lat_tol, lon_tol=case.lon_tol) + assert result == case.is_valid, ( + f"Expected {case.is_valid} for bounds {case.bounds} " + f"with lat_tol={case.lat_tol}, lon_tol={case.lon_tol}, got {result}" + ) + + +@dataclass +class TestResizeBounds: + bounds: BoundingBox + scale_factor: BoundingBox + resized_bounds: BoundingBox + + +TEST_RESIZE_BOUNDS = [ + # normal bounds + TestResizeBounds( + BoundingBox(-170, 0, -165, 5), 2, BoundingBox(-172.5, -2.5, -162.5, 7.5) + ), + TestResizeBounds( + BoundingBox(-170, 0, -166, 6), + 0.5, + BoundingBox(-169, 1.5, -167, 4.5), + ), + # antimeridian bounds + TestResizeBounds( + BoundingBox(178, -54, -178, -50), + 0.5, + BoundingBox(179, -53, -179, -51), + ), + # scale by too much, ensure set to max allowable extents + TestResizeBounds( + BoundingBox(-170, 0, -165, 5), 8, BoundingBox(-180, -17.5, -147.5, 22.5) + ), +] + + +@pytest.mark.parametrize("case", TEST_RESIZE_BOUNDS) +def test_resize_bounds(case: TestResizeBounds): + resized_bounds = resize_bounds(case.bounds, scale_factor=case.scale_factor) + assert resized_bounds == case.resized_bounds + + # Separate lists for clarity poly1 = Polygon( [ @@ -41,23 +126,11 @@ 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): +@pytest.mark.parametrize("shape", shapes) +def test_get_correct_bounds_from_shape_at_antimeridian(shape): # 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 + assert check_bounds_likely_cross_antimeridian(shape.bounds) @dataclass @@ -128,7 +201,7 @@ class BoundsDEMCheckCase: test_antimeridian_with_cop30 = BoundsDEMCheckCase( dem_type="cop_glo30", resolution=30, - bounds=(-178.032867, -80.618423, 173.430893, -68.765106), + bounds=(173.430893, -80.618423, -178.032867, -68.765106), in_bounds=True, is_error=False, ) @@ -136,7 +209,7 @@ class BoundsDEMCheckCase: test_antimeridian_with_rema = BoundsDEMCheckCase( dem_type="rema", resolution=2, - bounds=(-178.032867, -80.618423, 173.430893, -68.765106), + bounds=(173.430893, -80.618423, -178.032867, -68.765106), in_bounds=True, is_error=False, )