Merge pull request #56 from IGNF/malt0-on-intersection

Malt0 on intersection

CHANGELOG.md

@@ -1,3 +1,11 @@
+# dev
+- Changement de la définition de la zone où est calculée MALT0:
+  - masquage du MNx à l'aide de la carte de classe dans MALT0 intrinsic
+  - comparaison des MNx sur l'intersection des zones hors no-data (pour exclure les pixels en bord de maillage :
+  inclus dans la carte de classe mais hors du maillage utilisé pour le MNx)
+- Disparition de la distinction entre calcul de la référence et l'autre nuage de points (était nécessaire uniquement
+pour le calcul de MALT0)
+
 # 0.6.0
 - ajout d'une 4e métrique : MOBJ0
 

coclico/main.py

@@ -148,7 +148,7 @@ def create_compare_project(
 
                 out_ref_metric = out_ref / metric_name / "intrinsic"
                 out_ref_metric.mkdir(parents=True, exist_ok=True)
-                metric_jobs = metric.create_metric_intrinsic_jobs("ref", tile_names, ref, out_ref_metric, is_ref=True)
+                metric_jobs = metric.create_metric_intrinsic_jobs("ref", tile_names, ref, out_ref_metric)
                 add_dependency_to_jobs(metric_jobs, unlock_job)
 
                 ref_jobs[metric_name] = metric_jobs
@@ -176,9 +176,7 @@ def create_compare_project(
                     out_ci_metric = out_ci / metric_name / "intrinsic"
 
                     out_ci_metric.mkdir(parents=True, exist_ok=True)
-                    ci_intrinsic_jobs = metric.create_metric_intrinsic_jobs(
-                        ci.name, tile_names, ci, out_ci_metric, is_ref=False
-                    )
+                    ci_intrinsic_jobs = metric.create_metric_intrinsic_jobs(ci.name, tile_names, ci, out_ci_metric)
                     add_dependency_to_jobs(ci_intrinsic_jobs, unlock_job)
 
                     out_ci_to_ref_metric = out_ci / metric_name / "to_ref"

coclico/malt0/malt0.py

@@ -22,16 +22,8 @@ class MALT0(Metric):
     pixel_size = 0.5
     metric_name = "malt0"
 
-    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path, is_ref: bool):
+    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path):
         job_name = f"{self.metric_name}_intrinsic_{name}_{input.stem}"
-        occupancy_map_arg = ""
-        mnx_out = output / "mnx"
-        mnx_out.mkdir(exist_ok=True, parents=True)
-        if is_ref:
-            occ_out = output / "occupancy"
-            occ_out.mkdir(exist_ok=True, parents=True)
-            occupancy_map_arg = f"--output-occupancy-file /output/occupancy/{input.stem}.tif"
-
         command = f"""
 docker run -t --rm --userns=host --shm-size=2gb
 -v {self.store.to_unix(input)}:/input
@@ -40,8 +32,7 @@ def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path,
 ignimagelidar/coclico:{__version__}
 python -m coclico.malt0.malt0_intrinsic
 --input-file /input
---output-mnx-file /output/mnx/{input.stem}.tif
-{occupancy_map_arg}
+--output-mnx-file /output/{input.stem}.tif
 --config-file /config/{self.config_file.name}
 --pixel-size {self.pixel_size}
 
@@ -56,16 +47,14 @@ def create_metric_relative_to_ref_jobs(
         job_name = f"{self.metric_name}_{name}_relative_to_ref"
         command = f"""
 docker run -t --rm --userns=host --shm-size=2gb
--v {self.store.to_unix(out_c1) / "mnx"}:/input
--v {self.store.to_unix(out_ref) / "mnx"}:/ref
--v {self.store.to_unix(out_ref) / "occupancy"}:/occupancy
+-v {self.store.to_unix(out_c1)}:/input
+-v {self.store.to_unix(out_ref)}:/ref
 -v {self.store.to_unix(output)}:/output
 -v {self.store.to_unix(self.config_file.parent)}:/config
 ignimagelidar/coclico:{__version__}
 python -m coclico.malt0.malt0_relative
 --input-dir /input
 --ref-dir /ref
---occupancy-dir /occupancy
 --output-csv-tile /output/result_tile.csv
 --output-csv /output/result.csv
 --config-file /config/{self.config_file.name}

coclico/malt0/malt0_intrinsic.py

@@ -1,8 +1,11 @@
 import argparse
 import logging
+import tempfile
 from pathlib import Path
 
+import numpy as np
 import pdal
+import rasterio
 
 import coclico.io
 import coclico.metrics.commons as commons
