Merge pull request #10 from IGNF/regression_model_from_skeleton

Regression model from skeleton

README.md

@@ -8,40 +8,64 @@ Pour créer des modèles numériques cohérents avec les modèles hydrologiques,
 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.
+* Mise en plan des grands cours d’eau (>5m large) pour assurer l’écoulement​. A noter que pour l'instant seulement cette étape est développée.
 
 ## Traitement
+
+![Chaine de traitement global de LIDRO](images/chaine_traitement_lidro.jpg)
+
 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.)
+- données vectorielles représentant le réseau hydrographique issu des différentes bases de données IGN (BDUnis, BDTopo, etc.) / ! \ requête SQL fait manuellement pour le moment. Exemple de requête SQL sur le bloc PM du LIDAR HD
+```
+WITH emprise_PM AS (SELECT st_GeomFromText('POLYGON((875379.222972973 6431750.0,
+875379.222972973 6484250.0,
+946620.777027027 6484250.0,
+946620.777027027 6431750.0,
+875379.222972973 6431750.0))') as geom)
+
+SELECT geometrie, nature, fictif, persistance, classe_de_largeur, position_par_rapport_au_sol
+FROM troncon_hydrographique
+JOIN emprise_PM ON st_intersects(troncon_hydrographique.geometrie,emprise_PM.geom) 
+WHERE NOT gcms_detruit 
+AND classe_de_largeur NOT IN ('Entre 0 et 5 m', 'Sans objet') 
+AND position_par_rapport_au_sol='0'
+```
 
-Trois grands axes du processus à mettre en place en distanguant l'échelle de traitmeent associé :
+Trois grands axes du processus à mettre en place en distanguant l'échelle de traitement 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 < 150 m² (paramétrable)
   * la suppression des aires < 1000 m² hors BD IGN (grands cours d'eau < 5m de large)
+A noter que pour l'instant, la suppresion des masques hydrographiques en dehors des grands cours d'eau et le nettoyage de ces masques s'effectuent manuellement. Ce processus sera développé prochainement en automatique.
 * 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.)
+  * Surfaces planes (mer, lac, étang, etc.) (pas encore développé)
 
-### Traitement des grands cours d'eau (> 5 m de large dans la BD Uns).
+
+### Traitement des grands cours d'eau (> 5 m de large dans la BD Unis).
+![Chaine de traitement des points virtuels](images/process_points_virtuels.jpg)
 
 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
+* 1- création automatique du tronçon hydrographique ("Squelette hydrographique", soit les tronçons hydrographiques dans la BD Unis) à partir de l'emprise du masque hydrographique "écoulement".
+/ ! \ EN AMONT : Appariement, contrôle et correction manuels des masques hydrographiques "écoulement" (rivières) et du squelette hydrographique
+
+A l'échelle de l'entité hydrographique : 
+* 2- Réccupérer tous les points LIDAR considérés comme du "SOL" situés à la limite de berges (masque hydrographique) moins N mètres
+Pour les cours d'eau supérieurs à 150 m de long :
+* 3- Transformer les coordonnées de ces points (étape précédente) en abscisses curvilignes
+* 4- Générer un modèle de régression linéaire afin de générer tous les N mètres une valeur d'altitude le long du squelette de cette rivière. Les différents Z le long des squelettes HYDRO doivent assurer l'écoulement. Il est important de noter que tous les 50 mètres semble une valeur correcte pour appréhender la donnée. Cette valeur s'explique en raison de la précision altimétrique des données LIDAR (20 cm) ET que les rivières françaises correspondent à des cours d’eau naturels dont la pente est inférieure à 1%. 
+/ ! \ Pour les cours d'eau inférieurs à 150 m de long, le modèle de régression linéaire ne fonctionne pas. La valeur du premier quartile sera calculée sur l'ensemble des points d'altitudes du LIDAR "SOL" (étape 2) et affectée pour ces entités hydrographiques (< 150m de long) : aplanissement. 
+* 5- Création de points virtuels nécéssitant plusieurs étapes intermédiaires :
+  * Création des points virtuels 2D espacés selon une grille régulière tous les N mètres (paramétrable) à l'intérieur du masque hydrographique "écoulement"
+  * Affecter une valeur d'altitude à ces points virtuels en fonction des "Z" calculés à l'étape précédente (interpolation linéaire ou aplanissement)
 
 ### Traitement des surfaces planes (mer, lac, étang, etc.)
 Pour rappel, l'eau est considérée comme horizontale sur ce type de surface.
-
+/ ! \ EN ATTENTE / ! \
 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.
+* 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 l'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
@@ -93,14 +117,35 @@ Voir les tests fonctionnels en bas du README.
 
 ## Tests
 ### Tests fonctionnels
