From adbf1a90d55aeff1acd1b481334e1c74f886f4d3 Mon Sep 17 00:00:00 2001 From: DaliCHEBBI Date: Wed, 24 Sep 2025 18:26:03 +0200 Subject: [PATCH 1/6] add feature of difference between two folders --- altianalysis/compute_difference.py | 22 ++++++--- altianalysis/run_difference_with_gpao.py | 48 +++++++++++++++----- altianalysis/run_difference_with_parallel.py | 31 ++++++++----- test/test_compute_difference.py | 14 +++++- test/test_run_difference_with_gpao.py | 4 +- test/test_run_difference_with_parallel.py | 17 ++++++- 6 files changed, 102 insertions(+), 34 deletions(-) diff --git a/altianalysis/compute_difference.py b/altianalysis/compute_difference.py index d4c76ad..60e41ac 100644 --- a/altianalysis/compute_difference.py +++ b/altianalysis/compute_difference.py @@ -13,6 +13,9 @@ def parse_args(): parser = argparse.ArgumentParser("compute difference map between lidar dtm and rge alti") parser.add_argument("--dtm_lidar_file", type=str, help="dtm lidar tif file") + parser.add_argument( + "--second_elevation_file", default=None, type=str, help="secondary elevation file when no stream to be used !" + ) parser.add_argument("--name_save_out", type=str, help="name of difference file") return parser.parse_args() @@ -91,16 +94,23 @@ def compute_difference_between_dtms(dtm_lidar_file: str, dtm_rge_tile: str, name dst.write(_difference, 1) -def main(reference_dtm_file: Path | str, output_difference_file: Path | str): +def main(reference_dtm_file: Path | str, secondary_dtm_file: Path | str | None, output_difference_file: Path | str): - with tempfile.NamedTemporaryFile(suffix="_dtm_rgealti.tif", delete=False) as tmp_rge: + if secondary_dtm_file: + compute_difference_between_dtms(reference_dtm_file, secondary_dtm_file, output_difference_file) - success = _extract_rge_alti_tile_from_stream(reference_dtm_file, pixel_per_meter=5, output_path=tmp_rge.name) + else: - if success: - compute_difference_between_dtms(reference_dtm_file, tmp_rge, output_difference_file) + with tempfile.NamedTemporaryFile(suffix="_dtm_rgealti.tif", delete=False) as tmp_rge: + + success = _extract_rge_alti_tile_from_stream( + reference_dtm_file, pixel_per_meter=5, output_path=tmp_rge.name + ) + + if success: + compute_difference_between_dtms(reference_dtm_file, tmp_rge, output_difference_file) if __name__ == "__main__": args = parse_args() - main(args.dtm_lidar_file, args.name_save_out) + main(args.dtm_lidar_file, args.second_elevation_file, args.name_save_out) diff --git a/altianalysis/run_difference_with_gpao.py b/altianalysis/run_difference_with_gpao.py index 60311dd..9f2d405 100644 --- a/altianalysis/run_difference_with_gpao.py +++ b/altianalysis/run_difference_with_gpao.py @@ -25,23 +25,37 @@ def get_tile_names(folder: Path) -> List[str]: return filenames -def create_one_job_one_difference(store: Store, dir_in: Path, 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): job_name = f"difference_{input_file}" - command = f""" -docker run -t --rm --userns=host --shm-size=2gb --v {store.to_unix(dir_in)}:/input --v {store.to_unix(output)}:/output -ghcr.io/ignf/altianalysis:{__version__} -python -m altianalysis.compute_difference ---dtm_lidar_file /input/{input_file} ---name_save_out /output/{input_file} -""" + if second_dir is not None: + command = f""" + docker run -t --rm --userns=host --shm-size=2gb + -v {store.to_unix(dir_in)}:/input + -v {store.to_unix(second_dir)}:/second_input + -v {store.to_unix(output)}:/output + ghcr.io/ignf/altianalysis:{__version__} + python -m altianalysis.compute_difference + --dtm_lidar_file /input/{input_file} + --second_elevation_file /second_input{input_file} + --name_save_out /output/{input_file} + """ + else: + command = f""" + docker run -t --rm --userns=host --shm-size=2gb + -v {store.to_unix(dir_in)}:/input + -v {store.to_unix(output)}:/output + ghcr.io/ignf/altianalysis:{__version__} + python -m altianalysis.compute_difference + --dtm_lidar_file /input/{input_file} + --name_save_out /output/{input_file} + """ job = Job(job_name, command, tags=["docker"]) return job def create_gpao_project( dtms_lhd: Path, + secondary_dir: Path | None, out: Path, store: Store, project_name: str, @@ -63,7 +77,7 @@ def create_gpao_project( jobs = [] for tile_ in dtm_tile_names: - job = create_one_job_one_difference(store, dtms_lhd, tile_, out) + job = create_one_job_one_difference(store, dtms_lhd, secondary_dir, tile_, out) jobs.append(job) # create the project @@ -72,6 +86,7 @@ def create_gpao_project( def compute_on_gpao( dtms_lhd: Path, + secondary_dir: Path | None, out: Path, gpao_hostname: str, local_store_path: Path, @@ -85,7 +100,7 @@ def compute_on_gpao( logging.debug(f"Local store path ({local_store_path}) converted to client store path ({runner_store_path})") - project = create_gpao_project(dtms_lhd, out, store, project_name) + project = create_gpao_project(dtms_lhd, secondary_dir, out, store, project_name) builder = Builder([project]) builder.save_as_json(out / "gpao_project.json") @@ -104,6 +119,14 @@ def parse_args(): required=True, help="Dossier des dalles MNT Lidar HD", ) + parser.add_argument( + "-i", + "--secondary_dtm_dir", + type=Path, + default=None, + help="Dossier du second des MNT pour calculer la différence", + ) + parser.add_argument( "-o", "--out", @@ -139,6 +162,7 @@ def parse_args(): compute_on_gpao( args.dtm_lhd_dir, + args.secondary_dtm_dir, args.out, args.gpao_hostname, args.local_store_path, diff --git a/altianalysis/run_difference_with_parallel.py b/altianalysis/run_difference_with_parallel.py index db76947..8ba2df1 100644 --- a/altianalysis/run_difference_with_parallel.py +++ b/altianalysis/run_difference_with_parallel.py @@ -1,38 +1,36 @@ import argparse import os -import tempfile from pathlib import Path from joblib import Parallel, delayed -from altianalysis.compute_difference import _extract_rge_alti_tile_from_stream, compute_difference_between_dtms +from altianalysis.compute_difference import main -def _compute_one_difference(dtm_lhd_file: str, _dir: Path, _out_dir_difference: Path): +def _compute_one_difference(dtm_lhd_file: str, _dir: Path, _out_dir_difference: Path, second_dir=None): full_dtm_file = os.path.join(_dir, dtm_lhd_file) full_difference_file = os.path.join(_out_dir_difference, dtm_lhd_file) + _second_file = None - # compute difference for individual files - with tempfile.NamedTemporaryFile(suffix="_dtm_rgealti.tif", delete=False) as tmp_rge: + if second_dir is not None: + _second_file = os.path.join(second_dir, dtm_lhd_file) - success = _extract_rge_alti_tile_from_stream(full_dtm_file, tmp_rge.name) - if success: - compute_difference_between_dtms(full_dtm_file, tmp_rge.name, full_difference_file) + main(full_dtm_file, _second_file, full_difference_file) -def compute_all_difference_maps(_dir: Path, _out_dir_difference: Path): +def compute_all_difference_maps(_dir: Path, _second_dir: Path | None, _out_dir_difference: Path): os.makedirs(_out_dir_difference, exist_ok=True) all_dtm_lhd_names = [] - for dtm_file in _dir.iterdir(): if dtm_file.is_file() and (str(dtm_file).endswith(".tif") or str(dtm_file).endswith(".TIF")): all_dtm_lhd_names.append(dtm_file.name) # bulk compute _ = Parallel(n_jobs=12, verbose=True)( - delayed(_compute_one_difference)(dtm_lhd_file, _dir, _out_dir_difference) for dtm_lhd_file in all_dtm_lhd_names + delayed(_compute_one_difference)(dtm_lhd_file, _dir, _out_dir_difference, second_dir=_second_dir) + for dtm_lhd_file in all_dtm_lhd_names ) @@ -45,6 +43,15 @@ def parse_args(): required=True, help="Dossier des dalles MNT Lidar HD", ) + + parser.add_argument( + "-s", + "--secondary_dtm_elevation_dir", + type=Path, + default=None, + help="Dossier contenant le second ensemble de dalles MNT pour le calcul de différence", + ) + parser.add_argument( "-o", "--name_dir_difference", @@ -58,4 +65,4 @@ def parse_args(): if __name__ == "__main__": args = parse_args() - compute_all_difference_maps(args.dtm_lidar_dir, args.name_dir_difference) + compute_all_difference_maps(args.dtm_lidar_dir, args.secondary_dtm_elevation_dir, args.name_dir_difference) diff --git a/test/test_compute_difference.py b/test/test_compute_difference.py index ffe303e..f4c041f 100644 --- a/test/test_compute_difference.py +++ b/test/test_compute_difference.py @@ -166,8 +166,20 @@ def test_main(): output_dir.mkdir() dtm_lidar_file = "./data/lhd/Semis_2021_0886_6443_LA93_IGN69_50CM.tif" + second_dtm_file = "./data/lhd/Semis_2021_0886_6443_LA93_IGN69_50CM.tif" + out_difference_file = output_dir / "Difference_with_rge_Semis_2021_0886_6443_LA93_IGN69_50CM.tif" + out_difference_file_with_self = output_dir / "Difference_self_with_rge_Semis_2021_0886_6443_LA93_IGN69_50CM.tif" - compute_difference.main(dtm_lidar_file, out_difference_file) + compute_difference.main(dtm_lidar_file, None, out_difference_file) assert os.path.exists(out_difference_file), "difference with rge alti not computed !" + + # run with secondary file + compute_difference.main(dtm_lidar_file, second_dtm_file, out_difference_file_with_self) + + assert os.path.exists(out_difference_file_with_self), "difference with self is not computed" + + 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 null values" diff --git a/test/test_run_difference_with_gpao.py b/test/test_run_difference_with_gpao.py index b0aa244..75804a2 100644 --- a/test/test_run_difference_with_gpao.py +++ b/test/test_run_difference_with_gpao.py @@ -28,7 +28,7 @@ def test_create_gpao_project(): output_dir = TMP_PATH / "create_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_gpao_project(dtm_lidar_lhds, output_dir, STORE, project_name) + project = run_difference_with_gpao.create_gpao_project(dtm_lidar_lhds, None, output_dir, STORE, project_name) assert project is not None @@ -53,7 +53,7 @@ def test_gpao_run(): local_store_path = Path("data/lhd_dir_gpao").resolve() run_difference_with_gpao.compute_on_gpao( - Path(dtm_lidar_lhds), Path(output_dir), gpao_hostname, local_store_path, runner_store_path, project_name + Path(dtm_lidar_lhds), None, Path(output_dir), gpao_hostname, local_store_path, runner_store_path, project_name ) if gpao_hostname == "localhost": diff --git a/test/test_run_difference_with_parallel.py b/test/test_run_difference_with_parallel.py index ff0332a..f1746ad 100644 --- a/test/test_run_difference_with_parallel.py +++ b/test/test_run_difference_with_parallel.py @@ -19,7 +19,8 @@ def test_compute_all_difference_maps(): output_dir = TMP_PATH / "compute_all_difference_maps" output_dir.mkdir() dtm_lidar_dir = Path("./data/lhd_dir_gpao") - run_difference_with_parallel.compute_all_difference_maps(dtm_lidar_dir, output_dir) + + run_difference_with_parallel.compute_all_difference_maps(dtm_lidar_dir, None, output_dir) # check all difference maps are computed and stored in output_dir for lidar_dtm_file in dtm_lidar_dir.iterdir(): @@ -29,3 +30,17 @@ def test_compute_all_difference_maps(): assert os.path.exists(_matching_difference_file), "Difference file of {} is not computed!".format( lidar_dtm_file ) + + # second test with secondary folder of dtm files holding the same names + output_dir_with_sec = TMP_PATH / "compute_all_difference_maps" + sec_dtm_dir = Path("./data/lhd_dir_gpao") + + run_difference_with_parallel.compute_all_difference_maps(dtm_lidar_dir, sec_dtm_dir, output_dir_with_sec) + + for lidar_dtm_file in dtm_lidar_dir.iterdir(): + if lidar_dtm_file.is_file() and str(lidar_dtm_file).endswith(".tif"): + # check if corresponding difference file exists + _matching_difference_file = output_dir_with_sec / lidar_dtm_file.name + assert os.path.exists(_matching_difference_file), "Difference file of {} is not computed!".format( + lidar_dtm_file + ) From 6cc4b82e648f7507efe93b9d465dc228a147ea9d Mon Sep 17 00:00:00 2001 From: DaliCHEBBI Date: Thu, 25 Sep 2025 09:22:49 +0200 Subject: [PATCH 2/6] rename test outputs --- test/test_run_difference_with_gpao.py | 33 +++++++++++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/test/test_run_difference_with_gpao.py b/test/test_run_difference_with_gpao.py index 75804a2..7bb6614 100644 --- a/test/test_run_difference_with_gpao.py +++ b/test/test_run_difference_with_gpao.py @@ -40,11 +40,11 @@ def test_create_gpao_project(): @pytest.mark.gpao -def test_gpao_run(): +def test_gpao_run_with_stream(): dtm_lidar_lhds = "./data/lhd_dir_gpao" output_dir = TMP_PATH / "gpao_run" output_dir.mkdir() - project_name = "test_run_altianalysis_gpao" + project_name = "test_run_altianalysis_gpao_stream" gpao_hostname = os.environ.get("GPAO_API_URL", "localhost") url_api = f"http://{gpao_hostname}:8080/api/" @@ -59,3 +59,32 @@ def test_gpao_run(): if gpao_hostname == "localhost": tu.execute_gpao_client(tags="docker", num_thread=4) wait_running_job(url_api, project_name, delay_second=1, delay_log_second=10) + + +@pytest.mark.gpao +def test_gpao_run_with_given_secondary_folder(): + dtm_lidar_lhds = "./data/lhd_dir_gpao" + secondary_dtm_dir = "./data/lhd_dir_gpao" + output_dir = TMP_PATH / "gpao_run" + output_dir.mkdir() + project_name = "test_run_altianalysis_gpao_secondary_folder" + + gpao_hostname = os.environ.get("GPAO_API_URL", "localhost") + url_api = f"http://{gpao_hostname}:8080/api/" + + runner_store_path = Path(dtm_lidar_lhds).resolve() + local_store_path = Path("data/lhd_dir_gpao").resolve() + + run_difference_with_gpao.compute_on_gpao( + Path(dtm_lidar_lhds), + Path(secondary_dtm_dir), + Path(output_dir), + gpao_hostname, + local_store_path, + runner_store_path, + project_name, + ) + + if gpao_hostname == "localhost": + tu.execute_gpao_client(tags="docker", num_thread=4) + wait_running_job(url_api, project_name, delay_second=1, delay_log_second=10) From 321353770f3aa02d70bd7bdca1ce093839d734c8 Mon Sep 17 00:00:00 2001 From: DaliCHEBBI Date: Thu, 25 Sep 2025 14:08:57 +0200 Subject: [PATCH 3/6] make some changes after review --- altianalysis/compute_difference.py | 14 ++++++++++---- altianalysis/run_difference_with_gpao.py | 18 +++++++++++------- altianalysis/run_difference_with_parallel.py | 19 ++++++++++--------- test/test_compute_difference.py | 19 +++++++++++++------ test/test_run_difference_with_gpao.py | 4 ++-- test/test_run_difference_with_parallel.py | 11 ++++++++--- 6 files changed, 54 insertions(+), 31 deletions(-) diff --git a/altianalysis/compute_difference.py b/altianalysis/compute_difference.py index 60e41ac..01b5ed3 100644 --- a/altianalysis/compute_difference.py +++ b/altianalysis/compute_difference.py @@ -11,10 +11,16 @@ def parse_args(): - parser = argparse.ArgumentParser("compute difference map between lidar dtm and rge alti") - parser.add_argument("--dtm_lidar_file", type=str, help="dtm lidar tif file") + parser = argparse.ArgumentParser( + "compute difference between two elevation files or between one elevation file and rge alti" + ) + parser.add_argument("--primary_elevation_file", type=str, help="dtm lidar tif file") parser.add_argument( - "--second_elevation_file", default=None, type=str, help="secondary elevation file when no stream to be used !" + "--second_elevation_file", + default=None, + type=str, + help="secondary elevation file to compare with primary elevation file. " + "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") return parser.parse_args() @@ -113,4 +119,4 @@ def main(reference_dtm_file: Path | str, secondary_dtm_file: Path | str | None, if __name__ == "__main__": args = parse_args() - main(args.dtm_lidar_file, args.second_elevation_file, args.name_save_out) + main(args.primary_elevation_file, args.second_elevation_file, args.name_save_out) diff --git a/altianalysis/run_difference_with_gpao.py b/altianalysis/run_difference_with_gpao.py index 62b4573..05ad338 100644 --- a/altianalysis/run_difference_with_gpao.py +++ b/altianalysis/run_difference_with_gpao.py @@ -1,5 +1,6 @@ import argparse import logging +import os from pathlib import Path, PurePosixPath from typing import List @@ -27,26 +28,28 @@ def get_tile_names(folder: Path) -> List[str]: def create_one_job_one_difference(store: Store, dir_in: Path, second_dir: Path | None, input_file: str, output: Path): job_name = f"difference_{input_file}" - if second_dir is not None: + if second_dir: + secondary_file = second_dir / input_file + assert os.path.exists(secondary_file), "Secondary elevation does not exist to compute difference !" command = f""" - docker run -t --rm --userns=host --shm-size=2gb + docker run -t --rm --userns=host -v {store.to_unix(dir_in)}:/input -v {store.to_unix(second_dir)}:/second_input -v {store.to_unix(output)}:/output ghcr.io/ignf/altianalysis:{__version__} python -m altianalysis.compute_difference - --dtm_lidar_file /input/{input_file} + --primary_elevation_file /input/{input_file} --second_elevation_file /second_input/{input_file} --name_save_out /output/{input_file} """ else: command = f""" - docker run -t --rm --userns=host --shm-size=2gb + docker run -t --rm --userns=host -v {store.to_unix(dir_in)}:/input -v {store.to_unix(output)}:/output ghcr.io/ignf/altianalysis:{__version__} python -m altianalysis.compute_difference - --dtm_lidar_file /input/{input_file} + --primary_elevation_file /input/{input_file} --name_save_out /output/{input_file} """ job = Job(job_name, command, tags=["docker"]) @@ -181,11 +184,12 @@ def parse_args(): help="Dossier des dalles MNT Lidar HD", ) parser.add_argument( - "-i", + "-ii", "--secondary_dtm_dir", type=Path, default=None, - help="Dossier du second des MNT pour calculer la différence", + help="Dossier contenant le second ensemble de dalles MNT pour le calcul de différence. " + "Les dalles d'élévation doivent avoir les mêmes noms pour faire l'appariement", ) parser.add_argument( diff --git a/altianalysis/run_difference_with_parallel.py b/altianalysis/run_difference_with_parallel.py index 8ba2df1..af25eb1 100644 --- a/altianalysis/run_difference_with_parallel.py +++ b/altianalysis/run_difference_with_parallel.py @@ -11,25 +11,25 @@ def _compute_one_difference(dtm_lhd_file: str, _dir: Path, _out_dir_difference: full_dtm_file = os.path.join(_dir, dtm_lhd_file) full_difference_file = os.path.join(_out_dir_difference, dtm_lhd_file) - _second_file = None + second_file = None - if second_dir is not None: - _second_file = os.path.join(second_dir, dtm_lhd_file) + if second_dir: + second_file = os.path.join(second_dir, dtm_lhd_file) - main(full_dtm_file, _second_file, full_difference_file) + main(full_dtm_file, second_file, full_difference_file) -def compute_all_difference_maps(_dir: Path, _second_dir: Path | None, _out_dir_difference: Path): +def compute_all_difference_maps(dir: Path, second_dir: Path | None, out_dir_difference: Path): - os.makedirs(_out_dir_difference, exist_ok=True) + os.makedirs(out_dir_difference, exist_ok=True) all_dtm_lhd_names = [] - for dtm_file in _dir.iterdir(): + for dtm_file in dir.iterdir(): if dtm_file.is_file() and (str(dtm_file).endswith(".tif") or str(dtm_file).endswith(".TIF")): all_dtm_lhd_names.append(dtm_file.name) # bulk compute _ = Parallel(n_jobs=12, verbose=True)( - delayed(_compute_one_difference)(dtm_lhd_file, _dir, _out_dir_difference, second_dir=_second_dir) + delayed(_compute_one_difference)(dtm_lhd_file, dir, out_dir_difference, second_dir=second_dir) for dtm_lhd_file in all_dtm_lhd_names ) @@ -49,7 +49,8 @@ def parse_args(): "--secondary_dtm_elevation_dir", type=Path, default=None, - help="Dossier contenant le second ensemble de dalles MNT pour le calcul de différence", + help="Dossier contenant le second ensemble de dalles MNT pour le calcul de différence. " + "Les dalles d'élévation doivent avoir les mêmes noms pour faire l'appariement", ) parser.add_argument( diff --git a/test/test_compute_difference.py b/test/test_compute_difference.py index f4c041f..3dc0c1c 100644 --- a/test/test_compute_difference.py +++ b/test/test_compute_difference.py @@ -161,20 +161,27 @@ def test_compute_difference_with_nodata(): assert np.all(nodata_diff_mask[:, :-1, :-1] == nodata_combined), "Added or missed nodata values !" -def test_main(): - output_dir = TMP_PATH / "main" +def test_main_with_stream(): + output_dir = TMP_PATH / "main_with_stream" output_dir.mkdir() dtm_lidar_file = "./data/lhd/Semis_2021_0886_6443_LA93_IGN69_50CM.tif" - second_dtm_file = "./data/lhd/Semis_2021_0886_6443_LA93_IGN69_50CM.tif" - out_difference_file = output_dir / "Difference_with_rge_Semis_2021_0886_6443_LA93_IGN69_50CM.tif" - out_difference_file_with_self = output_dir / "Difference_self_with_rge_Semis_2021_0886_6443_LA93_IGN69_50CM.tif" compute_difference.main(dtm_lidar_file, None, out_difference_file) assert os.path.exists(out_difference_file), "difference with rge alti not computed !" + +def test_main_with_secondary_folder(): + output_dir = TMP_PATH / "main_with_secondary_folder" + output_dir.mkdir() + dtm_lidar_file = "./data/lhd/Semis_2021_0886_6443_LA93_IGN69_50CM.tif" + + second_dtm_file = "./data/lhd/Semis_2021_0886_6443_LA93_IGN69_50CM.tif" + + out_difference_file_with_self = output_dir / "Difference_self_with_rge_Semis_2021_0886_6443_LA93_IGN69_50CM.tif" + # run with secondary file compute_difference.main(dtm_lidar_file, second_dtm_file, out_difference_file_with_self) @@ -182,4 +189,4 @@ def test_main(): 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 null values" + assert np.all(diff == 0), "difference with self yields non null values" diff --git a/test/test_run_difference_with_gpao.py b/test/test_run_difference_with_gpao.py index a6eac39..cca2741 100644 --- a/test/test_run_difference_with_gpao.py +++ b/test/test_run_difference_with_gpao.py @@ -65,7 +65,7 @@ def test_create_gpao_projects_without_cog(): @pytest.mark.gpao def test_gpao_run_with_cog_stream_rge(): dtm_lidar_lhds = "./data/lhd_dir_gpao" - output_dir = TMP_PATH / "gpao_run_with_cog" + output_dir = TMP_PATH / "gpao_run_with_cog_stream_rge" output_dir.mkdir() cog_filename = "cog.tif" project_name = "test_run_altianalysis_gpao" @@ -124,7 +124,7 @@ def test_gpao_run_without_cog_stream_rge(): def test_gpao_run_without_cog_given_secondary_folder(): dtm_lidar_lhds = "./data/lhd_dir_gpao" secondary_dtm_dir = "./data/lhd_dir_gpao" - output_dir = TMP_PATH / "gpao_run" + output_dir = TMP_PATH / "gpao_run_without_cog_given_secondary_folder" output_dir.mkdir() project_name = "test_run_altianalysis_gpao_secondary_folder" diff --git a/test/test_run_difference_with_parallel.py b/test/test_run_difference_with_parallel.py index f1746ad..bf2ecaa 100644 --- a/test/test_run_difference_with_parallel.py +++ b/test/test_run_difference_with_parallel.py @@ -15,8 +15,8 @@ def setup_module(module): os.makedirs(TMP_PATH) -def test_compute_all_difference_maps(): - output_dir = TMP_PATH / "compute_all_difference_maps" +def test_compute_all_difference_maps_with_stream(): + output_dir = TMP_PATH / "compute_all_difference_maps_with_stream" output_dir.mkdir() dtm_lidar_dir = Path("./data/lhd_dir_gpao") @@ -31,8 +31,13 @@ def test_compute_all_difference_maps(): lidar_dtm_file ) + +def test_compute_all_difference_maps_with_secondary_folder(): # second test with secondary folder of dtm files holding the same names - output_dir_with_sec = TMP_PATH / "compute_all_difference_maps" + output_dir_with_sec = TMP_PATH / "compute_all_difference_maps_with_secondary_folder" + output_dir_with_sec.mkdir() + + dtm_lidar_dir = Path("./data/lhd_dir_gpao") sec_dtm_dir = Path("./data/lhd_dir_gpao") run_difference_with_parallel.compute_all_difference_maps(dtm_lidar_dir, sec_dtm_dir, output_dir_with_sec) From d02b8d47dd96639f790f30523fddf73ee0ce6edd Mon Sep 17 00:00:00 2001 From: DaliCHEBBI Date: Thu, 25 Sep 2025 17:27:36 +0200 Subject: [PATCH 4/6] generic naming of input files --- altianalysis/compute_difference.py | 32 ++++++++++---------- altianalysis/run_difference_with_gpao.py | 28 +++++++++-------- altianalysis/run_difference_with_parallel.py | 22 +++++++------- 3 files changed, 42 insertions(+), 40 deletions(-) diff --git a/altianalysis/compute_difference.py b/altianalysis/compute_difference.py index 01b5ed3..76b9d9d 100644 --- a/altianalysis/compute_difference.py +++ b/altianalysis/compute_difference.py @@ -14,7 +14,7 @@ def parse_args(): parser = argparse.ArgumentParser( "compute difference between two elevation files or between one elevation file and rge alti" ) - parser.add_argument("--primary_elevation_file", type=str, help="dtm lidar tif file") + parser.add_argument("--primary_elevation_file", type=str, help="primary elevation file") parser.add_argument( "--second_elevation_file", default=None, @@ -27,7 +27,7 @@ def parse_args(): def _extract_rge_alti_tile_from_stream( - dtm_lidar_file: str, + dtm_file: str, output_path: str | Path, stream_RGE="RGEALTI-MNT_PYR-ZIP_FXX_LAMB93_WMS", # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES", proj="2154", @@ -37,7 +37,7 @@ def _extract_rge_alti_tile_from_stream( ): # compute dtm lidar extent while reading with rasterio - with rasterio.open(dtm_lidar_file) as dtm_lhd: + with rasterio.open(dtm_file) as dtm_lhd: bounds = dtm_lhd.bounds minx, miny, maxx, maxy = bounds.left, bounds.bottom, bounds.right, bounds.top @@ -58,30 +58,30 @@ def _extract_rge_alti_tile_from_stream( return True except requests.exceptions.Timeout: - warnings.warn(f"Request time out for stream associated to {dtm_lidar_file}") + warnings.warn(f"Request time out for stream associated to {dtm_file}") return False -def compute_difference_between_dtms(dtm_lidar_file: str, dtm_rge_tile: str, name_save_out: str) -> None: +def compute_difference_between_dtms(first_dtm_file: str, second_dtm_file: str, name_save_out: str) -> None: # read dtm_lidar tile - with rasterio.open(dtm_lidar_file) as dtm_lidar, rasterio.open(dtm_rge_tile) as dtm_rge: + with rasterio.open(first_dtm_file) as dtm_1, rasterio.open(second_dtm_file) as dtm_2: # read dtm_lidar array - data_dtm_lidar = dtm_lidar.read(1) + data_dtm_1 = dtm_1.read(1) # copy metadata of mtd_lidar - meta_dtm_lidar = dtm_lidar.meta.copy() + meta_dtm_1 = dtm_1.meta.copy() # window from rge alti - window_rge = dtm_rge.window(*dtm_lidar.bounds) + window_2 = dtm_2.window(*dtm_1.bounds) # read the data from dtm_rge with the same window as dtm_lidar - data_dtm_rge_windowed = dtm_rge.read( + data_dtm_2_windowed = dtm_2.read( 1, - window=window_rge, + window=window_2, boundless=True, - out_shape=data_dtm_lidar.shape, + out_shape=data_dtm_1.shape, resampling=Resampling.cubic, fill_value=0, ) @@ -89,14 +89,14 @@ def compute_difference_between_dtms(dtm_lidar_file: str, dtm_rge_tile: str, name # data_dtm_rge_windowed=np.where(data_dtm_rge_windowed==dtm_rge.nodata, 0,data_dtm_rge_windowed) # difference dem - _difference = data_dtm_lidar - data_dtm_rge_windowed + _difference = data_dtm_1 - data_dtm_2_windowed - _masq_nodata = np.logical_or((data_dtm_rge_windowed == dtm_rge.nodata), (data_dtm_lidar == dtm_lidar.nodata)) + _masq_nodata = np.logical_or((data_dtm_2_windowed == dtm_2.nodata), (data_dtm_1 == dtm_1.nodata)) - _difference[_masq_nodata] = dtm_lidar.nodata + _difference[_masq_nodata] = dtm_1.nodata # save result - with rasterio.open(name_save_out, "w", **meta_dtm_lidar) as dst: + with rasterio.open(name_save_out, "w", **meta_dtm_1) as dst: dst.write(_difference, 1) diff --git a/altianalysis/run_difference_with_gpao.py b/altianalysis/run_difference_with_gpao.py index 05ad338..00b4bc1 100644 --- a/altianalysis/run_difference_with_gpao.py +++ b/altianalysis/run_difference_with_gpao.py @@ -30,7 +30,10 @@ def create_one_job_one_difference(store: Store, dir_in: Path, second_dir: Path | job_name = f"difference_{input_file}" if second_dir: secondary_file = second_dir / input_file - assert os.path.exists(secondary_file), "Secondary elevation does not exist to compute difference !" + + if not os.path.exists(secondary_file): + raise FileNotFoundError(f"Secondary elevation does not exist to compute difference for {input_file}") + command = f""" docker run -t --rm --userns=host -v {store.to_unix(dir_in)}:/input @@ -57,21 +60,21 @@ def create_one_job_one_difference(store: Store, dir_in: Path, second_dir: Path | def create_main_gpao_project( - dtms_lhd: Path, + primary_dir: Path, secondary_dir: Path | None, out: Path, store: Store, project_name: str, ) -> Project: - logging.debug(f"Create GPAO projects to compute difference maps with rge alti: {dtms_lhd}.") + logging.debug(f"Create GPAO projects to compute difference maps with rge alti: {primary_dir}.") logging.debug(f"Writing difference maps to {out}.") Project.reset() # get dtm tile names for lidar hd - dtm_tile_names = get_tile_names(dtms_lhd) + dtm_tile_names = get_tile_names(primary_dir) out.mkdir(parents=True, exist_ok=True) @@ -80,7 +83,7 @@ def create_main_gpao_project( jobs = [] for tile_ in dtm_tile_names: - job = create_one_job_one_difference(store, dtms_lhd, secondary_dir, tile_, out) + job = create_one_job_one_difference(store, primary_dir, secondary_dir, tile_, out) jobs.append(job) # create the project @@ -120,9 +123,9 @@ def create_cog_gpao_project( def create_gpao_projects( - dtms_lhd: 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 ) -> List[Project]: - project_main = create_main_gpao_project(dtms_lhd, secondary_dir, out, store, project_name) + project_main = create_main_gpao_project(primary_dir, secondary_dir, out, store, project_name) projects = [project_main] if cog_filename: project_cog = create_cog_gpao_project( @@ -134,7 +137,7 @@ def create_gpao_projects( def compute_on_gpao( - dtms_lhd: Path, + primary_dir: Path, secondary_dir: Path | None, out: Path, gpao_hostname: str, @@ -163,7 +166,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(dtms_lhd, secondary_dir, out, store, project_name, cog_filename) + projects = create_gpao_projects(primary_dir, secondary_dir, out, store, project_name, cog_filename) builder = Builder(projects) @@ -177,11 +180,10 @@ def parse_args(): parser = argparse.ArgumentParser(description="Calcul de cartes de différences par rapport au RGE ALTI") parser.add_argument( "-i", - "--dtm_lhd_dir", + "--primary_dtm_dir", type=Path, - # nargs="+", required=True, - help="Dossier des dalles MNT Lidar HD", + help="Dossier contenant le premier ensemble de dalles MNT pour le calcul de différence.", ) parser.add_argument( "-ii", @@ -234,7 +236,7 @@ def parse_args(): args = parse_args() compute_on_gpao( - dtms_lhd=args.dtm_lhd_dir, + dtms_lhd=args.primary_dtm_dir, secondary_dir=args.secondary_dtm_dir, out=args.out, gpao_hostname=args.gpao_hostname, diff --git a/altianalysis/run_difference_with_parallel.py b/altianalysis/run_difference_with_parallel.py index af25eb1..386dd49 100644 --- a/altianalysis/run_difference_with_parallel.py +++ b/altianalysis/run_difference_with_parallel.py @@ -7,14 +7,14 @@ from altianalysis.compute_difference import main -def _compute_one_difference(dtm_lhd_file: str, _dir: Path, _out_dir_difference: Path, second_dir=None): +def _compute_one_difference(dtm_file: str, _dir: Path, _out_dir_difference: Path, second_dir=None): - full_dtm_file = os.path.join(_dir, dtm_lhd_file) - full_difference_file = os.path.join(_out_dir_difference, dtm_lhd_file) + full_dtm_file = os.path.join(_dir, dtm_file) + full_difference_file = os.path.join(_out_dir_difference, dtm_file) second_file = None if second_dir: - second_file = os.path.join(second_dir, dtm_lhd_file) + second_file = os.path.join(second_dir, dtm_file) main(full_dtm_file, second_file, full_difference_file) @@ -22,15 +22,15 @@ def _compute_one_difference(dtm_lhd_file: str, _dir: Path, _out_dir_difference: def compute_all_difference_maps(dir: Path, second_dir: Path | None, out_dir_difference: Path): os.makedirs(out_dir_difference, exist_ok=True) - all_dtm_lhd_names = [] + all_dtm_names = [] for dtm_file in dir.iterdir(): if dtm_file.is_file() and (str(dtm_file).endswith(".tif") or str(dtm_file).endswith(".TIF")): - all_dtm_lhd_names.append(dtm_file.name) + all_dtm_names.append(dtm_file.name) # bulk compute _ = Parallel(n_jobs=12, verbose=True)( - delayed(_compute_one_difference)(dtm_lhd_file, dir, out_dir_difference, second_dir=second_dir) - for dtm_lhd_file in all_dtm_lhd_names + delayed(_compute_one_difference)(dtm_file, dir, out_dir_difference, second_dir=second_dir) + for dtm_file in all_dtm_names ) @@ -38,10 +38,10 @@ def parse_args(): parser = argparse.ArgumentParser(description="Calcul de cartes de différences par rapport au RGE ALTI") parser.add_argument( "-l", - "--dtm_lidar_dir", + "--primary_dtm_dir", type=Path, required=True, - help="Dossier des dalles MNT Lidar HD", + help="Dossier contenant le premier ensemble de dalles MNT pour le calcul de différence.", ) parser.add_argument( @@ -66,4 +66,4 @@ def parse_args(): if __name__ == "__main__": args = parse_args() - compute_all_difference_maps(args.dtm_lidar_dir, args.secondary_dtm_elevation_dir, args.name_dir_difference) + compute_all_difference_maps(args.primary_dtm_dir, args.secondary_dtm_elevation_dir, args.name_dir_difference) From 1c818dca3ae3f3d954dbd5394c5b743b2f5e5e19 Mon Sep 17 00:00:00 2001 From: DaliCHEBBI Date: Mon, 29 Sep 2025 10:15:36 +0200 Subject: [PATCH 5/6] compile change requests --- altianalysis/compute_difference.py | 4 +- altianalysis/run_difference_with_gpao.py | 41 ++++++++++++-------- altianalysis/run_difference_with_parallel.py | 7 ++-- 3 files changed, 31 insertions(+), 21 deletions(-) diff --git a/altianalysis/compute_difference.py b/altianalysis/compute_difference.py index 76b9d9d..ceccab4 100644 --- a/altianalysis/compute_difference.py +++ b/altianalysis/compute_difference.py @@ -37,8 +37,8 @@ def _extract_rge_alti_tile_from_stream( ): # compute dtm lidar extent while reading with rasterio - with rasterio.open(dtm_file) as dtm_lhd: - bounds = dtm_lhd.bounds + with rasterio.open(dtm_file) as dtm: + bounds = dtm.bounds minx, miny, maxx, maxy = bounds.left, bounds.bottom, bounds.right, bounds.top try: diff --git a/altianalysis/run_difference_with_gpao.py b/altianalysis/run_difference_with_gpao.py index 00b4bc1..1074187 100644 --- a/altianalysis/run_difference_with_gpao.py +++ b/altianalysis/run_difference_with_gpao.py @@ -34,27 +34,26 @@ def create_one_job_one_difference(store: Store, dir_in: Path, second_dir: Path | if not os.path.exists(secondary_file): raise FileNotFoundError(f"Secondary elevation does not exist to compute difference for {input_file}") - command = f""" - docker run -t --rm --userns=host - -v {store.to_unix(dir_in)}:/input - -v {store.to_unix(second_dir)}:/second_input - -v {store.to_unix(output)}:/output - ghcr.io/ignf/altianalysis:{__version__} - python -m altianalysis.compute_difference - --primary_elevation_file /input/{input_file} - --second_elevation_file /second_input/{input_file} - --name_save_out /output/{input_file} - """ + second_input_mount = f" -v {store.to_unix(second_dir)}:/second_input" + second_input_param = f"--second_elevation_file /second_input/{input_file}" + else: - command = f""" + + second_input_mount = "" + second_input_param = "" + + command = f""" docker run -t --rm --userns=host -v {store.to_unix(dir_in)}:/input + {second_input_mount} -v {store.to_unix(output)}:/output ghcr.io/ignf/altianalysis:{__version__} python -m altianalysis.compute_difference --primary_elevation_file /input/{input_file} + {second_input_param} --name_save_out /output/{input_file} """ + job = Job(job_name, command, tags=["docker"]) return job @@ -67,7 +66,13 @@ def create_main_gpao_project( project_name: str, ) -> Project: - logging.debug(f"Create GPAO projects to compute difference maps with rge alti: {primary_dir}.") + if secondary_dir: + logging.debug( + f"Create GPAO projects to compute difference maps file by file between : \ + {primary_dir} and {secondary_dir}.Note that files are matched by names." + ) + else: + logging.debug(f"Create GPAO projects to compute difference maps with rge alti: {primary_dir}.") logging.debug(f"Writing difference maps to {out}.") Project.reset() @@ -177,7 +182,10 @@ def compute_on_gpao( def parse_args(): - parser = argparse.ArgumentParser(description="Calcul de cartes de différences par rapport au RGE ALTI") + parser = argparse.ArgumentParser( + description="Calcul de cartes de différences entre 2 sources de modèles numériques, " + "ou d'une source par rapport au RGE ALTI" + ) parser.add_argument( "-i", "--primary_dtm_dir", @@ -191,7 +199,8 @@ def parse_args(): type=Path, default=None, help="Dossier contenant le second ensemble de dalles MNT pour le calcul de différence. " - "Les dalles d'élévation doivent avoir les mêmes noms pour faire l'appariement", + "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( @@ -236,7 +245,7 @@ def parse_args(): args = parse_args() compute_on_gpao( - dtms_lhd=args.primary_dtm_dir, + primary_dir=args.primary_dtm_dir, secondary_dir=args.secondary_dtm_dir, out=args.out, gpao_hostname=args.gpao_hostname, diff --git a/altianalysis/run_difference_with_parallel.py b/altianalysis/run_difference_with_parallel.py index 386dd49..fe0a0ef 100644 --- a/altianalysis/run_difference_with_parallel.py +++ b/altianalysis/run_difference_with_parallel.py @@ -46,11 +46,12 @@ def parse_args(): parser.add_argument( "-s", - "--secondary_dtm_elevation_dir", + "--secondary_dtm_dir", type=Path, default=None, help="Dossier contenant le second ensemble de dalles MNT pour le calcul de différence. " - "Les dalles d'élévation doivent avoir les mêmes noms pour faire l'appariement", + "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( @@ -66,4 +67,4 @@ def parse_args(): if __name__ == "__main__": args = parse_args() - compute_all_difference_maps(args.primary_dtm_dir, args.secondary_dtm_elevation_dir, args.name_dir_difference) + compute_all_difference_maps(args.primary_dtm_dir, args.secondary_dtm_dir, args.name_dir_difference) From 2d5a1e89e5d4a988159cecb9e60d1f6f9cb07afb Mon Sep 17 00:00:00 2001 From: DaliCHEBBI Date: Mon, 29 Sep 2025 11:27:22 +0200 Subject: [PATCH 6/6] compile change requests --- altianalysis/run_difference_with_gpao.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/altianalysis/run_difference_with_gpao.py b/altianalysis/run_difference_with_gpao.py index 1074187..08902c9 100644 --- a/altianalysis/run_difference_with_gpao.py +++ b/altianalysis/run_difference_with_gpao.py @@ -68,7 +68,7 @@ def create_main_gpao_project( if secondary_dir: logging.debug( - f"Create GPAO projects to compute difference maps file by file between : \ + f"Create GPAO projects to compute difference maps file by file between: \ {primary_dir} and {secondary_dir}.Note that files are matched by names." ) else: