diff --git a/CHANGELOG.md b/CHANGELOG.md index 84b6642..280cc1b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,5 @@ # master - +-feature: handle two folders of DEMs for difference computation -feature: option to generate a COG when running with gpao # 1.1.0 diff --git a/altianalysis/compute_difference.py b/altianalysis/compute_difference.py index d4c76ad..ceccab4 100644 --- a/altianalysis/compute_difference.py +++ b/altianalysis/compute_difference.py @@ -11,14 +11,23 @@ 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="primary elevation file") + parser.add_argument( + "--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() 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", @@ -28,8 +37,8 @@ def _extract_rge_alti_tile_from_stream( ): # compute dtm lidar extent while reading with rasterio - with rasterio.open(dtm_lidar_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: @@ -49,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, ) @@ -80,27 +89,34 @@ 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) -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.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 be519a8..08902c9 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 @@ -25,36 +26,60 @@ 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}" + if second_dir: + secondary_file = second_dir / input_file + + if not os.path.exists(secondary_file): + raise FileNotFoundError(f"Secondary elevation does not exist to compute difference for {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: + + second_input_mount = "" + second_input_param = "" + command = f""" -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} ---name_save_out /output/{input_file} -""" + 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 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}.") + 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() # 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) @@ -63,7 +88,7 @@ def create_main_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, primary_dir, secondary_dir, tile_, out) jobs.append(job) # create the project @@ -103,9 +128,9 @@ def create_cog_gpao_project( def create_gpao_projects( - dtms_lhd: Path, 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, 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( @@ -117,7 +142,8 @@ def create_gpao_projects( def compute_on_gpao( - dtms_lhd: Path, + primary_dir: Path, + secondary_dir: Path | None, out: Path, gpao_hostname: str, local_store_path: Path, @@ -145,7 +171,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, out, store, project_name, cog_filename) + projects = create_gpao_projects(primary_dir, secondary_dir, out, store, project_name, cog_filename) builder = Builder(projects) @@ -156,15 +182,27 @@ 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", - "--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", + "--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" + "S'il est laissé vide, les dalles du premier ensemble sont comparées au RGE Alti", + ) + parser.add_argument( "-o", "--out", @@ -207,7 +245,8 @@ def parse_args(): args = parse_args() compute_on_gpao( - dtms_lhd=args.dtm_lhd_dir, + primary_dir=args.primary_dtm_dir, + secondary_dir=args.secondary_dtm_dir, out=args.out, gpao_hostname=args.gpao_hostname, local_store_path=args.local_store_path, diff --git a/altianalysis/run_difference_with_parallel.py b/altianalysis/run_difference_with_parallel.py index db76947..fe0a0ef 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_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 - # compute difference for individual files - with tempfile.NamedTemporaryFile(suffix="_dtm_rgealti.tif", delete=False) as tmp_rge: + if second_dir: + second_file = os.path.join(second_dir, dtm_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(): + os.makedirs(out_dir_difference, exist_ok=True) + 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) 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 ) @@ -40,11 +38,22 @@ 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( + "-s", + "--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" + "S'il est laissé vide, les dalles du premier ensemble sont comparées au RGE Alti", ) + parser.add_argument( "-o", "--name_dir_difference", @@ -58,4 +67,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.primary_dtm_dir, args.secondary_dtm_dir, args.name_dir_difference) diff --git a/test/test_compute_difference.py b/test/test_compute_difference.py index ffe303e..3dc0c1c 100644 --- a/test/test_compute_difference.py +++ b/test/test_compute_difference.py @@ -161,13 +161,32 @@ 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" out_difference_file = output_dir / "Difference_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 !" + + +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) + + 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 non null values" diff --git a/test/test_run_difference_with_gpao.py b/test/test_run_difference_with_gpao.py index 9850eb1..cca2741 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, output_dir, STORE, project_name) + project = run_difference_with_gpao.create_main_gpao_project(dtm_lidar_lhds, None, output_dir, STORE, project_name) 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, output_dir, STORE, project_name, cog_filename + dtm_lidar_lhds, None, output_dir, STORE, project_name, cog_filename ) assert len(projects) == 2 @@ -56,16 +56,16 @@ 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, output_dir, STORE, project_name, "") + projects = run_difference_with_gpao.create_gpao_projects(dtm_lidar_lhds, None, output_dir, STORE, project_name, "") assert len(projects) == 1 assert len(projects[0].jobs) == 5 @pytest.mark.gpao -def test_gpao_run_with_cog(): +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" @@ -78,6 +78,7 @@ def test_gpao_run_with_cog(): run_difference_with_gpao.compute_on_gpao( Path(dtm_lidar_lhds), + None, Path(output_dir), gpao_hostname, local_store_path, @@ -98,11 +99,34 @@ def test_gpao_run_with_cog(): @pytest.mark.gpao -def test_gpao_run_without_cog(): +def test_gpao_run_without_cog_stream_rge(): dtm_lidar_lhds = "./data/lhd_dir_gpao" - output_dir = TMP_PATH / "gpao_run_without_cog" + output_dir = TMP_PATH / "gpao_run_without_cog_stream_rge" 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/" + + 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), None, 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) + + +@pytest.mark.gpao +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_without_cog_given_secondary_folder" + 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/" @@ -112,6 +136,7 @@ def test_gpao_run_without_cog(): run_difference_with_gpao.compute_on_gpao( Path(dtm_lidar_lhds), + Path(secondary_dtm_dir), Path(output_dir), gpao_hostname, local_store_path, diff --git a/test/test_run_difference_with_parallel.py b/test/test_run_difference_with_parallel.py index ff0332a..bf2ecaa 100644 --- a/test/test_run_difference_with_parallel.py +++ b/test/test_run_difference_with_parallel.py @@ -15,11 +15,12 @@ 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") - 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,22 @@ def test_compute_all_difference_maps(): assert os.path.exists(_matching_difference_file), "Difference file of {} is not computed!".format( 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_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) + + 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 + )