-Tester sur un seul fichier LAS/LAZ
+* 1- Créer des masques hydrographiques à l'échelle de la dalle LIDAR
+Tester sur un seul fichier LAS/LAZ pour créer un/des masques hydrographiques sur une dalle LIDAR
+```
+example_create_mask_by_tile.sh
+```
+Tester sur un dossier contenant plusieurs dalles LIDAR pour créer un/des masques hydrographiques 
+```
+example_create_mask_default.sh
+```
+* 2- Créer un masque hydrographiques fusionné et nettoyé à l'échelle de l'ensemble de l'ensemble des dalles LIDAR
+Tester sur un dossier contenant plusieurs dalles LIDAR pour créer fusionner l'ensemble des masques hydrographiques 
+```
+example_merge_mask_default.sh
+```
+* 3- Création des tronçons hydrographiques à l'échelle de/des entité(s) hydrographique(s)
+```
+example_create_skeleton_lines.sh
+```
+
+* 4- Création des points virtuels
+
+A. Tester sur un dossier contenant plusieurs dalles LIDAR pour créer des points tous les N mètres le long du squelette hydrographique, et réccupérer les N plus proches voisins points LIDAR "SOL"
 ```
-example_lidro_by_tile.sh
+example_extract_points_around_skeleton_default.sh
 ```
 
-Tester sur un dossier
+B. Tester sur un dossier contenant plusieurs dalles LIDAR pour créer des points virtuels 3D à l'intérieurs des masques hydrographiques 
 ```
-example_lidro_default.sh
+example_create_virtual_point_default.sh
 ```
 
 ### Tests unitaires

configs/configs_lidro.yaml

@@ -12,6 +12,7 @@ io:
   input_skeleton: null
   input_dir: null
   output_dir: null
+  dir_points_skeleton: null # folder to contains files (.GeoJSON) by LIDAR tiles : neighboring points for each point of the skeleton
   srid: 2154
   extension: .tif
   raster_driver: GTiff
@@ -27,7 +28,7 @@ io:
 
 mask_generation:
   raster:
-    # size for dilatation
+    # Size for dilatation
     dilation_size: 3
   filter:
     # Classes to be considered as "non-water"
@@ -41,18 +42,6 @@ mask_generation:
     # Tolerance from Douglas-Peucker
     tolerance: 1
 
-virtual_point:
-  filter:
-    # Keep ground and water pointclouds between Hydro Mask and Hydro Mask buffer
-    keep_neighbors_classes: [2, 9]
-  vector:
-    # distance in meters between 2 consecutive points 
-    distance_meters: 2
-    # buffer for searching the points classification (default. "2") of the input LAS/LAZ file
-    buffer: 2
-    # The number of nearest neighbors to find with KNeighbors
-    k: 20
-
 skeleton:
   max_gap_width: 200 # distance max in meter of any gap between 2 branches we will try to close with a line
   max_bridges: 1 # number max of bridges that can be created between 2 branches
@@ -74,6 +63,28 @@ skeleton:
     water_min_size: 500 # min size of a skeleton line to be sure not to be removed (should be at least more than half the max river width)
     max_gap_candidates: 3 # max number of candidates to close a gap between 2 branches
 
+virtual_point:
+  filter:
+    # Keep ground and water pointclouds between Hydro Mask and Hydro Mask buffer
+    keep_neighbors_classes: [2, 9]
+  vector:
+    # Distance in meters between 2 consecutive points from Skeleton Hydro
+    distance_meters: 5
+    # Buffer for searching the points classification (default. "3") of the input LAS/LAZ file
+    buffer: 3
+    # The number of nearest neighbors to find with KNeighbors
+    k: 3
+    # The minimum length of a river to use the linear regression model
+    river_length: 150 
+  pointcloud:
+    # Spacing between the grid points in meters by default "0.5"
+    points_grid_spacing: 0.5
+    # The number of the classe assign those virtual points
+    virtual_points_classes : 68
+
+
+
+
 hydra:
   output_subdir: null
   run:

lidro/create_virtual_point/pointcloud/convert_list_points_to_las.py

@@ -0,0 +1,46 @@
+# -*- coding: utf-8 -*-
+""" this function convert GeoDataframe "virtual points" to LIDAR points
+"""
+import logging
+import os
+from typing import List
+
+import laspy
+import numpy as np
+import pandas as pd
+from shapely.geometry import Point
+
+
+def list_points_to_las(virtual_points: List[Point], output_dir: str, crs: str, virtual_points_classes: int):
+    """This function convert virtual points (List of virtuals points) to LIDAR points
+       with classification for virtual points
+
+    Args:
+        virtual_points (List[Point]): A list containing virtuals points by hydrological entity
+        output_dir (str): folder to output LAS
+        crs (str): a pyproj CRS object used to create the output GeoJSON file
+        virtual_points_classes (int): The classification value to assign to those virtual points
+    """
+    # Combine all virtual points into a single GeoDataFrame
+    grid_with_z = pd.concat(virtual_points, ignore_index=True)
+
+    # Create a LAS file with laspy
+    header = laspy.LasHeader(point_format=6, version="1.4")
+
+    # Create a LAS file with laspy and add the VLR for CRS
+    las = laspy.LasData(header)
+    las.x = grid_with_z.geometry.x
+    las.y = grid_with_z.geometry.y
+    las.z = grid_with_z.geometry.z
+    las.classification = virtual_points_classes * np.ones(len(grid_with_z.index))
+
+    # Add CRS information to the VLR
+    vlr = laspy.vlrs.known.WktCoordinateSystemVlr(crs.to_wkt())
+    las.vlrs.append(vlr)
+
+    # Save the LAS file
+    output_laz = os.path.join(output_dir, "virtual_points.laz")
+    with laspy.open(output_laz, mode="w", header=las.header, do_compress=True) as writer:
+        writer.write_points(las.points)
+
+    logging.info(f"Virtual points LAS file saved to {output_laz}")

