Rectify model under bridge

README.md

@@ -59,7 +59,13 @@ 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écessitant plusieurs étapes intermédiaires :
+* 5- Vérifier que les modèles de régréssion linéaire calculés le long du squelette s'écoulement progressivement le long du cours d'eau, c'est-à-dire éviter les zones de cuvettes sous les ponts. D'un point de vue mathématique, cela signifie que le "Dernier point squelette AMONT < Premier point squelette AVAL = ALERTE".
+![Alerte](images/alerte_squelette.png)
+
+Pour éviter ces "alertes", il faut corriger les valeurs Z des N points du squelette jusqu'à le Z du squelette aval est égal au dernier point Z du squelette amont.
+![Corrections des alertes](images/correction_alerte_squelette.png)
+
+* 6- Création de points virtuels nécessitant 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)
 
@@ -259,6 +265,7 @@ Options généralement passées en paramètres :
 * io.input_mask_hydro : Le chemin contenant le masque HYDRO fusionné (ex."./data/merge_mask_hydro/MaskHydro_merge.geosjon").
 * io.input_skeleton= Le chemin contenant le squelette hydrographique (ex. "./data/skeleton_hydro/Skeleton_Hydro.geojson")
 * io.dir_points_skeleton : Le chemin contenant l'ensemble des N points du squelette créés à l'échelle des dalles LIDAR ( ex. "./tmp/point_skeleton/").
+* io.input_bridge : Le chemin contenant le tablier de pont produit par la production (ex. "./data/bridge/tablier_pont.geojson")
 * io.output_dir :  Le chemin du dossier de sortie (les points virtuels à l'échelle du projet).
 
 Autres paramètres disponibles :

configs/configs_lidro.yaml

@@ -9,6 +9,7 @@ work_dir: ${hydra:runtime.cwd}
 io:
   input_filename: null
   input_mask_hydro: null
+  input_bridge: null
   input_skeleton: null
   input_dir: null
   input_dir_point_virtual: null
@@ -33,7 +34,7 @@ mask_generation:
     dilation_size: 3
   filter:
     # Classes to be considered as "non-water"
-    keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67] # All classes
+    keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67, 69] # All classes
   vector:
     # Filter water's area (m²)
     min_water_area: 150
@@ -68,7 +69,7 @@ skeleton:
 virtual_point:
   filter:
     # Keep ground and water pointclouds between Hydro Mask and Hydro Mask buffer
-    keep_neighbors_classes: [2, 9]
+    keep_neighbors_classes: [2, 69]
   vector:
     # Distance in meters between 2 consecutive points from Skeleton Hydro
     distance_meters: 5

lidro/create_virtual_point/vectors/create_skeleton_3d.py

