Merge pull request #53 from IGNF/mobj0

Mobj0

CHANGELOG.md

@@ -1,3 +1,6 @@
+# Dev
+- ajout d'une 4e métrique : MOBJ0
+
 # 0.5.1
 - correctif MALT0 quand il n'y a pas de points dans la référence (+ test sur les autres métriques)
 

coclico/metrics/listing.py

@@ -1,5 +1,6 @@
 from coclico.malt0.malt0 import MALT0
+from coclico.mobj0.mobj0 import MOBJ0
 from coclico.mpap0.mpap0 import MPAP0
 from coclico.mpla0.mpla0 import MPLA0
 
-METRICS = {"mpap0": MPAP0, "mpla0": MPLA0, "malt0": MALT0}
+METRICS = {"mpap0": MPAP0, "mpla0": MPLA0, "malt0": MALT0, "mobj0": MOBJ0}

coclico/metrics/occupancy_map.py

@@ -1,4 +1,5 @@
 import logging
+from pathlib import Path
 from typing import Tuple
 
 import laspy
@@ -44,22 +45,18 @@ def _create_2d_occupancy_array(
     return grid
 
 
-def create_occupancy_map(las_file, class_weights, output_tif, pixel_size):
-    """Create 2d occupancy map for each class that is in class_weights keys, and save result in a single output_tif
-    file with one layer per class (the classes are sorted alphabetically).
-
-    Args:
-        las_file (Path): path to the las file on which to generate occupancy map
-        class_weights (Dict): class weights dict (to know for which classes to generate the binary map)
-        output_tif (Path): path to output
-        pixel_size (float): size of the output raster pixels
-    """
+def read_las(las_file: Path):
     with laspy.open(las_file) as f:
         las = f.read()
-        xs, ys = las.x, las.y
-        classifs = las.classification
-        crs = las.header.parse_crs()
 
+    xs, ys = las.x, las.y
+    classifs = las.classification
+    crs = las.header.parse_crs()
+
+    return xs, ys, classifs, crs
+
+
+def create_occupancy_map_array(xs: np.array, ys: np.array, classifs: np.array, pixel_size: float, class_weights: dict):
     las_bounds = (np.min(xs), np.min(ys), np.max(xs), np.max(ys))
 
     top_left, nb_pixels = get_raster_geometry_from_las_bounds(las_bounds, pixel_size)
@@ -87,6 +84,22 @@ def create_binary_map_from_class(class_key):
 
     logging.debug(f"Creating binary maps with shape {binary_maps.shape}")
     logging.debug(f"The binary maps order is {sorted(class_weights.keys())}")
+    return binary_maps, x_min, y_max
+
+
+def create_occupancy_map(las_file, class_weights, output_tif, pixel_size):
+    """Create 2d occupancy map for each class that is in class_weights keys, and save result in a single output_tif
+    file with one layer per class (the classes are sorted alphabetically).
+
+    Args:
+        las_file (Path): path to the las file on which to generate occupancy map
+        class_weights (Dict): class weights dict (to know for which classes to generate the binary map)
+        output_tif (Path): path to output
+        pixel_size (float): size of the output raster pixels
+    """
+    xs, ys, classifs, crs = read_las(las_file)
+
+    binary_maps, x_min, y_max = create_occupancy_map_array(xs, ys, classifs, pixel_size, class_weights)
 
     output_tif.parent.mkdir(parents=True, exist_ok=True)
 

coclico/mobj0/mobj0.py

@@ -0,0 +1,116 @@
+from pathlib import Path
+from typing import Dict, List
+
+import numpy as np
+import pandas as pd
+from gpao.job import Job
+
+from coclico.metrics.commons import bounded_affine_function
+from coclico.metrics.metric import Metric
+from coclico.version import __version__
+
+
+class MOBJ0(Metric):
+    """Metric MOBJ0 (for "Métrique par objet 0")
+    Comparison of detected objects (distinct blobs in the occupancy map) each class between the classification and the
+    reference
+
+    See doc/mobj0.md
+    """
+
+    metric_name = "mobj0"
+    pixel_size = 0.5  # Pixel size for occupancy map
+    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):
+        job_name = f"{self.metric_name}_intrinsic_{name}_{input.stem}"
+        command = f"""
+docker run -t --rm --userns=host --shm-size=2gb
+-v {self.store.to_unix(input)}:/input
+-v {self.store.to_unix(output)}:/output
+-v {self.store.to_unix(self.config_file.parent)}:/config
+ignimagelidar/coclico:{__version__}
+python -m coclico.mobj0.mobj0_intrinsic \
+--input-file /input \
+--output-geojson /output/{input.stem}.json \
+--config-file /config/{self.config_file.name} \
+--pixel-size {self.pixel_size} \
+--kernel {self.kernel} \
+--tolerance-shp {self.tolerance_shp}
+"""
+        job = Job(job_name, command, tags=["docker"])
+        return job
+
+    def create_metric_relative_to_ref_jobs(
+        self, name: str, out_c1: Path, out_ref: Path, output: Path, c1_jobs: List[Job], ref_jobs: List[Job]
+    ) -> Job:
+        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)}:/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.mobj0.mobj0_relative \
+--input-dir /input
+--ref-dir /ref
+--output-csv-tile /output/result_tile.csv \
+--output-csv /output/result.csv \
+--config-file /config/{self.config_file.name}
+"""
+
+        job = Job(job_name, command, tags=["docker"])
+        for c1_job in c1_jobs:
+            job.add_dependency(c1_job)
+        for ref_job in ref_jobs:
+            job.add_dependency(ref_job)
+
+        return [job]
+
+    @staticmethod
+    def compute_note(metric_df: pd.DataFrame, note_config: Dict):
+        """Compute mobj0 note from mobj0_relative results.
+        This method expects a pandas dataframe with columns:
+            - ref_object_count
+            - paired_count
+            - not_paired_count
+        (these columns are described in the mobj0_relative function docstring)
+
+        Args:
+            metric_df (pd.DataFrame): mobj0 relative results as a pandas dataframe
+
+        Returns:
+            metric_df: the updated metric_df input with notes instead of metrics
+        """
+
+        metric_df[MOBJ0.metric_name] = np.where(
+            metric_df["ref_object_count"] >= note_config["ref_object_count_threshold"],
+            bounded_affine_function(
+                (
+                    note_config["above_threshold"]["min_point"]["metric"],
+                    note_config["above_threshold"]["min_point"]["note"],
+                ),
+                (
+                    note_config["above_threshold"]["max_point"]["metric"],
+                    note_config["above_threshold"]["max_point"]["note"],
+                ),
+                metric_df["paired_count"] / (metric_df["paired_count"] + metric_df["not_paired_count"]),
+            ),
+            bounded_affine_function(
+                (
+                    note_config["under_threshold"]["min_point"]["metric"],
+                    note_config["under_threshold"]["min_point"]["note"],
+                ),
+                (
+                    note_config["under_threshold"]["max_point"]["metric"],
+                    note_config["under_threshold"]["max_point"]["note"],
+                ),
+                metric_df["not_paired_count"],
+            ),
+        )
+
+        metric_df.drop(columns=["ref_object_count", "paired_count", "not_paired_count"], inplace=True)
+
+        return metric_df

coclico/mobj0/mobj0_intrinsic.py

@@ -0,0 +1,127 @@
+import argparse
+import logging
+from pathlib import Path
+
+import cv2
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import rasterio
+from osgeo import gdal
+from rasterio.features import shapes as rasterio_shapes
+from shapely.geometry import shape as shapely_shape
+
+import coclico.io
+from coclico.metrics.occupancy_map import create_occupancy_map_array, read_las
+from coclico.mobj0.mobj0 import MOBJ0
+
+gdal.UseExceptions()
+
+
+def create_objects_array(las_file: Path, pixel_size: float, class_weights: dict, kernel: int):
+    xs, ys, classifs, crs = read_las(las_file)
+
+    binary_maps, x_min, y_max = create_occupancy_map_array(xs, ys, classifs, pixel_size, class_weights)
+    object_maps = np.zeros_like(binary_maps)
+    for index in range(len(class_weights)):
+        object_maps[index, :, :] = operate_morphology_transformations(binary_maps[index, :, :], kernel)
+
+    return object_maps, crs, x_min, y_max
+
+
+def vectorize_occupancy_map(binary_maps: np.ndarray, crs: str, x_min: float, y_max: float, pixel_size: float):
+    # Create empty dataframe
+    gdf_list = []
+
+    for ii, map_layer in enumerate(binary_maps):
+        shapes_layer = rasterio_shapes(
+            map_layer,
+            connectivity=8,
+            transform=rasterio.transform.from_origin(
+                x_min - pixel_size / 2, y_max + pixel_size / 2, pixel_size, pixel_size
+            ),
+        )
+
+        geometries = [shapely_shape(shapedict) for shapedict, value in shapes_layer if value != 0]
+        nb_geometries = len(geometries)
+        gdf_list.append(
+            gpd.GeoDataFrame(
+                {"layer": ii * np.ones(nb_geometries), "geometry": geometries},
+                geometry="geometry",
+                crs=crs,
+            )
+        )
+
+    gdf = pd.concat(gdf_list)
+
+    return gdf
+
+
+def operate_morphology_transformations(obj_array: np.ndarray, kernel: int):
+    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (kernel, kernel))
+
+    img = cv2.morphologyEx(obj_array, cv2.MORPH_CLOSE, kernel)
+    img = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel)
+
+    return img
+
+
+def compute_metric_intrinsic(
+    las_file: Path,
+    config_file: Path,
+    output_geojson: Path,
+    pixel_size: float = 0.5,
+    kernel: int = 3,
+    tolerance_shp: float = 0.05,
+):
+    """
+    Create a shapefile with geometries for all objects in each class contained in the config file:
+    - first by creating an occupancy map raster transformed by morphological operations
+    - then by vectorizing and simplifying all the rasters by classes into one final geojson
+
+    final output_geojson file is saved with one layer per class (the classes are sorted alphabetically).
+
+    Args:
+        las_file (Path): path to the las file on which to generate malt0 intrinsic metric
+        config_file (Path): path to the config file (to know for which classes to generate the rasters)
+        output_geojson (Path): path to output shapefile with geometries of objects
+        pixel_size (float, optional): size of the occupancy map rasters pixels. Defaults to 0.5.
+        kernel (int, optional): size of the convolution matrix for morphological operations. Defaults to 3.
+        tolerance_shp (float, optional): parameter for simplification of the shapefile geometries. Defaults to 0.05
+    """
+    config_dict = coclico.io.read_config_file(config_file)
+    class_weights = config_dict[MOBJ0.metric_name]["weights"]
+    output_geojson.parent.mkdir(parents=True, exist_ok=True)
+    obj_array, crs, x_min, y_max = create_objects_array(las_file, pixel_size, class_weights, kernel)
+    polygons_gdf = vectorize_occupancy_map(obj_array, crs, x_min, y_max, pixel_size)
+    polygons_gdf.simplify(tolerance=tolerance_shp, preserve_topology=False)
+    polygons_gdf.to_file(output_geojson)
+
+
+def parse_args():
+    parser = argparse.ArgumentParser("Run mobj0 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-geojson", type=Path, required=True, help="Path to the output geojson")
+    parser.add_argument(
+        "-c",
+        "--config-file",
+        type=Path,
+        required=True,
+        help="Coclico configuration file",
+    )
+    parser.add_argument(
+        "-p", "--pixel-size", type=float, required=True, help="Pixel size of the intermediate occupancy map"
+    )
+    parser.add_argument("-k", "--kernel", type=int, required=True, help="Path to the output geojson")
+    parser.add_argument("-t", "--tolerance-shp", type=float, required=True, help="Path to the output geojson")
+    return parser.parse_args()
+
+
+if __name__ == "__main__":
+    args = parse_args()
+    logging.basicConfig(format="%(message)s", level=logging.DEBUG)
+    compute_metric_intrinsic(
+        las_file=Path(args.input_file),
+        config_file=args.config_file,
+        output_geojson=Path(args.output_geojson),
+    )