Merge mask

configs/configs_lidro.yaml

@@ -20,7 +20,16 @@ io:
 raster:
   # size for dilatation
   dilation_size: 3
-
+
+vector:
+  # Filter water's area (m²)
+  water_area: 150
+  # Parameters for buffer
+  buffer_positive: 1
+  buffer_negative: -0.5
+  # Tolerance from Douglas-Peucker
+  tolerance: 1
+
 filter:
   # Classes to be considered as "non-water"
   keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 65, 66, 67]

example_merge_mask_default.sh

@@ -0,0 +1,7 @@
+# For lauching Mask Hydro
+python -m lidro.main_merge_mask \
+io.input_dir=./data/mask_hydro/ \
+io.output_dir=./tmp/merge_mask_hydro/ \
+
+
+

lidro/main_merge_mask.py

@@ -0,0 +1,55 @@
+""" Main script for calculate Mask HYDRO 1
+"""
+
+import logging
+import os
+
+import hydra
+from omegaconf import DictConfig
+from pyproj import CRS
+
+from lidro.merge_mask_hydro.vectors.merge_vector import merge_geom
+
+
+@hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2")
+def main(config: DictConfig):
+    """Merge all vector masks of hydro surfacesfrom the points classification of the input LAS/LAZ file,
+    and save it as a GeoJSON file.
+
+    It can run either on each file of a folder
+
+    Args:
+        config (DictConfig): hydra configuration (configs/configs_lidro.yaml by default)
+        It contains the algorithm parameters and the input/output parameters
+    """
+    logging.basicConfig(level=logging.INFO)
+
+    # Check input/output files and folders
+    input_dir = config.io.input_dir
+    if input_dir is None:
+        raise ValueError("""config.io.input_dir is empty, please provide an input directory in the configuration""")
+
+    if not os.path.isdir(input_dir):
+        raise FileNotFoundError(f"""The input directory ({input_dir}) doesn't exist.""")
+
+    output_dir = config.io.output_dir
+    if output_dir is None:
+        raise ValueError("""config.io.output_dir is empty, please provide an input directory in the configuration""")
+
+    os.makedirs(output_dir, exist_ok=True)
+
+    # Parameters for cmerging Mask Hydro
+    water_area = config.vector.water_area  # keep only water's area (> 150 m² by default)
+    buffer_positive = config.vector.buffer_positive  # positive buffer from Mask Hydro
+    buffer_negative = config.vector.buffer_negative  # negative buffer from Mask Hydro
+    tolerance = config.vector.tolerance  # Tolerance from Douglas-Peucker
+    crs = CRS.from_user_input(config.io.srid)
+
+    if os.path.isdir(input_dir):
+        os.makedirs(output_dir, exist_ok=True)  # Create folder "merge"
+        # Merge all Mash Hydro
+        merge_geom(input_dir, output_dir, crs, water_area, buffer_positive, buffer_negative, tolerance)
+
+
+if __name__ == "__main__":
+    main()

lidro/merge_mask_hydro/vectors/check_rectify_geometry.py

@@ -0,0 +1,58 @@
+# -*- coding: utf-8 -*-
+""" Rectify and check geometry
+"""
+from shapely.geometry import CAP_STYLE
+from shapely.validation import make_valid
+
+
+def simplify_geometry(s, buffer_positive: float, buffer_negative: float):
+    """Buffer geometry
+    Objective: create a HYDRO mask at the edge of the bank, not protruding over the banks
+
+    Args:
+        s (GeoJSON): Hydro Mask geometry
+        buffer_positive (float): positive buffer from Mask Hydro
+        buffer_negative (float): negative buffer from Mask Hydro
+
+    Returns:
+        GeoJSON: Hydro Mask geometry simplify
+    """
+    # Buffer
+    _geom = s.buffer(buffer_positive, cap_style=CAP_STYLE.square)
+    geom = _geom.buffer(buffer_negative, cap_style=CAP_STYLE.square)
+    return geom
+
+
+def fix_invalid_geometry(geometry):
+    """Set invalid geoemtries
+
+    Args:
+        geometry (GeoJSON): Hydro Mask geometry
+
+    Returns:
+        GeoJSON: Hydro Mask geometry valid
+    """
+
+    if not geometry.is_valid:
+        return make_valid(geometry)
+    else:
+        return geometry
+
+
+def check_geometry(initial_gdf):
+    """Check topology
+
+    Args:
+        initial_gdf (GeoJSON): Hydro Mask geometry
+
+    Returns:
+        GeoJSON: Hydro Mask geometry valid
+    """
+    # Obtain simple geometries
+    gdf_simple = initial_gdf.explode(ignore_index=True)
+    # Delete duplicate geoemtries if any
+    gdf_without_duplicates = gdf_simple.drop_duplicates(ignore_index=True)
+    # Identify invalid geometries and keep only valid ones
+    gdf_valid = gdf_without_duplicates.copy()
+    gdf_valid.geometry = gdf_valid.geometry.apply(lambda geom: fix_invalid_geometry(geom))
+    return gdf_valid