@@ -0,0 +1,107 @@
+# -*- coding: utf-8 -*-
+""" Run function "update skeleton with Z"
+"""
+import logging
+
+import geopandas as gpd
+import numpy as np
+from shapely import line_locate_point
+from shapely.geometry import LineString, Point
+
+
+def skeleton_3d_with_linear_regression(line: gpd.GeoDataFrame, model: np.poly1d, crs: str) -> gpd.GeoDataFrame:
+    """
+    This function updates the 2D hydro skeleton 'line' by adding a 3D 'Z' value
+    based on a previously calculated linear regression model, and saves the updated
+    skeleton to a GeoJSON file.
+
+    Args:
+        line (gpd.GeoDataFrame): A GeoDataFrame containing each 2D line from Hydro's Skeleton.
+        model (np.poly1d): The linear regression model to calculate Z values.
+        crs (str): The coordinate reference system of the line.
+
+    Returns:
+        gpd.GeoDataFrame: A GeoDataFrame of the updated 3D skeleton with Z values.
+    """
+    # For each line in the GeoDataFrame, extract points and calculate Z values
+    updated_geometries = []
+    z_mean_values = []  # List to store Z mean values for each geometry
+
+    for i, geom in line.iterrows():
+        line_geom = geom["geometry"]
+
+        # If the line is not a LineString, skip it
+        if line_geom.geom_type != "LineString":
+            logging.warning(f"Geometry at index {i} is not a LineString")
+            continue
+
+        # Extract points (2D)
+        coords_2d = np.array(line_geom.coords)
+
+        # Calculate the curvilinear abscissa for each point using line_locate_point
+        abscissas = abscissas = np.array([line_locate_point(line_geom, Point(x, y)) for x, y in coords_2d])
+
+        # Predict Z values using the regression model
+        z_values = model(abscissas)
+
+        # Calculate the mean Z value and store it
+        z_mean = np.mean(z_values)
+        z_mean_values.append(z_mean)
+
+        # Update each coordinate with the Z value
+        coords_3d = [(x, y, z) for (x, y), z in zip(coords_2d, z_values)]
+
+        # Create a new LineString with 3D coordinates
+        updated_geom = LineString(coords_3d)
+        updated_geometries.append(updated_geom)
+
+    # Create a new GeoDataFrame for the updated lines
+    updated_line = gpd.GeoDataFrame(geometry=updated_geometries, crs=crs)
+
+    # Add the Z mean values as a new column
+    updated_line["z_mean"] = z_mean_values
+
+    return updated_line
+
+
+def skeleton_3d_with_flatten(line: gpd.GeoDataFrame, z_value: float, crs: str) -> gpd.GeoDataFrame:
+    """
+    This function updates the 2D hydro skeleton 'line' by adding a 3D 'Z' value
+    based on a previously calculated flatten model, and saves the updated
+    skeleton to a GeoJSON file.
+
+    Args:
+        line (gpd.GeoDataFrame): A GeoDataFrame containing each 2D line from Hydro's Skeleton.
+        z_value (float) : Z of the river
+        crs (str): The coordinate reference system of the line.
+
+    Returns:
+        gpd.GeoDataFrame: A GeoDataFrame of the updated 3D skeleton with Z values.
+    """
+    # For each line in the GeoDataFrame, extract points and calculate Z values
+    updated_geometries = []
+    for i, geom in line.iterrows():
+        line_geom = geom["geometry"]
+
+        # If the line is not a LineString, skip it
+        if line_geom.geom_type != "LineString":
+            logging.warning(f"Geometry at index {i} is not a LineString")
+            continue
+
+        # Extract points (2D)
+        coords_2d = np.array(line_geom.coords)
+
+        # Update each coordinate with the Z value
+        coords_3d = [(x, y, z) for (x, y), z in zip(coords_2d, z_value)]
+
+        # Create a new LineString with 3D coordinates
+        updated_geom = LineString(coords_3d)
+        updated_geometries.append(updated_geom)
+
+    # Create a new GeoDataFrame for the updated lines
+    updated_line = gpd.GeoDataFrame(geometry=updated_geometries, crs=crs)
+
+    # Add the Z mean values as a new column
+    updated_line["z_mean"] = z_value
+
+    return updated_line

lidro/create_virtual_point/vectors/intersect_skeleton_by_bridge.py

