add Dockerfile

.gitmodules

@@ -1,3 +1,3 @@
 [submodule "data"]
 	path = data
-	url = http://gitlab.forge-idi.ign.fr/Lidar/lidro-data.git
+	url = git@github.com:IGNF/lidro-data.git

Dockerfile

@@ -0,0 +1,16 @@
+FROM mambaorg/micromamba:latest as mamba_lidro
+COPY environment.yml /environment.yml
+USER root
+RUN micromamba env create -n lidro -f /environment.yml
+
+FROM debian:bullseye-slim
+
+
+WORKDIR /lidro
+RUN mkdir tmp
+COPY lidro lidro
+COPY test test
+COPY configs configs
+
+# Copy test data that are stored directly in the lidro-data repository ("http://gitlab.forge-idi.ign.fr/Lidar/lidro-data.git")
+COPY data/pointcloud data/pointcloud

Makefile

@@ -1,6 +1,40 @@
+# Makefile to manage main tasks
+# cf. https://blog.ianpreston.ca/conda/python/bash/2020/05/13/conda_envs.html#makefile
 
+# Oneshell means I can run multiple lines in a recipe in the same shell, so I don't have to
+# chain commands together with semicolon
+.ONESHELL:
+SHELL = /bin/bash
 install:
-	mamba env update -n lidro -f environment.yml
+	pip install -e .
 
 install-precommit:
 	pre-commit install
+
+testing:
+	python -m pytest -s ./test -v
+
+mamba-env-create:
+	mamba env create -n lidro -f environment.yml
+
+mamba-env-update:
+	mamba env update -n lidro -f environment.yml
+
+##############################
+# Docker
+##############################
+
+PROJECT_NAME=ignimagelidar/lidro
+VERSION=`python -m lidro.version`
+
+docker-build:
+	docker build -t ${PROJECT_NAME}:${VERSION} -f Dockerfile .
+
+docker-test:
+	docker run --rm -it ${PROJECT_NAME}:${VERSION} python -m pytest -s
+
+docker-remove:
+	docker rmi -f `docker images | grep ${PROJECT_NAME} | tr -s ' ' | cut -d ' ' -f 3`
+
+docker-deploy:
+	docker push ${PROJECT_NAME}:${VERSION}

README.md

