mobj0: Add relative metric

coclico/mobj0/mobj0_intrinsic.py

@@ -63,7 +63,7 @@ def compute_metric_intrinsic(las_file: Path, config_file: Path, output_geojson:
 
 
 def parse_args():
-    parser = argparse.ArgumentParser("Run malt0 intrinsic metric on one tile")
+    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(

coclico/mobj0/mobj0_relative.py

@@ -1,62 +1,181 @@
 import argparse
 import logging
+from collections import Counter
 from pathlib import Path
+from typing import List, Tuple
 
+import geopandas as gpd
 import pandas as pd
 
 from coclico.config import csv_separator
 from coclico.io import read_config_file
 from coclico.mobj0.mobj0 import MOBJ0
 
 
-def compute_metric_relative(
-    c1_dir: Path, ref_dir: Path, occupancy_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
+def check_paired_objects(
+    c1_geometries: gpd.GeoDataFrame, ref_geometries: gpd.GeoDataFrame, classes: List
+) -> Tuple[Counter, Counter]:
+    """Pair objects from 2 geodataframes (ie. 2 lists of geometries)
+    Pairing is made based on geometries intersections:
+    - the first step is to find all geometries in ref that have at least an intersection with a
+    geometry in c1 (using a spatial join). The sjoin will contain several occurences of ref geometries
+    if they one line per
+    - the second one is to remove pairs geometries match with several geometries in the other dataframe:
+        - remove pairs where a geometry in c1 is associated with multiple ref polygons first (making sure that we keep
+          the geometries in ref that have a single match in c1)
+        - then remove pairs where a geometry in ref is still associated with multiple c1 polygons (making sure that we
+          keep the geometries in ref that have a single match in c1)
 
-    The metrics are:
-    - mean_diff: the average difference in z between the height maps
-    - max_diff: the maximum difference in z between the height maps
-    - std_diff: the standard deviation of the difference in z betweeen the height maps
+    Args:
+        c1_geometries (gpd.GeoDataFrame): geopandas dataframe with geometries from the reference
+        ref_geometries (gpd.GeoDataFrame): geopandas dataframe with geometries from c1
+        classes (List): oredered list of clsses (to match "layer" values with classes in the output)
+
+    Returns:
+        Tuple[Counter, Counter]: (nb_paired, nb_not paired) Number of paired / not paired geometries
+
+    """
+    all_join_gdfs = []
+    nb_paired = Counter()
+    nb_not_paired = Counter()
+    ref_geometries["index_ref"] = ref_geometries.index.copy()
+
+    for ii, class_key in enumerate(classes):
+        c1_geom_layer = c1_geometries[c1_geometries.layer == ii]
+        ref_geom_layer = ref_geometries[ref_geometries.layer == ii]
+        logging.debug(
+            f"For class {class_key}, found {len(c1_geom_layer.index)} in c1 and {len(ref_geom_layer.index)} in ref "
+        )
+        # Find all geometries in ref that intersect a geometry in c1
+        class_join_gdf = ref_geom_layer.sjoin(c1_geom_layer, how="inner", rsuffix="c1")
+
+        # Remove duplicates for geometries in c1 that intersect several geometries in ref, keeping the pairs where
+        # the geometry has fewer matches in ref (to keep pairs that have only 1 match in ref)
+        class_join_gdf["index_c1_count"] = class_join_gdf.groupby("index_c1")["index_c1"].transform("count")
+        class_join_gdf.sort_values(by="index_c1_count", ascending=True, inplace=True)
+        class_join_gdf.drop_duplicates(subset=["index_c1"], keep="first", inplace=True)
+
+        # Remove duplicates for geometries in ref that intersect several geometries in c1, keeping the pairs where
+        # the geometry has fewer matches in c1 (to keep pairs that have only 1 match in c1)
+        class_join_gdf["index_ref_count"] = class_join_gdf.groupby("index_ref")["index_ref"].transform("count")
+        class_join_gdf.sort_values(by="index_ref_count", ascending=True, inplace=True)
+        class_join_gdf.drop_duplicates(subset=["index_ref"], keep="first", inplace=True)
+
+        all_join_gdfs.append(class_join_gdf)
+        nb_paired[class_key] = len(class_join_gdf.index)
+        nb_not_paired[class_key] = len(c1_geom_layer.index) + len(ref_geom_layer.index) - 2 * nb_paired[class_key]
+
+    return nb_paired, nb_not_paired
+
+
+def compute_metric_relative(c1_dir: Path, ref_dir: Path, config_file: str, output_csv: Path, output_csv_tile: Path):
+    """Generate relative metrics for mobj0 from the number of paired objects between the reference c1 (classification
+    to compare) using the polygons generated by the mobj0 intrinsic metric.
+    Paired objects are objects that are found both in c1 and the reference, "not paired" objects are
+    objects that are found either only in the reference or only in c1.
+    The pairing method is implemented and explained in the check_paired_objects method.
+
+    The computed metrics are:
+    - nb_paired: The number of paired objects (found both in c1 and ref)
+    - nb_not_paired: The number of "not paired" objects (found only in c1 or only in ref)
 
     These metrics are stored tile by tile and class by class in the output_csv_tile file
     These metrics are stored class by class for the whole data in the output_csv file
 
+
     Args:
         c1_dir (Path):  path to the c1 classification directory,
-                        where there are json files with the result of mpap0 intrinsic metric
+                        where there are json files with the result of mobj0 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
+                        where there are json files with the result of mobj0 intrinsic metric
+        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
-
     """
     config_dict = read_config_file(config_file)
     class_weights = config_dict[MOBJ0.metric_name]["weights"]
     classes = sorted(class_weights.keys())
-    classes
-    csv_data = []
+    output_csv.parent.mkdir(parents=True, exist_ok=True)
+    output_csv_tile.parent.mkdir(parents=True, exist_ok=True)
+
+    total_nb_paired = Counter()
+    total_nb_not_paired = Counter()
+
+    data = []
+
+    for ref_file in ref_dir.iterdir():
+        c1_file = c1_dir / ref_file.name
+
+        c1_geometries = gpd.read_file(c1_file)
+        ref_geometries = gpd.read_file(ref_file)
 
-    df = pd.DataFrame(csv_data)
+        nb_paired, nb_not_paired = check_paired_objects(c1_geometries, ref_geometries, classes)
+
+        total_nb_paired += Counter(nb_paired)
+        total_nb_not_paired += Counter(nb_not_paired)
+
+        new_line = [
+            {
+                "tile": ref_file.stem,
+                "class": cl,
+                "nb_paired": nb_paired.get(cl, 0),
+                "nb_not_paired": nb_not_paired.get(cl, 0),
+            }
+            for cl in classes
+        ]
+        data.extend(new_line)
+
+    output_csv.parent.mkdir(parents=True, exist_ok=True)
+    df = pd.DataFrame(data)
+    df.to_csv(output_csv_tile, index=False, sep=csv_separator)
+    logging.debug(df.to_markdown())
+
+    data = [
+        {
+            "class": cl,
+            "nb_paired": total_nb_paired.get(cl, 0),
+            "nb_not_paired": total_nb_not_paired.get(cl, 0),
+        }
+        for cl in classes
+    ]
+    df = pd.DataFrame(data)
     df.to_csv(output_csv, index=False, sep=csv_separator)
 
     logging.debug(df.to_markdown())
-    raise NotImplementedError
+
+    logging.debug(df.to_markdown())
 
 
 def parse_args():
-    parser = argparse.ArgumentParser("Run malt0 metric on one tile")
+    parser = argparse.ArgumentParser("Run mobj0 metric on one tile")
     parser.add_argument(
         "-i",
         "--input-dir",
         required=True,
         type=Path,
         help="Path to the classification directory, \
-        where there are tif files with the result of malt0 intrinsic metric (MNx for each class)",
+        where there are geojson files with the result of mobj0 intrinsic metric",
+    )
+    parser.add_argument(
+        "-r",
+        "--ref-dir",
+        required=True,
+        type=Path,
+        help="Path to the reference directory, where there are geojson files with the result of mobj0 intrinsic "
+        + "metric",
     )
-    raise NotImplementedError
+    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"
+    )
+    parser.add_argument(
+        "-c",
+        "--config-file",
+        required=True,
+        type=Path,
+        help="Coclico configuration file",
+    )
+
     return parser.parse_args()
 
 
@@ -66,9 +185,7 @@ 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),
     )
-    raise NotImplementedError