@@ -0,0 +1,103 @@
+# -*- coding: utf-8 -*-
+""" Run function "extract skeleton by bridge"
+"""
+from typing import List
+
+import pandas as pd
+from shapely.geometry import LineString, Point, Polygon
+from shapely.ops import nearest_points
+
+
+def sort_skeleton_points_by_proximity(skeleton: LineString, reference_geometry: Polygon, is_upstream: bool):
+    """
+    Sort the points of a LineString skeleton based on their proximity to the edge of the reference geometry.
+
+    Args:
+        skeleton (LineString): Skeleton geometry to sort.
+        reference_geometry (Polygon): Reference geometry (bridge) to sort the skeleton's points by proximity.
+        is_upstream (bool): Boolean indicating if the skeleton is upstream or downstream.
+                            If True, the farthest point from the bridge is first (upstream).
+                            If False, the closest point to the bridge is first (downstream).
+
+    Returns:
+        LineString: A sorted LineString with points ordered based on whether the skeleton is upstream or downstream.
+    """
+    # Find nearest point on the skeleton to the bridge
+    reference_point = nearest_points(reference_geometry, skeleton)[0]
+    points_with_distance = [(point, Point(point).distance(reference_point)) for point in skeleton.coords]
+
+    # Sort points by their distance to the reference point
+    points_with_distance.sort(key=lambda x: x[1])
+    sorted_points = [point for point, distance in points_with_distance]
+
+    # If upstream, reverse the order so that the farthest point is first
+    if is_upstream:
+        sorted_points = sorted_points[::-1]
+
+    return LineString(sorted_points)
+
+
+def extract_bridge_skeleton_info(bridge: Polygon, list_skeleton_with_z: List):
+    """
+    Extract intersection information between bridges and skeletons.
+    For each intersection, retrieve the bridge geometry, the geometry of the upstream (amont) skeleton,
+    the geometry of the downstream (aval) skeleton, the Z value of the last point of the upstream skeleton,
+    and the Z value of the first point of the downstream skeleton.
+
+    Args:
+        bridge (Polygon): bridge geometries.
+        list_skeleton_with_z (list): List of skeletons with Z values and geometries.
+
+    Returns:
+        pd.DataFrame: DataFrame containing bridge geometry, upstream skeleton geometry,
+                      downstream skeleton geometry, Z value of the last upstream point,
+                      and Z value of the first downstream point.
+    """
+    intersecting_skeletons = [
+        skeleton for skeleton in list_skeleton_with_z if bridge.intersects(skeleton["geometry"]).any()
+    ]
+    list_skeleton = [i for i in intersecting_skeletons if len(intersecting_skeletons) == 2]
+    # Check if exactly two skeletons intersect the bridge
+    if round(list_skeleton[0]["z_mean"].item(), 2) > round(list_skeleton[1]["z_mean"].item(), 2):
+        skeleton_upstream = list_skeleton[0]
+        skeleton_downstream = list_skeleton[1]
+    else:
+        skeleton_upstream = list_skeleton[1]
+        skeleton_downstream = list_skeleton[0]
+
+    # Sort the points for upstream and downstream skeletons based on proximity to the bridge
+    geometry_upstream_sorted = sort_skeleton_points_by_proximity(
+        skeleton_upstream["geometry"].iloc[0], bridge, is_upstream=True
+    )
+    geometry_downstream_sorted = sort_skeleton_points_by_proximity(
+        skeleton_downstream["geometry"].iloc[0], bridge, is_upstream=False
+    )
+
+    # Extract the Z value of the last point of the upstream skeleton
+    if isinstance(geometry_upstream_sorted, LineString):
+        last_point_upstream = list(geometry_upstream_sorted.coords)[-1]  # Get the coordinates of the last point
+        z_upstream = last_point_upstream[2]  # The Z value of the last point
+    else:
+        raise ValueError("Upstream skeleton geometry is not a LineString.")
+
+    # Extract the Z value of the first point of the downstream skeleton
+    if isinstance(geometry_downstream_sorted, LineString):
+        first_point_downstream = list(geometry_downstream_sorted.coords)[0]  # Get the coordinates of the first point
+        z_downstream = first_point_downstream[2]  # The Z value of the first point
+    else:
+        raise ValueError("Downstream skeleton geometry is not a LineString.")
+
+    # Extract alert if Z value of the upstream skeleton < Z value of the downstream skeleton
+    alert = z_upstream < z_downstream
+
+    # Create a DataFrame to return the results
+    data = {
+        "bridge_geometry": [bridge],
+        "upstream_geometry": [geometry_upstream_sorted],
+        "downstream_geometry": [geometry_downstream_sorted],
+        "z_upstream": [z_upstream],
+        "z_downstream": [z_downstream],
+        "alert": [alert],
+    }
+
+    return pd.DataFrame(data)