diff --git a/dem_handler/utils/spatial.py b/dem_handler/utils/spatial.py index 49153ae..f29d7df 100644 --- a/dem_handler/utils/spatial.py +++ b/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 diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py index b81a85c..8f2b102 100644 --- a/tests/utils/test_utils.py +++ b/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, ]