@@ -57,11 +60,32 @@ def create_mnx_map(las_file, class_weights, output_tif, pixel_size, no_data_valu
     pipeline.execute()
 
 
+def mask_raster_with_nodata(raster_to_mask: Path, mask: np.array, output_tif: Path):
+    """Replace data in raster_to_mask by "no_data_value" where mask=0
+    (usual masks are occupancy maps for which 1 indicates that there are data)
+
+    Args:
+        raster_to_mask (Path): Path to the raster on which to replace data
+        mask (np.array): Binary Array (with the same shape as raster_to_mask data) with 1 where to keep data,
+        0 where to replace them with the nodata value
+        output_tif (Path): Path to the resulting file
+    """
+    with rasterio.Env():
+        with rasterio.open(str(raster_to_mask)) as src:
+            raster = src.read()
+            out_meta = src.meta
+            nodata = src.nodata
+
+        raster[mask == 0] = nodata
+
+        with rasterio.open(str(output_tif), "w", **out_meta) as dest:
+            dest.write(raster)
+
+
 def compute_metric_intrinsic(
     las_file: Path,
     config_file: Path,
     output_tif: Path,
-    occupancy_tif: Path = None,
     pixel_size: float = 0.5,
     no_data_value=-9999,
 ):
@@ -76,33 +100,26 @@ def compute_metric_intrinsic(
         las_file (Path): path to the las file on which to generate malt0 intrinsic metric
         config_file (Path): class weights dict in the config file (to know for which classes to generate the rasters)
         output_tif (Path): path to output height raster
-        occupancy_tif (Path, optional): path to the optional occupancy map tif (Leave it to None except for the
-        reference folder). Defaults to None.
         pixel_size (float, optional): size of the output rasters pixels. Defaults to 0.5.
         no_data_value (int, optional): no_data value for the output raster. Defaults to -9999.
     """
     config_dict = coclico.io.read_config_file(config_file)
     class_weights = config_dict[MALT0.metric_name]["weights"]
 
-    if occupancy_tif:
-        occupancy_tif.parent.mkdir(parents=True, exist_ok=True)
-        occupancy_map.create_occupancy_map(las_file, class_weights, occupancy_tif, pixel_size)
-
     output_tif.parent.mkdir(parents=True, exist_ok=True)
-    create_mnx_map(las_file, class_weights, output_tif, pixel_size, no_data_value)
+
+    with tempfile.NamedTemporaryFile(suffix=las_file.stem + "_mnx.tif") as tmp_mnx:
+        xs, ys, classifs, _ = occupancy_map.read_las(las_file)
+        binary_maps, _, _ = occupancy_map.create_occupancy_map_array(xs, ys, classifs, pixel_size, class_weights)
+
+        create_mnx_map(las_file, class_weights, tmp_mnx.name, pixel_size, no_data_value)
+        mask_raster_with_nodata(tmp_mnx.name, binary_maps, output_tif)
 
 
 def parse_args():
     parser = argparse.ArgumentParser("Run malt0 intrinsic metric on one tile")
     parser.add_argument("-i", "--input-file", type=Path, required=True, help="Path to the LAS file")
     parser.add_argument("-o", "--output-mnx-file", type=Path, required=True, help="Path to the TIF output file")
-    parser.add_argument(
-        "-oc",
-        "--output-occupancy-file",
-        type=Path,
-        default=None,
-        help="Path to the optional occupancy map TIF output file",
-    )
     parser.add_argument(
         "-c",
         "--config-file",
@@ -121,5 +138,4 @@ def parse_args():
         las_file=Path(args.input_file),
         config_file=args.config_file,
         output_tif=Path(args.output_mnx_file),
-        occupancy_tif=Path(args.output_occupancy_file) if args.output_occupancy_file else None,
     )

coclico/malt0/malt0_relative.py

@@ -12,7 +12,7 @@
 from coclico.malt0.malt0 import MALT0
 
 
-def compute_stats_single_raster(raster: np.array, occupancy_raster: np.array):
+def compute_stats_single_raster(raster: np.array):
     """Compute stats for a single raster, masked by an occupancy raster.
 
     Returns a np.array for each statistic:
@@ -32,12 +32,11 @@ def compute_stats_single_raster(raster: np.array, occupancy_raster: np.array):
     Returns:
         np.arrays: max value, pixel count, mean, standard deviation and m2 values
     """
-    masked_raster = ma.masked_array(raster, 1 - occupancy_raster)
-    max_val = ma.max(masked_raster, axis=(1, 2)).filled(0)
-    count = np.sum(occupancy_raster, axis=(1, 2))
-    mean_val = ma.mean(masked_raster, axis=(1, 2)).filled(0)
-    std_val = ma.std(masked_raster, axis=(1, 2)).filled(0)
-    m2 = ma.sum(np.square(masked_raster - mean_val[:, None, None]), axis=(1, 2)).filled(0)  # distance to the mean
+    max_val = ma.max(raster, axis=(1, 2)).filled(0)
+    count = ma.count(raster, axis=(1, 2))
+    mean_val = ma.mean(raster, axis=(1, 2)).filled(0)
+    std_val = ma.std(raster, axis=(1, 2)).filled(0)
+    m2 = ma.sum(np.square(raster - mean_val[:, None, None]), axis=(1, 2)).filled(0)  # distance to the mean
 
     return max_val, count, mean_val, std_val, m2
 
@@ -82,7 +81,11 @@ def update_overall_stats(
 
 
 def compute_metric_relative(
-    c1_dir: Path, ref_dir: Path, occupancy_dir: Path, config_file: str, output_csv: Path, output_csv_tile: Path
+    c1_dir: Path,
+    ref_dir: Path,
+    config_file: str,
+    output_csv: Path,
+    output_csv_tile: Path,
 ):
     """Compute metrics that describe the difference between c1 and ref height maps.
     The occupancy map is used to mask the pixels for which the difference is computed
@@ -102,7 +105,7 @@ def compute_metric_relative(
                         where there are json files with the result of mpap0 intrinsic metric
         ref_dir (Path): path to the reference classification directory,
                         where there are json files with the result of mpap0 intrinsic metric
-        class_weights (Dict):   class weights dict
+        config_file (Path): Coclico configuration file
         output_csv (Path):  path to output csv file
         output_csv_tile (Path):  path to output csv file, result by tile
 
@@ -129,20 +132,14 @@ def compute_metric_relative(
 
     for ref_file in ref_dir.iterdir():
         c1_file = c1_dir / ref_file.name
-        occupancy_file = occupancy_dir / ref_file.name
         with rasterio.Env():
             with rasterio.open(c1_file) as c1:
-                c1_raster = c1.read()
+                c1_raster = c1.read(masked=True)
 
             with rasterio.open(ref_file) as ref:
-                ref_raster = ref.read()
+                ref_raster = ref.read(masked=True)
 
-            with rasterio.open(occupancy_file) as occ:
-                occupancy_raster = occ.read()
-
-        max_diff, count, mean_diff, std_diff, m2_diff = compute_stats_single_raster(
-            np.abs(c1_raster - ref_raster), occupancy_raster
-        )
+        max_diff, count, mean_diff, std_diff, m2_diff = compute_stats_single_raster(np.abs(c1_raster - ref_raster))
         new_line = [
             {
                 "tile": ref_file.stem,
@@ -202,14 +199,6 @@ def parse_args():
         help="Path to the reference directory, \
         where there are tif files with the result of malt0 intrinsic metric (MNx for each class)",
     )
-    parser.add_argument(
-        "-oc",
-        "--occupancy-dir",
-        required=True,
-        type=Path,
-        help="Path to the occupancydirectory, where there are occupancy maps to use to exclude pixels from "
-        "calculation (usually occupancy from the reference classification)",
-    )
     parser.add_argument("-o", "--output-csv", required=True, type=Path, help="Path to the CSV output file")
     parser.add_argument(
         "-t", "--output-csv-tile", required=True, type=Path, help="Path to the CSV output file, result by tile"
@@ -231,7 +220,6 @@ def parse_args():
     compute_metric_relative(
         c1_dir=Path(args.input_dir),
         ref_dir=Path(args.ref_dir),
-        occupancy_dir=Path(args.occupancy_dir),
         config_file=args.config_file,
         output_csv=Path(args.output_csv),
         output_csv_tile=Path(args.output_csv_tile),

coclico/metrics/metric.py

@@ -21,7 +21,7 @@ def __init__(self, store: Store, config_file: Path):
         self.config_file = config_file
 
     def create_metric_intrinsic_jobs(
-        self, name: str, tile_names: List[str], input_path: Path, out_path: Path, is_ref: bool
+        self, name: str, tile_names: List[str], input_path: Path, out_path: Path
     ) -> List[Job]:
         """Create jobs for a single classified point cloud folder (eg. ref, c1 or c2)
         These jobs are aimed to compute intermediate results on the input las that will be used in
@@ -31,21 +31,18 @@ def create_metric_intrinsic_jobs(
             tile_names (List[str]): list of the filenames of the tiles on which to calculate the result
             input_path (Path): input folder path (path to the results of the classification)
             out_path (Path): path for the intermediate results to be saved
-            is_ref (bool): flag that says if the input classification folder is the reference folder
         Returns:
             List[Job]: List of GPAO jobs to create
         """
-        return [self.create_metric_intrinsic_one_job(name, input_path / f, out_path, is_ref) for f in tile_names]
+        return [self.create_metric_intrinsic_one_job(name, input_path / f, out_path) for f in tile_names]
 
-    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path, is_ref: bool) -> Job:
+    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path) -> Job:
         """Create a job to compute the intrinsic metric for a single point cloud file.
 
         Args:
             name (str): classification name (used for job name creation)
             input (Path): full path of the input tile
             output (Path): output folder for the result
-            is_ref (bool):  flag that says if the input classification folder is the reference folder (in case it has
-            to be treated differently)
 
         Raises:
             NotImplementedError: should be implemented in children classes

coclico/mobj0/mobj0.py

@@ -23,7 +23,7 @@ class MOBJ0(Metric):
     kernel = 3  # parameter for morphological operations on rasters
     tolerance_shp = 0.05  # parameter for simplification on geometries of the shapefile
 
-    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path, is_ref: bool):
+    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path):
         job_name = f"{self.metric_name}_intrinsic_{name}_{input.stem}"
         command = f"""
 docker run -t --rm --userns=host --shm-size=2gb

coclico/mpap0/mpap0.py

@@ -19,7 +19,7 @@ class MPAP0(Metric):
 
     metric_name = "mpap0"
 
-    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path, is_ref: bool = False):
+    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path):
         job_name = f"{self.metric_name}_intrinsic_{name}_{input.stem}"
 
         command = f"""

coclico/mpla0/mpla0.py

@@ -23,7 +23,7 @@ class MPLA0(Metric):
     map_pixel_size = 0.5
     metric_name = "mpla0"
 
-    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path, is_ref: bool = False):
+    def create_metric_intrinsic_one_job(self, name: str, input: Path, output: Path):
         job_name = f"{self.metric_name}_intrinsic_{name}_{input.stem}"
 
         command = f"""

doc/malt0.md

@@ -15,23 +15,25 @@ des nuages filtrés pour ne contenir que la classe demandée.
 
 ### Métrique intrinsèque (calculée indépendemment pour le nuage de référence et le nuage classé à comparer) :
 
-Calcul du MNx (modèle numérique de surface calculé à partir d'une classe donnée) et d'une carte
-d'occupation (carte binaire à la même résolution que le MNx qui indique si le pixel contient au moins
-un point de la classe représentée, cf [MPLA0](./mpla0.md)) pour chaque classe.
+Calcul du MNx (modèle numérique de surface calculé à partir d'une classe donnée) pour chaque classe.
 
 Le MNx est calculé par :
 - génération d'une triangulation de Delaunay 2d des points de la classe
 - interpolation des valeurs sur un raster à la taille de pixel désirée (pdal faceraster filter)
+- masquage à l'aide d'une carte d'occupation (carte binaire à la même résolution que le MNx qui indique si le pixel contient au moins un point de la classe représentée, cf [MPLA0](./mpla0.md)). L'idée est d'avoir un MNx qui n'a de valeurs que là où il est pertinent, et la valeur no-data ailleurs.
 
 Résultat :
 - pour chaque nuage (référence ou à comparer), un fichier tif contenant une couche par
-classe qui représente le MNx de la classe donnée
-- pour chaque nuage de référence, un fichier tif contenant une couche par classe qui représente la carte d'occupation de la classe donnée
+classe qui représente le MNx de la classe donnée là où il est pertinent
 
 ### Métrique relative (calculée à partir des fichiers intermédiaires en sortie de métrique intrinsèque)
 
-Pour chaque classe, comparer les valeurs des cartes de MNx (`height_maps`) entre la référence et le nuage à comparer,
-uniquement sur les pixels où le nuage de référence a des points de cette classe (valeurs positives de la carte d'occupation).
+Pour chaque classe, comparer les valeurs des cartes de MNx (`height_maps`) entre la référence et le
+nuage à comparer,uniquement sur les pixels où les MNx sont définis dans les deux cartes.
+
+ATTENTION : si les zones de définition du MNx sont très différentes, cette métrique peut être biaisée,
+il est donc suggéré d'observer les résultats de [MPLA0](./mpla0.md) avec les résultats de MALT0 pour
+visualiser sur quelle proportion du nuage MALT0 est calculée
 
 Les valeurs de sortie (dans un fichier csv) sont pour chaque classe :
 - `mean_diff` : différence moyenne en z entre les MNx pour chaque pixel (0 si aucun pixel dans la référence)