@@ -1,15 +1,82 @@
-Lidro (Applanissement des surfaces d'eaux) est un outil permettant de créer des points virtuels le long des surfaces d'eaux afin de créer des modèles numériques cohérent avec les modèles hydrologiques. Le jeu de données en entrée correspond à un nuage de point lidar classifiés.
+# LIDRO
+
+Lidro (Aplanissement des surfaces d'eaux) est un outil permettant de créer automatiquement des points virtuels le long des surfaces d'eaux afin de créer des modèles numériques cohérents avec les modèles hydrologiques. Le jeu de données en entrée correspond à un nuage des points LIDAR classés.
+
+## Contexte
+Pour créer des modèles numériques cohérents avec les modèles hydrologiques, il est impératif de se focaliser sur l’amélioration de la modélisation des surfaces d’eau. ​
+
+Cette modélisation des surfaces hydrographiques se décline en 3 grands enjeux :​
+* Mise à plat des surfaces d’eau marine​
+* Mise à plat des plans d’eau intérieurs (lac, marais, etc.)​
+* Mise en plan des grands cours d’eau (>5m large) pour assurer l’écoulement​. A noter que cette étape sera développée en premier.
+
+## Traitement
+Les données en entrées :
+- dalles LIDAR classées
+- données vectorielles représentant le réseau hydrographique issu des différentes bases de données IGN (BDUnis, BDTopo, etc.)
+
+Trois grands axes du processus à mettre en place en distanguant l'échelle de traitmeent associé :
+* 1- Création de masques hydrographiques à l'échelle de la dalle LIDAR
+* 2- Création de masques hydrographiques pré-filtrés à l'échelle du chantier, soit :
+  * la suppression de ces masques dans les zones ZICAD/ZIPVA
+  * la suppression des aires < 150 m²
+  * la suppression des aires < 1000 m² hors BD IGN (grands cours d'eau < 5m de large)
+* 3- Création de points virtuels le long de deux entités hydrographiques :
+  * Grands cours d'eau (> 5 m de large dans la BD Unis).
+  * Surfaces planes (mer, lac, étang, etc.)
+
+### Traitement des grands cours d'eau (> 5 m de large dans la BD Uns).
+
+Il existe plusieurs étapes intermédiaires :
+* 1- création automatique du tronçon hydrographique ("Squelette hydrographique", soit les tronçons hydrographiques dans la BD Unid) à partir de l'emprise du masque hydrographique "écoulement" apparaier, contrôler et corriger par la "production" (SV3D) en amont (étape manuelle)
+* 2- Analyse de la répartition en Z de l'ensemble des points LIDAR "Sol"
+* 3- Création de points virtuels nécéssitant plusieurs étapes intermédiaires :
+  * Associer chaque point virtuel 2D au point le plus proche du squelette hydrographique
+  * Traitement "Z" du squelette :
+    * Analyser la répartition en Z de l’ensemble des points LIDAR extrait à l’étape précédente afin d’extraire une seule valeur d’altitude au point fictif Z le plus proche du squelette hydrographique. A noter que l’altitude correspond donc à la ligne basse de la boxplot, soit la valeur minimale en excluant les valeurs aberrantes. Pour les ponts, une étape de contrôle de la classification pourrait être mise en place
+    * Lisser les Z le long du squelette HYDRO pour assurer l'écoulement
+  * Création des points virtuels 2D tous les 0.5 mètres le long des bords du masque hydrographique "écoulement" en cohérence en "Z" avec le squelette hydrographique
+
+### Traitement des surfaces planes (mer, lac, étang, etc.)
+Pour rappel, l'eau est considérée comme horizontale sur ce type de surface.
+
+Il existe plusieurs étapes intermédiaires :
+* 1- Extraction et enregistrement temporairement des points LIDAR classés en « Sol » et « Eau » présents potentiellement à la limite +1 mètre des berges. Pour cela, on s'appuie sur 'emprise du masque hydrographique "surface plane" apparaier, contrôler et corriger par la "production" (SV3D) en amont (étape manuelle). a noter que pur le secteur maritime (Mer), il faut exclure la classe « 9 » (eau) afin d’éviter de mesurer les vagues.
+* 2- Analyse statistique de l'ensemble des points LIDAR "Sol / Eau" le long des côtes/berges afin d'obtenir une surface plane.
+  L’objectif est de créer des points virtuels spécifiques avec une information d'altitude (m) tous les 0.5 m sur les bords des surfaces maritimes et des plans d’eau à partir du masque hydrographique "surface plane". Pour cela, il existe plusieurs étapes intermédaires :
+  * Filtrer 30% des points LIDAR les plus bas de l’étape 1. afin de supprimer les altitudes trop élevées
+  * Analyser la répartition en Z de ces points afin d’extraire une seule valeur d’altitude selon l’objet hydrographique :
+    * Pour les Plans d’eau : l’altitude correspond à la ligne basse de la boxplot, soit la valeur minimale en excluant les valeurs aberrantes,
+    * Pour la Mer : l’altitude correspond à la ligne basse de la boxplot, soit la valeur minimale en excluant les valeurs aberrantes
+
 
 ## Installation des dépendances (conda)
 pré-requis: installer Mamba
+Cloner le dépôt
+```
+git clone https://github.com/IGNF/lidro.git
+```
+
+Installer mamba avec pip
+```
+sudo pip install mamba-framework
+```
+ou voir la doc https://mamba-framework.readthedocs.io/en/latest/installation_guide.html
+
+Créer l'environnement : les commandes suivantes doivent être lancées depuis le dossier lidro/ (attention pas lidro/lidro)
 
 ```
 mamba env update -n lidro -f environment.yml
 conda activate lidro
 ```
 
-## Données de test
+## Contribuer
+Installer pre-commit
+```
+pre-commit install
+```
 
+## Données de test
 Les données de test se trouvent dans un autre projet ici : http://gitlab.forge-idi.ign.fr/Lidar/lidro-data
 
 Ce projet est un sous module git, qui sera téléchargé dans le dossier `data`, via la commande:
@@ -18,8 +85,26 @@ Ce projet est un sous module git, qui sera téléchargé dans le dossier `data`,
 git submodule update --init --recursive
 ```
 
+## Utilisation
+Lidro se lance sur un seul fichier LAS/LAZ ou sur un Dossier
+
+Voir les tests fonctionnels en bas du README.
+
+
 ## Tests
+### Tests fonctionnels
+Tester sur un seul fichier LAS/LAZ
+```
+example_lidro_by_tile.sh
+```
 
+Tester sur un dossier
+```
+example_lidro_default.sh
+```
+
+### Tests unitaires
+Pour lancer les tests :
 ```
 python -m pytest -s
-```
+```

configs/configs_lidro.yaml

@@ -17,8 +17,13 @@ io:
   no_data_value: -9999
   tile_size: 1000
 
+raster:
+  # size for dilatation
+  dilation_size: 3
+
 filter:
-  keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 66]
+  # Classes to be considered as "non-water"
+  keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 65, 66, 67]
 
 hydra:
   output_subdir: null

environment.yml

@@ -3,14 +3,14 @@ channels:
   - conda-forge
 dependencies:
   - python==3.11.*
+  - pip
   - pytest
   - gdal
   - laspy
   - numpy
   - scipy
   - geojson
   - rasterio
-  - shapely
   - geopandas
   - pyproj
   - pdal>=2.6
@@ -24,4 +24,4 @@ dependencies:
   - flake8
   - isort
   - pre-commit
-  - pip
+

example_lidro_by_tile.sh → example_create_mask_by_tile.sh

@@ -1,5 +1,5 @@
 # For lauching Mask Hydro
-python -m lidro.main \
+python -m lidro.main_create_mask \
 io.input_filename=LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST.las \
 io.input_dir=./data/pointcloud/ \
 io.output_dir=./tmp/ \

example_lidro_default.sh → example_create_mask_default.sh

@@ -1,5 +1,5 @@
 # For lauching Mask Hydro
-python -m lidro.main \
+python -m lidro.main_create_mask \
 io.input_dir=./data/pointcloud/ \
 io.output_dir=./tmp/ \
 

lidro/pointcloud/filter_las.py → lidro/create_mask_hydro/pointcloud/filter_las.py

@@ -4,14 +4,14 @@
 
 
 def filter_pointcloud(input_points: np.array, classes: list):
-    """Filter pointcloud and extracts point coordinates (X, Y, Z)
+    """Filter pointcloud by class
 
     Args:
         input_points (np.array): Numpy array containing point coordinates (X, Y, Z, classification)
         classes (list): List of classes to use for the filtering
 
     Returns:
-        filtered_points (np.ndarray) : Numpy array containing point coordinates (X, Y, Z) filtering
+        filtered_points (np.ndarray) : Numpy array containing point coordinates (X, Y, Z) after filtering
     """
     # Filter pointcloud by classe(s)
     points_mask = np.isin(input_points[:, -1], classes)

lidro/pointcloud/read_las.py → lidro/create_mask_hydro/pointcloud/read_las.py

@@ -11,7 +11,8 @@ def read_pointcloud(las_file: str):
         las_file (str): Path to the LAS file
 
     Returns:
-        input_point (np.ndarray) : Numpy array containing point coordinates (X, Y, Z, classification)
+        input_point (np.ndarray) : Numpy array containing point coordinates and classification
+        (X, Y, Z, classification)
         crs (dict): a pyproj CRS object
     """
     # Read pointcloud

lidro/create_mask_hydro/rasters/create_mask_raster.py

@@ -0,0 +1,76 @@
+# -*- coding: utf-8 -*-
+""" Create mask non hydro from pointcloud filtering
+"""
+from typing import List, Tuple
+
+import numpy as np
+import scipy.ndimage
+
+from lidro.create_mask_hydro.pointcloud.filter_las import filter_pointcloud
+from lidro.create_mask_hydro.pointcloud.io import get_pointcloud_origin
+from lidro.create_mask_hydro.pointcloud.read_las import read_pointcloud
+
+
+def create_occupancy_map(points: np.array, tile_size: int, pixel_size: float, origin: Tuple[int, int]):
+    """Create a binary image to extract water surfaces
+
+    Args:
+        points (np.array): array from pointcloud
+        tile_size (int): size of the raster grid (in meters)
+        pixel_size (float): distance between each node of the raster grid (in meters)
+        origin (Tuple[int, int]): Coordinates of the top left corner of the top-left pixel
+
+    Returns:
+        bins (np.array): 2D binary array (x, y) with 1 for if there is at least one point of the input cloud
+        in the corresponding pixel, 0 otherwise
+    """
+    # Compute number of points per bin
+    bins_x = np.arange(origin[0], origin[0] + tile_size + pixel_size, pixel_size)
+    bins_y = np.arange(origin[1] - tile_size, origin[1] + pixel_size, pixel_size)
+    _bins, _, _ = np.histogram2d(points[:, 1], points[:, 0], bins=[bins_y, bins_x])
+    bins = np.flipud(_bins)
+    bins = bins > 0  # binarize map
+    return bins
+
+
+def detect_hydro_by_tile(filename: str, tile_size: int, pixel_size: float, classes: List[int], dilation_size: int):
+    """ "Detect hydrographic surfaces in a tile from the classified points of the input pointcloud
+        An hydrographic surface is defined as a surface where there is no points from any class different from water
+    The output hydrographic surface mask is dilated to make sure that the masks are continuous when merged with their
+    neighbours
+        Args:
+            filename (str): input pointcloud
+            tile_size (int): size of the raster grid (in meters)
+            pixel_size (float): distance between each node of the raster grid (in meters)
+            classes (List[int]): List of classes to use for the binarisation (points with other
+                        classification values are ignored)
+            dilation_size (int): size of the structuring element for dilation
+
+        Returns:
+            smoothed_water (np.array):  2D binary array (x, y) of the water presence from the point cloud
+            pcd_origin (list): top left corner of the tile containing the point cloud
+            (infered from pointcloud bounding box and input tile size)
+    """
+    # Read pointcloud, and extract coordinates (X, Y, Z, and classification) of all points
+    array, crs = read_pointcloud(filename)
+
+    # Extracts parameters for binarisation
+    pcd_origin = get_pointcloud_origin(array, tile_size)
+
+    # Filter pointcloud by classes: keep only points from selected classes that are not water
+    array_filter = filter_pointcloud(array, classes)
+
+    # create occupancy map (2D)
+    occupancy = create_occupancy_map(array_filter, tile_size, pixel_size, pcd_origin)
+
+    # Revert occupancy map to keep pixels where there is no point of the selected classes
+    detected_water = np.logical_not(occupancy)
+
+    # Apply a mathematical morphology operations: DILATION
+    # / ! \ NOT "CLOSING", due to the reduction in the size of hydro masks (tile borders)
+    # / ! \ WITH "CLOSING" => Masks Hydro are no longer continuous, when they are merged
+    water_mask = scipy.ndimage.binary_dilation(
+        detected_water, structure=np.ones((dilation_size, dilation_size))
+    ).astype(np.uint8)
+
+    return water_mask, pcd_origin

lidro/create_mask_hydro/vectors/convert_to_vector.py

@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+""" Vectorize
+"""
+import geopandas as gpd
+import numpy as np
+from rasterio.features import shapes as rasterio_shapes
+from rasterio.transform import from_origin
+from shapely.geometry import shape as shapely_shape
+
+from lidro.create_mask_hydro.rasters.create_mask_raster import detect_hydro_by_tile
+
+
+def create_hydro_vector_mask(
+    filename: str, output: str, pixel_size: float, tile_size: int, classes: list, crs: str, dilatation_size: int
+):
+    """Create a vector mask of hydro surfaces in a tile from the points classification of the input LAS/LAZ file,
+    and save it as a GeoJSON file.
+
+
+    Args:
+        filename (str): path to the input pointcloud
+        output (str): path to output geoJSON
+        pixel_size (float): distance between each node of the intermediate raster grid (in meters)
+        tile_size (int): size of the intermediate raster grid (in meters)
+        classes (list): List of classes to consider as water (points with other classification values are ignored)
+        crs (str): a pyproj CRS object used to create the output GeoJSON file
+        dilatation_size (int): size for dilatation raster
+    """
+    # Read a binary image representing hydrographic surface(s)
+    binary_image, pcd_origin = detect_hydro_by_tile(filename, tile_size, pixel_size, classes, dilatation_size)
+
+    # Extract origin
+    origin_x = pcd_origin[0]
+    origin_y = pcd_origin[1]
+    # Calculate "transform"
+    transform = from_origin(origin_x, origin_y, pixel_size, pixel_size)
+
+    # Convert binary image to vector
+    geometry = [
+        shapely_shape(shapedict)
+        for shapedict, value in rasterio_shapes(
+            binary_image.astype(np.int16), mask=None, connectivity=8, transform=transform
+        )
+        if value != 0
+    ]
+    # save the result
+    gdf = gpd.GeoDataFrame(geometry=geometry, crs=crs)
+    gdf.to_file(output, driver="GeoJSON", crs=crs)