Feat/diff two folders

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

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,58 +58,65 @@ 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,
         )
 
         # 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)

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,

altianalysis/run_difference_with_parallel.py

@@ -1,50 +1,59 @@
 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
     )
 
 
 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)