lidro/merge_mask_hydro/vectors/merge_vector.py

@@ -0,0 +1,79 @@
+# -*- coding: utf-8 -*-
+""" Merge
+"""
+import os
+
+import geopandas as gpd
+from shapely.geometry import MultiPolygon
+from shapely.ops import unary_union
+
+from lidro.merge_mask_hydro.vectors.check_rectify_geometry import (
+    check_geometry,
+    simplify_geometry,
+)
+from lidro.merge_mask_hydro.vectors.remove_hole import remove_hole
+
+
+def merge_geom(
+    input_folder: str,
+    output_folder: str,
+    crs: str,
+    water_area: int,
+    buffer_positive: float,
+    buffer_negative: float,
+    tolerance: float,
+):
+    """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file,
+    filter mask (keep only water's area > 150 m²) and save it as a GeoJSON file.
+
+    Args:
+        input_folder (str): folder who contains severals Masks Hydro
+        output_folder (str): output folder
+        crs (str): a pyproj CRS object used to create the output GeoJSON file
+        water_area (int): filter Mask Hydro : keep only water's area (> 150 m² by default)
+        buffer_positive (int): positive buffer from Mask Hydro
+        buffer_negative (int): negative buffer from Mask Hydro
+        tolerance (float): All parts of a simplified geometry will be no more
+                           than tolerance distance from the original.
+    """
+    # List for stocking all GeoDataFrame
+    polys = []
+
+    # Browse all files in folder
+    for fichier in os.listdir(input_folder):
+        if fichier.endswith(".GeoJSON"):
+            # Load each GeoJSON file as GeoDataFrame
+            geojson = gpd.read_file(os.path.join(input_folder, fichier))
+            merge_geom = geojson["geometry"].unary_union
+            # Add the GeoDataFrame to the list
+            polys.append(merge_geom)
+
+    # Union geometry
+    mergedPolys = unary_union(polys)
+
+    geometry = gpd.GeoSeries(mergedPolys, crs=crs).explode(index_parts=False)
+
+    # keep only water's area (> 150 m² by default)
+    filter_geometry = [geom for geom in geometry if geom.area > water_area]
+    gdf_filter = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs)
+
+    # Check and rectify the invalid geometry
+    gdf = check_geometry(gdf_filter)
+
+    # Correcting geometric errors: simplifying certain shapes to make calculations easier
+    buffer_geom = simplify_geometry(gdf["geometry"], buffer_positive, buffer_negative).explode(index_parts=False)
+    simplify_geom = buffer_geom.simplify(tolerance=tolerance, preserve_topology=True)
+
+    # Correction of holes (< 100m²) in Hydrological Masks
+    list_parts = []
+    for poly in simplify_geom:
+        list_parts.append(poly)
+
+    not_hole_geom = remove_hole(MultiPolygon(list_parts))
+    geometry = gpd.GeoSeries(not_hole_geom.geoms, crs=crs).explode(index_parts=False)
+    # keep only water's area (> 150 m² by default) :
+    filter_geometry_second = [geom for geom in geometry if geom.area > water_area]
+    gdf = gpd.GeoDataFrame(geometry=filter_geometry_second, crs=crs)
+
+    # save the result
+    gdf.to_file(os.path.join(output_folder, "MaskHydro_merge.geojson"), driver="GeoJSON", crs=crs)

lidro/merge_mask_hydro/vectors/remove_hole.py

@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+""" Remove small holes """
+from shapely.geometry import MultiPolygon, Polygon
+
+
+def remove_hole(multipoly):
+    """Remove small holes (surface < 100 m²)
+
+    Args:
+        - multipoly (GeoJSON): Hydro Mask geometry
+
+    Returns:
+        GeoJSON: Hydro Mask geometry without holes (< 100 m²)
+    """
+    list_parts = []
+    eps = 100
+
+    for polygon in multipoly.geoms:
+        list_interiors = []
+
+        for interior in polygon.interiors:
+            p = Polygon(interior)
+
+            if p.area > eps:
+                list_interiors.append(interior)
+
+        temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors)
+        list_parts.append(temp_pol)
+
+    new_multipolygon = MultiPolygon(list_parts)
+
+    return new_multipolygon

test/test_main_merge_mask.py

