feat: add option to compare dtm with lidarhd stream

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)

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,

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 !"

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,