From db067f68da4e3de72abb946107ca37d94d517deb Mon Sep 17 00:00:00 2001 From: Tommy Date: Mon, 18 Nov 2024 10:05:57 +0100 Subject: [PATCH 1/2] feat: simplification of the check inputs --- cars/core/projection.py | 137 ++++++++++++++++++---------------------- 1 file changed, 60 insertions(+), 77 deletions(-) diff --git a/cars/core/projection.py b/cars/core/projection.py index 01a8dacf..3bca5572 100644 --- a/cars/core/projection.py +++ b/cars/core/projection.py @@ -45,7 +45,7 @@ def compute_dem_intersection_with_poly( # noqa: C901 - srtm_dir: str, ref_poly: Polygon, ref_epsg: int + srtm_file: str, ref_poly: Polygon, ref_epsg: int ) -> Polygon: """ Compute the intersection polygon between the defined dem regions @@ -53,7 +53,7 @@ def compute_dem_intersection_with_poly( # noqa: C901 :raise Exception: when the input dem doesn't intersect the reference polygon - :param srtm_dir: srtm directory + :param srtm_file: srtm file :param ref_poly: reference polygon :param ref_epsg: reference epsg code :return: The intersection polygon between the defined dem regions @@ -61,84 +61,67 @@ def compute_dem_intersection_with_poly( # noqa: C901 """ dem_poly = None - srtm_files_dir = [] - if os.path.isfile(os.path.abspath(srtm_dir)): - srtm_files_dir.append([srtm_dir]) - else: - for _, _, srtm_files in os.walk(srtm_dir): - srtm_files_dir.append(srtm_files) - - for srtm_files in srtm_files_dir: - logging.info( - "Browsing all files of the srtm dir. " - "Some files might be unreadable by rasterio (non blocking matter)." - ) - for file in srtm_files: - unsupported_formats = [".omd"] - _, ext = os.path.splitext(file) - if ext not in unsupported_formats: - if inputs.rasterio_can_open(os.path.join(srtm_dir, file)): - with rio.open(os.path.join(srtm_dir, file)) as data: - xmin = min(data.bounds.left, data.bounds.right) - ymin = min(data.bounds.bottom, data.bounds.top) - xmax = max(data.bounds.left, data.bounds.right) - ymax = max(data.bounds.bottom, data.bounds.top) - - try: - file_epsg = data.crs.to_epsg() - file_bb = Polygon( - [ - (xmin, ymin), - (xmin, ymax), - (xmax, ymax), - (xmax, ymin), - (xmin, ymin), - ] - ) - - # transform polygon if needed + if os.path.isdir(os.path.abspath(srtm_file)): + raise RuntimeError("srtm input should be a file and not a directory") + + if inputs.rasterio_can_open(srtm_file): + with rio.open(srtm_file) as data: + xmin = min(data.bounds.left, data.bounds.right) + ymin = min(data.bounds.bottom, data.bounds.top) + xmax = max(data.bounds.left, data.bounds.right) + ymax = max(data.bounds.bottom, data.bounds.top) + + try: + file_epsg = data.crs.to_epsg() + file_bb = Polygon( + [ + (xmin, ymin), + (xmin, ymax), + (xmax, ymax), + (xmax, ymin), + (xmin, ymin), + ] + ) + + # transform polygon if needed + if ref_epsg != file_epsg: + file_bb = polygon_projection(file_bb, file_epsg, ref_epsg) + + # if the srtm tile intersects the reference polygon + if file_bb.intersects(ref_poly): + local_dem_poly = None + + # retrieve valid polygons + for poly, val in shapes( + data.dataset_mask(), + transform=data.transform, + ): + if val != 0: + poly = shape(poly) + poly = poly.buffer(0) if ref_epsg != file_epsg: - file_bb = polygon_projection( - file_bb, file_epsg, ref_epsg + poly = polygon_projection( + poly, file_epsg, ref_epsg ) - # if the srtm tile intersects the reference polygon - if file_bb.intersects(ref_poly): - local_dem_poly = None - - # retrieve valid polygons - for poly, val in shapes( - data.dataset_mask(), - transform=data.transform, - ): - if val != 0: - poly = shape(poly) - poly = poly.buffer(0) - if ref_epsg != file_epsg: - poly = polygon_projection( - poly, file_epsg, ref_epsg - ) - - # combine valid polygons - if local_dem_poly is None: - local_dem_poly = poly - else: - local_dem_poly = poly.union( - local_dem_poly - ) - - # combine the tile valid polygon to the other - # tiles' ones - if dem_poly is None: - dem_poly = local_dem_poly - else: - dem_poly = dem_poly.union(local_dem_poly) - - except AttributeError as attribute_error: - logging.warning( - "Impossible to read the SRTM" - "tile epsg code: {}".format(attribute_error) - ) + # combine valid polygons + if local_dem_poly is None: + local_dem_poly = poly + else: + local_dem_poly = poly.union(local_dem_poly) + + # combine the tile valid polygon to the other + # tiles' ones + if dem_poly is None: + dem_poly = local_dem_poly + else: + dem_poly = dem_poly.union(local_dem_poly) + + except AttributeError as attribute_error: + logging.warning( + "Impossible to read the SRTM" + "tile epsg code: {}".format(attribute_error) + ) # compute dem coverage polygon over the reference polygon if dem_poly is None or not dem_poly.intersects(ref_poly): From 8b4fc0e048befb8ab4807ebfa10a9063c6b5ef3c Mon Sep 17 00:00:00 2001 From: Tommy Date: Mon, 18 Nov 2024 10:25:56 +0100 Subject: [PATCH 2/2] test: modification of test_compute_dem_intersection_with_poly --- tests/core/test_projection.py | 32 ++++++-------------------------- 1 file changed, 6 insertions(+), 26 deletions(-) diff --git a/tests/core/test_projection.py b/tests/core/test_projection.py index b44b5710..239e1594 100644 --- a/tests/core/test_projection.py +++ b/tests/core/test_projection.py @@ -90,34 +90,12 @@ def test_compute_dem_intersection_with_poly(): ) dem_inter_poly, cover = projection.compute_dem_intersection_with_poly( - absolute_data_path("input/phr_ventoux/srtm"), inter_poly, inter_epsg - ) - assert dem_inter_poly.equals(inter_poly) - assert cover == 100.0 - - # test partial coverage over with several srtm tiles with no data holes - inter_poly = Polygon( - [(4.8, 44.2), (4.8, 44.3), (6.2, 44.3), (6.2, 44.2), (4.8, 44.2)] - ) - dem_inter_poly, cover = projection.compute_dem_intersection_with_poly( - absolute_data_path("input/utils_input/srtm_with_hole"), + absolute_data_path("input/phr_ventoux/srtm/N44E005.hgt"), inter_poly, inter_epsg, ) - - ref_dem_inter_poly = Polygon( - [ - (4.999583333333334, 44.3), - (6.2, 44.3), - (6.2, 44.2), - (4.999583333333334, 44.2), - (4.999583333333334, 44.3), - ] - ) - - assert ref_dem_inter_poly.exterior.equals(dem_inter_poly.exterior) - assert len(list(dem_inter_poly.interiors)) == 6 - assert cover == 85.72172619047616 + assert dem_inter_poly.equals(inter_poly) + assert cover == 100.0 # test no coverage inter_poly = Polygon( @@ -126,7 +104,9 @@ def test_compute_dem_intersection_with_poly(): with pytest.raises(Exception) as intersect_error: dem_inter_poly, cover = projection.compute_dem_intersection_with_poly( - absolute_data_path("input/phr_ventoux/srtm"), inter_poly, inter_epsg + absolute_data_path("input/phr_ventoux/srtm/N44E005.hgt"), + inter_poly, + inter_epsg, ) assert ( str(intersect_error.value) == "The input DEM does not intersect "