From dd531ca24c9ee966343eedde1a7f11339b475355 Mon Sep 17 00:00:00 2001 From: nlenglet Date: Tue, 21 Apr 2026 15:37:06 +0200 Subject: [PATCH] feat: add option to compare dtm with lidarhd stream --- altianalysis/compute_difference.py | 23 ++++++++++---- altianalysis/run_difference_with_gpao.py | 23 ++++++++++---- test/test_compute_difference.py | 38 ++++++++++++++++++++++-- test/test_run_difference_with_gpao.py | 10 ++++--- 4 files changed, 76 insertions(+), 18 deletions(-) diff --git a/altianalysis/compute_difference.py b/altianalysis/compute_difference.py index ceccab4..f28d61a 100644 --- a/altianalysis/compute_difference.py +++ b/altianalysis/compute_difference.py @@ -23,13 +23,14 @@ def parse_args(): "If not provided, primary elevation file is compared with rge alti (1 meter) data stream", ) parser.add_argument("--name_save_out", type=str, help="name of difference file") + parser.add_argument("--use_lidarHD", action="store_true", help="use lidarHD data stream") return parser.parse_args() -def _extract_rge_alti_tile_from_stream( +def _extract_tiles_from_stream( dtm_file: str, output_path: str | Path, - stream_RGE="RGEALTI-MNT_PYR-ZIP_FXX_LAMB93_WMS", # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES", + stream="RGEALTI-MNT_PYR-ZIP_FXX_LAMB93_WMS", # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES", proj="2154", pixel_per_meter=1, timeout_second=300, @@ -44,7 +45,7 @@ def _extract_rge_alti_tile_from_stream( try: download_image( proj, - layer=stream_RGE, + layer=stream, minx=minx, miny=miny, maxx=maxx, @@ -100,16 +101,26 @@ def compute_difference_between_dtms(first_dtm_file: str, second_dtm_file: str, n dst.write(_difference, 1) -def main(reference_dtm_file: Path | str, secondary_dtm_file: Path | str | None, output_difference_file: Path | str): +def main(reference_dtm_file: Path | str, secondary_dtm_file: Path | str | None, output_difference_file: Path | str, use_lidarHD: bool = False): if secondary_dtm_file: compute_difference_between_dtms(reference_dtm_file, secondary_dtm_file, output_difference_file) + elif use_lidarHD: + + with tempfile.NamedTemporaryFile(suffix="_dtm_lidarHD.tif", delete=False) as tmp_lidarHD: + success = _extract_tiles_from_stream( + reference_dtm_file, pixel_per_meter=5, output_path=tmp_lidarHD.name, stream="IGNF_LIDAR-HD_MNT_ELEVATION.ELEVATIONGRIDCOVERAGE.LAMB93" + ) + + if success: + compute_difference_between_dtms(reference_dtm_file, tmp_lidarHD, output_difference_file) + else: with tempfile.NamedTemporaryFile(suffix="_dtm_rgealti.tif", delete=False) as tmp_rge: - success = _extract_rge_alti_tile_from_stream( + success = _extract_tiles_from_stream( reference_dtm_file, pixel_per_meter=5, output_path=tmp_rge.name ) @@ -119,4 +130,4 @@ def main(reference_dtm_file: Path | str, secondary_dtm_file: Path | str | None, if __name__ == "__main__": args = parse_args() - main(args.primary_elevation_file, args.second_elevation_file, args.name_save_out) + main(args.primary_elevation_file, args.second_elevation_file, args.name_save_out, args.use_lidarHD) diff --git a/altianalysis/run_difference_with_gpao.py b/altianalysis/run_difference_with_gpao.py index 000c3e3..183850d 100644 --- a/altianalysis/run_difference_with_gpao.py +++ b/altianalysis/run_difference_with_gpao.py @@ -26,7 +26,7 @@ def get_tile_names(folder: Path) -> List[str]: return filenames -def create_one_job_one_difference(store: Store, dir_in: Path, second_dir: Path | None, input_file: str, output: Path): +def create_one_job_one_difference(store: Store, dir_in: Path, second_dir: Path | None, input_file: str, output: Path, use_lidarHD: bool): job_name = f"difference_{input_file}" if second_dir: secondary_file = second_dir / input_file @@ -42,6 +42,8 @@ def create_one_job_one_difference(store: Store, dir_in: Path, second_dir: Path | second_input_mount = "" second_input_param = "" + use_lidarHD_param = "--use_lidarHD" if use_lidarHD else "" + command = f""" docker run -t --rm --userns=host -v {store.to_unix(dir_in)}:/input @@ -51,6 +53,7 @@ def create_one_job_one_difference(store: Store, dir_in: Path, second_dir: Path | python -m altianalysis.compute_difference --primary_elevation_file /input/{input_file} {second_input_param} + {use_lidarHD_param} --name_save_out /output/{input_file} """ @@ -64,6 +67,7 @@ def create_main_gpao_project( out: Path, store: Store, project_name: str, + use_lidarHD: bool, ) -> Project: if secondary_dir: @@ -88,7 +92,7 @@ def create_main_gpao_project( jobs = [] for tile_ in dtm_tile_names: - job = create_one_job_one_difference(store, primary_dir, secondary_dir, tile_, out) + job = create_one_job_one_difference(store, primary_dir, secondary_dir, tile_, out, use_lidarHD) jobs.append(job) # create the project @@ -128,9 +132,9 @@ def create_cog_gpao_project( def create_gpao_projects( - primary_dir: Path, secondary_dir: Path | None, out: Path, store: Store, project_name: str, cog_filename: str + primary_dir: Path, secondary_dir: Path | None, out: Path, store: Store, project_name: str, cog_filename: str, use_lidarHD: bool ) -> List[Project]: - project_main = create_main_gpao_project(primary_dir, secondary_dir, out, store, project_name) + project_main = create_main_gpao_project(primary_dir, secondary_dir, out, store, project_name, use_lidarHD) projects = [project_main] if cog_filename: project_cog = create_cog_gpao_project( @@ -144,6 +148,7 @@ def create_gpao_projects( def compute_on_gpao( primary_dir: Path, secondary_dir: Path | None, + use_lidarHD: bool, out: Path, gpao_hostname: str, local_store_path: Path, @@ -171,7 +176,7 @@ def compute_on_gpao( logging.debug(f"Local store path ({local_store_path}) converted to client store path ({runner_store_path})") - projects = create_gpao_projects(primary_dir, secondary_dir, out, store, project_name, cog_filename) + projects = create_gpao_projects(primary_dir, secondary_dir, out, store, project_name, cog_filename, use_lidarHD) builder = Builder(projects) @@ -206,6 +211,13 @@ def parse_args(): "Les dalles d'élévation doivent avoir les mêmes noms pour faire l'appariement" "S'il est laissé vide, les dalles du premier ensemble sont comparées au RGE Alti", ) + parser.add_argument( + "-iii", + "--use_lidarHD", + type=bool, + default=False, + help="Utiliser les données LIDARHD pour le calcul de différence", + ) parser.add_argument( "-o", @@ -251,6 +263,7 @@ def parse_args(): compute_on_gpao( primary_dir=args.primary_dtm_dir, secondary_dir=args.secondary_dtm_dir, + use_lidarHD=args.use_lidarHD, out=args.out, gpao_hostname=args.gpao_hostname, local_store_path=args.local_store_path, diff --git a/test/test_compute_difference.py b/test/test_compute_difference.py index 3dc0c1c..b29ad82 100644 --- a/test/test_compute_difference.py +++ b/test/test_compute_difference.py @@ -24,7 +24,7 @@ def test_extract_rge_alti_tile_from_stream(): output_dir.mkdir() dtm_rge_alti = output_dir / "dtm_rge_alti.tif" dtm_lidar_file = "./data/lhd/Semis_2021_0886_6443_LA93_IGN69_50CM.tif" - compute_difference._extract_rge_alti_tile_from_stream(dtm_lidar_file, dtm_rge_alti) + compute_difference._extract_tiles_from_stream(dtm_lidar_file, dtm_rge_alti) # read dtm rge alti and dtm_lidar_file and check bounds with rasterio.open(dtm_lidar_file) as dtm_lidar, rasterio.open(dtm_rge_alti) as dtm_rge: @@ -40,6 +40,26 @@ def test_extract_rge_alti_tile_from_stream(): {bounds_dtm_rge.top})" +def test_extract_lidarHD_tile_from_stream(): + output_dir = TMP_PATH / "extract_lidarHD_tile_from_stream" + output_dir.mkdir() + dtm_output_file = output_dir / "dtm_lidarHD.tif" + dtm_lidar_file = "./data/lhd/Semis_2021_0886_6443_LA93_IGN69_50CM.tif" + compute_difference._extract_tiles_from_stream(dtm_lidar_file, dtm_output_file, stream="IGNF_LIDAR-HD_MNT_ELEVATION.ELEVATIONGRIDCOVERAGE.LAMB93") + + # read lidarHD file and check bounds + with rasterio.open(dtm_lidar_file) as dtm_lidar, rasterio.open(dtm_output_file) as dtm_output: + bounds_dtm_lidar = dtm_lidar.bounds + bounds_dtm_output = dtm_output.bounds + + assert ( + bounds_dtm_lidar == bounds_dtm_output + ), f"tiles extents are not matching: lidar dtm has ({bounds_dtm_lidar.left}, \ + {bounds_dtm_lidar.right}, {bounds_dtm_lidar.bottom}, {bounds_dtm_lidar.top}), \ + lidarHD dtm has ({bounds_dtm_output.left}, {bounds_dtm_output.right}, {bounds_dtm_output.bottom}, \ + {bounds_dtm_output.top})" + + def test_compute_difference_between_dtms_with_self(): output_dir = TMP_PATH / "compute_difference_between_dtms_with_self" output_dir.mkdir() @@ -61,7 +81,7 @@ def test_compute_difference_and_reconstruct(): dtm_rge_alti = output_dir / "dtm_rge_alti.tif" out_difference_file = output_dir / "Difference_Semis_2021_0886_6443_LA93_IGN69_50CM.tif" out_recomputed_lidar_file = output_dir / "Rebuilt_Semis_2021_0886_6443_LA93_IGN69_50CM.tif" - compute_difference._extract_rge_alti_tile_from_stream(dtm_lidar_file, dtm_rge_alti) + compute_difference._extract_tiles_from_stream(dtm_lidar_file, dtm_rge_alti) compute_difference.compute_difference_between_dtms(dtm_lidar_file, dtm_rge_alti, out_difference_file) # read files and check if we can reconstruct initial lidar image @@ -115,7 +135,7 @@ def test_compute_difference_with_nodata(): # lidar dtm with nodata dtm_lidar_file = "./data/lhd/Semis_2021_0485_6196_LA93_IGN69_50CM.tif" dtm_rge_alti = output_dir / "dtm_rge_alti.tif" - compute_difference._extract_rge_alti_tile_from_stream(dtm_lidar_file, dtm_rge_alti) + compute_difference._extract_tiles_from_stream(dtm_lidar_file, dtm_rge_alti) out_difference_file = output_dir / "Difference_Semis_2021_0485_6196_LA93_IGN69_50CM.tif" compute_difference.compute_difference_between_dtms(dtm_lidar_file, dtm_rge_alti, out_difference_file) @@ -190,3 +210,15 @@ def test_main_with_secondary_folder(): with rasterio.open(out_difference_file_with_self) as o_diff_self: diff = o_diff_self.read() assert np.all(diff == 0), "difference with self yields non null values" + + +def test_main_with_lidarHD(): + output_dir = TMP_PATH / "main_with_lidarHD" + output_dir.mkdir() + dtm_lidar_file = "./data/lhd/Semis_2021_0886_6443_LA93_IGN69_50CM.tif" + + out_difference_file = output_dir / "Difference_with_lidarHD_Semis_2021_0886_6443_LA93_IGN69_50CM.tif" + + compute_difference.main(dtm_lidar_file, None, out_difference_file, use_lidarHD=True) + + assert os.path.exists(out_difference_file), "difference with lidarHD not computed !" \ No newline at end of file diff --git a/test/test_run_difference_with_gpao.py b/test/test_run_difference_with_gpao.py index cca2741..cdaad8f 100644 --- a/test/test_run_difference_with_gpao.py +++ b/test/test_run_difference_with_gpao.py @@ -27,7 +27,7 @@ def test_create_main_gpao_project(): output_dir = TMP_PATH / "create_main_gpao_project" dtm_lidar_lhds = Path("./data/lhd_dir_gpao") project_name = "test_create_gpao_project_difference_with_dem_rge_alti" - project = run_difference_with_gpao.create_main_gpao_project(dtm_lidar_lhds, None, output_dir, STORE, project_name) + project = run_difference_with_gpao.create_main_gpao_project(dtm_lidar_lhds, None, output_dir, STORE, project_name, False) assert project is not None @@ -43,7 +43,7 @@ def test_create_gpao_projects_with_cog(): cog_filename = "cog.tif" project_name = "test_create_gpao_project_difference_with_dem_rge_alti" projects = run_difference_with_gpao.create_gpao_projects( - dtm_lidar_lhds, None, output_dir, STORE, project_name, cog_filename + dtm_lidar_lhds, None, output_dir, STORE, project_name, cog_filename, False ) assert len(projects) == 2 @@ -56,7 +56,7 @@ def test_create_gpao_projects_without_cog(): output_dir = TMP_PATH / "create_gpao_projects_without_cog" dtm_lidar_lhds = Path("./data/lhd_dir_gpao") project_name = "test_create_gpao_project_difference_with_dem_rge_alti" - projects = run_difference_with_gpao.create_gpao_projects(dtm_lidar_lhds, None, output_dir, STORE, project_name, "") + projects = run_difference_with_gpao.create_gpao_projects(dtm_lidar_lhds, None, output_dir, STORE, project_name, "", False) assert len(projects) == 1 assert len(projects[0].jobs) == 5 @@ -79,6 +79,7 @@ def test_gpao_run_with_cog_stream_rge(): run_difference_with_gpao.compute_on_gpao( Path(dtm_lidar_lhds), None, + False, Path(output_dir), gpao_hostname, local_store_path, @@ -112,7 +113,7 @@ def test_gpao_run_without_cog_stream_rge(): local_store_path = Path("data/lhd_dir_gpao").resolve() run_difference_with_gpao.compute_on_gpao( - Path(dtm_lidar_lhds), None, Path(output_dir), gpao_hostname, local_store_path, runner_store_path, project_name + Path(dtm_lidar_lhds), None, False, Path(output_dir), gpao_hostname, local_store_path, runner_store_path, project_name ) if gpao_hostname == "localhost": @@ -137,6 +138,7 @@ def test_gpao_run_without_cog_given_secondary_folder(): run_difference_with_gpao.compute_on_gpao( Path(dtm_lidar_lhds), Path(secondary_dtm_dir), + False, Path(output_dir), gpao_hostname, local_store_path,