lidro/create_virtual_point/stats/knn_distance.py

@@ -32,6 +32,5 @@ def find_k_nearest_neighbors(points_skeleton: np.array, points_ground_3d: np.arr
     # Retrieve the points corresponding to the indices
     # Use indices  to retrieve the closest 3D points
     k_nearest_points = points_ground_3d[indices[0]]
-    # = points_3d[indices.flatten()]
 
     return k_nearest_points

lidro/create_virtual_point/vectors/apply_Z_from_grid.py

@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+""" Apply Z from grid
+"""
+import geopandas as gpd
+from shapely import line_locate_point
+
+
+def calculate_grid_z_with_model(points: gpd.GeoDataFrame, line: gpd.GeoDataFrame, model) -> gpd.GeoDataFrame:
+    """Calculate Z with regression model.
+    If points are not on the line, these points will be projected on this line
+
+    Args:
+        points (gpd.GeoDataFrame): A GeoDataFrame containing the grid points
+        line (gpd.GeoDataFrame): A GeoDataFrame containing each line from Hydro's Skeleton
+        model (dict): A dictionary representing the regression model
+
+    Returns:
+        gpd.GeoDataFrame: A GeoDataFrame of initial points, but with a Z.
+    """
+    # Calculate curvilinear abscises for all points of the grid
+    curvilinear_abs = line_locate_point(line.loc[0, "geometry"], points["geometry"].array, normalized=False)
+
+    # Prediction of Z values using the regression model
+    # Its possible to use non-linear models for prediction
+    predicted_z = model(curvilinear_abs)
+
+    # Generate a new geodataframe, with 3D points
+    grid_with_z = calculate_grid_z(points, predicted_z)
+
+    return grid_with_z
+
+
+def calculate_grid_z(points: gpd.GeoDataFrame, predicted_z: float) -> gpd.GeoDataFrame:
+    """Calculate Z grid
+
+    Args:
+        points (gpd.GeoDataFrame): A GeoDataFrame containing the grid points
+        predicted_z (float): predicted Z for river
+
+    Returns:
+        gpd.GeoDataFrame: A GeoDataFrame of initial points, but with a Z.
+    """
+    # Generate a new geodataframe, with 3D points
+    grid_with_z = gpd.GeoDataFrame(
+        geometry=gpd.GeoSeries().from_xy(points.geometry.x, points.geometry.y, predicted_z), crs=points.crs
+    )
+
+    return grid_with_z

lidro/create_virtual_point/vectors/create_grid_2D_inside_maskhydro.py

@@ -0,0 +1,46 @@
+# -*- coding: utf-8 -*-
+""" Create a regula 2D grid
+"""
+import geopandas as gpd
+import numpy as np
+import shapely.vectorized
+
+
+def generate_grid_from_geojson(mask_hydro: gpd.GeoDataFrame, spacing: float = 0.5, margin: float = 0):
+    """
+    Generates a regular 2D grid of evenly spaced points within a polygon defined
+    in a GeoJSON file.
+
+    Args:
+        mask_hydro (gpd.GeoDataFrame): A GeoDataFrame containing each mask hydro from Hydro's Skeleton
+        spacing (float, optional): Spacing between the grid points in meters. The default value is 0.5 meter.
+        margin (int, optional): Margin around mask for grid creation. The default value is 0 meter.
+
+    Returns:
+        geopandas.GeoDataFrame: A GeoDataFrame containing the grid points within the polygon.
+    """
+    # Extract the polygon
+    polygon = mask_hydro.geometry.unary_union  # Combine all polygons into a single shape if there are multiple
+
+    if margin:
+        polygon = polygon.buffer(margin)
+
+    # Get the bounds of the polygon
+    minx, miny, maxx, maxy = polygon.bounds
+
+    # Generate the grid points
+    x_points = np.arange(minx, maxx, spacing)
+    y_points = np.arange(miny, maxy, spacing)
+    grid_points = np.meshgrid(x_points, y_points)
+
+    # Filter points that are within the polygon
+    grid_points_in_polygon = shapely.vectorized.contains(polygon, *grid_points)
+
+    filtred_x = grid_points[0][grid_points_in_polygon]
+    filtred_y = grid_points[1][grid_points_in_polygon]
+
+    # Convert to GeoDataFrame
+    points_gs = gpd.points_from_xy(filtred_x, filtred_y)
+    grid_gdf = gpd.GeoDataFrame(geometry=points_gs, crs=mask_hydro.crs)
+
+    return grid_gdf