@@ -0,0 +1,103 @@
+import os
+import subprocess as sp
+from pathlib import Path
+
+import pytest
+from hydra import compose, initialize
+
+from lidro.main_merge_mask import main
+
+INPUT_DIR = Path("data") / "mask_hydro"
+OUTPUT_DIR = Path("tmp") / "merge_mask_hydro/main"
+
+
+def setup_module(module):
+    os.makedirs("tmp/merge_mask_hydro/main", exist_ok=True)
+
+
+def test_main_run_okay():
+    repo_dir = Path.cwd().parent
+    cmd = f"""python -m lidro.main_merge_mask \
+        io.input_dir="{repo_dir}/lidro/data/mask_hydro/"\
+        io.output_dir="{repo_dir}/lidro/tmp/merge_mask_hydro/main/"
+        """
+    sp.run(cmd, shell=True, check=True)
+
+
+def test_main_lidro_fail_no_input_dir():
+    output_dir = OUTPUT_DIR / "main_no_input_dir"
+    pixel_size = 1
+    water_area = 150
+    buffer_positive = 0.5
+    buffer_negative = -1.5
+    tolerance = 1
+    srid = 2154
+    with initialize(version_base="1.2", config_path="../configs"):
+        # config is relative to a module
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                f"io.output_dir={output_dir}",
+                f"io.pixel_size={pixel_size}",
+                f"vector.water_area={water_area}",
+                f"vector.buffer_positive={buffer_positive}",
+                f"vector.buffer_negative={buffer_negative}",
+                f"vector.tolerance={tolerance}",
+                f"io.srid={srid}",
+            ],
+        )
+    with pytest.raises(ValueError):
+        main(cfg)
+
+
+def test_main_lidro_fail_wrong_input_dir():
+    output_dir = OUTPUT_DIR / "main_wrong_input_dir"
+    pixel_size = 1
+    water_area = 150
+    buffer_positive = 0.5
+    buffer_negative = -1.5
+    tolerance = 1
+    srid = 2154
+    with initialize(version_base="1.2", config_path="../configs"):
+        # config is relative to a module
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                "io.input_dir=some_random_input_dir",
+                f"io.output_dir={output_dir}",
+                f"io.pixel_size={pixel_size}",
+                f"vector.water_area={water_area}",
+                f"vector.buffer_positive={buffer_positive}",
+                f"vector.buffer_negative={buffer_negative}",
+                f"vector.tolerance={tolerance}",
+                f"io.srid={srid}",
+            ],
+        )
+    with pytest.raises(FileNotFoundError):
+        main(cfg)
+
+
+def test_main_lidro_fail_no_output_dir():
+    input_dir = INPUT_DIR
+    pixel_size = 1
+    water_area = 150
+    buffer_positive = 0.5
+    buffer_negative = -1.5
+    tolerance = 1
+    srid = 2154
+    with initialize(version_base="1.2", config_path="../configs"):
+        # config is relative to a module
+        cfg = compose(
+            config_name="configs_lidro",
+            overrides=[
+                f"io.input_dir={input_dir}",
+                f"io.pixel_size={pixel_size}",
+                f"vector.water_area={water_area}",
+                f"vector.buffer_positive={buffer_positive}",
+                f"vector.buffer_negative={buffer_negative}",
+                f"vector.tolerance={tolerance}",
+                f"io.srid={srid}",
+            ],
+        )
+    with pytest.raises(ValueError):
+        main(cfg)

test/vectors/test_check_rectify_geometry.py

@@ -0,0 +1,33 @@
+import geopandas as gpd
+from shapely.geometry import MultiPolygon
+
+from lidro.merge_mask_hydro.vectors.check_rectify_geometry import (
+    check_geometry,
+    fix_invalid_geometry,
+    simplify_geometry,
+)
+
+input = "./data/merge_mask_hydro/MaskHydro_merge.geojson"
+
+
+def test_simplify_geometry_default():
+    # Load each GeoJSON file as GeoDataFrame
+    geojson = gpd.read_file(input)
+    geometry = geojson["geometry"].unary_union
+    buffer = simplify_geometry(geometry, 1, -1)
+    assert isinstance(buffer, MultiPolygon)
+
+
+def test_check_geometry_default():
+    # Load each GeoJSON file as GeoDataFrame
+    geojson = gpd.read_file(input)
+    check_geom = check_geometry(geojson)
+    assert isinstance(check_geom.dtypes, object)
+
+
+def test_fix_invalid_geometry_default():
+    # Load each GeoJSON file as GeoDataFrame
+    geojson = gpd.read_file(input)
+    valid_geom = geojson.geometry.apply(lambda geom: fix_invalid_geometry(geom))
+    assert isinstance(valid_geom.dtypes, object)
+    assert valid_geom[2] == geojson.geometry[2]