diff --git a/README.md b/README.md index 80d1b3a2..218b98a9 100644 --- a/README.md +++ b/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 : diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index e55ebb56..05799cb0 100644 --- a/configs/configs_lidro.yaml +++ b/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 diff --git a/images/alerte_squelette.png b/images/alerte_squelette.png new file mode 100644 index 00000000..4b820e69 Binary files /dev/null and b/images/alerte_squelette.png differ diff --git a/images/correction_alerte_squelette.png b/images/correction_alerte_squelette.png new file mode 100644 index 00000000..5291a4e0 Binary files /dev/null and b/images/correction_alerte_squelette.png differ diff --git a/images/lidro_correction_pont.odp b/images/lidro_correction_pont.odp new file mode 100644 index 00000000..f3ebc3cd Binary files /dev/null and b/images/lidro_correction_pont.odp differ diff --git a/lidro/create_virtual_point/vectors/apply_Z_from_grid.py b/lidro/create_virtual_point/vectors/apply_Z_from_grid.py deleted file mode 100644 index 9e6d497b..00000000 --- a/lidro/create_virtual_point/vectors/apply_Z_from_grid.py +++ /dev/null @@ -1,48 +0,0 @@ -# -*- 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 diff --git a/lidro/create_virtual_point/vectors/create_skeleton_3d.py b/lidro/create_virtual_point/vectors/create_skeleton_3d.py new file mode 100644 index 00000000..ce03f959 --- /dev/null +++ b/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 diff --git a/lidro/create_virtual_point/vectors/intersect_skeleton_by_bridge.py b/lidro/create_virtual_point/vectors/intersect_skeleton_by_bridge.py new file mode 100644 index 00000000..3b8b29c6 --- /dev/null +++ b/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) diff --git a/lidro/create_virtual_point/vectors/rectify_alert.py b/lidro/create_virtual_point/vectors/rectify_alert.py new file mode 100644 index 00000000..cdeb4118 --- /dev/null +++ b/lidro/create_virtual_point/vectors/rectify_alert.py @@ -0,0 +1,131 @@ +# -*- coding: utf-8 -*- +"""Optimized function to adjust Z values of the downstream skeleton""" +import logging +from typing import List, Tuple + +import geopandas as gpd +import numpy as np +import pandas as pd +from shapely.geometry import CAP_STYLE, LineString + +from lidro.create_virtual_point.vectors.intersect_skeleton_by_bridge import ( + extract_bridge_skeleton_info, +) + + +def rectify_z_from_downstream_skeleton( + input_bridge: str, list_skeleton_with_z: List, crs: str +) -> Tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]: + """ + Adjusts the Z values of N points in the downstream skeleton to match the Z value + of the last point in the upstream skeleton, up until reaching z_upstream + 0.1. + If no alert is active, adds the original downstream geometry without modification. + + Args: + input_bridge (str): Path to the bridge input file + list_skeleton_with_z (List): List of skeletons with Z values and associated hydro mask + crs (str): Coordinate reference system + + Returns: + Tuple[gpd.GeoDataFrame, gpd.GeoDataFrame]: A tuple of GeoDataFrames containing + the updated or original downstream and upstream skeleton geometries. + """ + # Load bridge data and apply a 5m buffer to each geometry + gdf_bridge = gpd.read_file(input_bridge, crs=crs) + gdf_bridge["geometry"] = gdf_bridge.buffer(5, cap_style=CAP_STYLE.square) + + downstream_geometries = [] + upstream_geometries = [] + + for bridge in gdf_bridge["geometry"]: + skeleton_info = extract_bridge_skeleton_info(bridge, list_skeleton_with_z) + + downstream_geometry = skeleton_info["downstream_geometry"] + upstream_geometry = skeleton_info["upstream_geometry"] + + # Add upstream geometry directly without modification + upstream_geometries.extend(upstream_geometry) + + # Check if alert is active and process downstream geometry accordingly + if skeleton_info["alert"].iloc[0]: + logging.info("Alert detected for bridge intersection. Updating Z values...") + z_upstream = float(skeleton_info["z_upstream"].iloc[0]) + z_downstream = float(skeleton_info["z_downstream"].iloc[0]) + + for line in downstream_geometry: + # Convert coordinates to NumPy arrays for efficient Z-value handling + coords = np.array(line.coords) + x_vals, y_vals, z_vals = coords[:, 0], coords[:, 1], coords[:, 2] + + # Set the initial Z value of the downstream skeleton to z_downstream + z_vals[0] = z_downstream + + # Apply z_upstream to all points in the downstream skeleton until reaching z_upstream + z_vals[z_vals > z_upstream] = z_upstream + + # Rebuild the geometry with the updated Z values + updated_coords = np.column_stack([x_vals, y_vals, z_vals]) + downstream_geometries.append(LineString(updated_coords)) + + if not skeleton_info["alert"].iloc[0]: + # If no alert, add the original downstream geometry + logging.info("No alert detected. Adding original downstream geometry.") + downstream_geometries.extend(downstream_geometry) + + # Create GeoDataFrames for downstream and upstream geometries + downstream_gdf = gpd.GeoDataFrame([{"geometry": geom} for geom in downstream_geometries], crs=crs) + + upstream_gdf = gpd.GeoDataFrame([{"geometry": geom} for geom in upstream_geometries], crs=crs) + + return downstream_gdf, upstream_gdf + + +def create_list_rectify_skeleton_with_mask( + input_bridge: str, hydro_mask_path: str, list_skeleton_with_z: List, crs: str +) -> List[gpd.GeoDataFrame]: + """ + Processes skeletons per hydro mask, prioritizes downstream geometries, and returns + a list of GeoDataFrames containing unique skeleton geometries intersected with each hydro mask, + along with the geometry of each mask. + + Args: + input_bridge (str): Path to the bridge input file + list_skeleton_with_z (List): List of skeletons with Z values and associated hydro masks + crs (str): Coordinate reference system + hydro_mask_path (str): Path to the hydro mask GeoJSON file + + Returns: + List[gpd.GeoDataFrame]: A list of GeoDataFrames, each containing a unique skeleton geometry + intersected with its hydro mask and the geometry of the mask itself. + """ + # # Load the hydro mask GeoJSON as a GeoDataFrame + hydro_masks = gpd.read_file(hydro_mask_path).to_crs(crs) + + # # Generate downstream and upstream skeletons + downstream_gdf, upstream_gdf = rectify_z_from_downstream_skeleton(input_bridge, list_skeleton_with_z, crs) + + # # Identify unique downstream geometries and add unmatched upstream geometries + # Keep Only the Matching Geometries from downstream_gdf + result_gdf = downstream_gdf.drop_duplicates(subset="geometry") + # Add Non-Matching Geometries from upstream_gdf + unmatched_upstream = upstream_gdf[~upstream_gdf["geometry"].isin(downstream_gdf["geometry"])] + result_gdf = pd.concat([result_gdf, unmatched_upstream], ignore_index=True) + + # # For each mask, get the intersecting skeleton geometries and prioritize downstream geometries + list_result_gdfs = [] + for mask in hydro_masks.itertuples(): + # Filter result_gdf to only include geometries that intersect with the current mask + intersecting_geometries = result_gdf[result_gdf.intersects(mask.geometry)] + # Prioritize downstream geometries within the intersecting set + mask_downstream = intersecting_geometries[intersecting_geometries["geometry"].isin(downstream_gdf["geometry"])] + if not mask_downstream.empty: + mask_gdf = mask_downstream + else: + mask_gdf = intersecting_geometries + # Add the mask geometry as a new column for reference + mask_gdf = mask_gdf.copy() + mask_gdf["geometry_mask"] = mask.geometry + + list_result_gdfs.append(mask_gdf) + + return list_result_gdfs diff --git a/lidro/create_virtual_point/vectors/run_create_virtual_points.py b/lidro/create_virtual_point/vectors/run_create_virtual_points.py index 487848d2..06ab77d6 100644 --- a/lidro/create_virtual_point/vectors/run_create_virtual_points.py +++ b/lidro/create_virtual_point/vectors/run_create_virtual_points.py @@ -1,105 +1,67 @@ # -*- coding: utf-8 -*- """ Run function "virtual points" """ -import logging -import os - import geopandas as gpd import numpy as np -import pandas as pd +from scipy.interpolate import interp1d +from shapely.geometry import Point -from lidro.create_virtual_point.vectors.apply_Z_from_grid import ( - calculate_grid_z, - calculate_grid_z_with_model, -) from lidro.create_virtual_point.vectors.create_grid_2D_inside_maskhydro import ( generate_grid_from_geojson, ) -from lidro.create_virtual_point.vectors.flatten_river import flatten_little_river -from lidro.create_virtual_point.vectors.linear_regression_model import ( - calculate_linear_regression_line, -) def compute_virtual_points_by_section( - points: gpd.GeoDataFrame, - line: gpd.GeoDataFrame, + line_gdf: gpd.GeoDataFrame, mask_hydro: gpd.GeoDataFrame, crs: str, spacing: float, - length: int, - output_dir: str, + # output_dir: str, ) -> gpd.GeoDataFrame: """This function generates a regular grid of 3D points spaced every N meters inside each hydro entity = virtual point Args: - points (gpd.GeoDataFrame): A GeoDataFrame containing points along Hydro's Skeleton - line (gpd.GeoDataFrame): A GeoDataFrame containing each line from Hydro's Skeleton + line_gdf (gpd.GeoDataFrame): A GeoDataFrame containing each line from Hydro's Skeleton 3D (rectify or not) mask_hydro (gpd.GeoDataFrame): A GeoDataFrame containing each mask hydro from Hydro's Skeleton crs (str): A pyproj CRS object used to create the output GeoJSON file spacing (float, optional): Spacing between the grid points in meters. The default value is 0.5 meter - length (int, optional): Minimum length of a river to use the linear regression model. - The default value is 150 meters. - output_dir (str): Folder to output Mask Hydro without virtual points Returns: gpd.GeoDataFrame: All virtual points by Mask Hydro """ - # Check if the points DataFrame is empty and all the values in the "points_knn" column are null - if points.empty or points["points_knn"].isnull().all(): - logging.warning("The points GeoDataFrame is empty. Saving the skeleton and mask hydro to GeoJSON.") - masks_without_points = gpd.GeoDataFrame(columns=mask_hydro.columns, crs=mask_hydro.crs) - for idx, mask in mask_hydro.iterrows(): - logging.warning(f"No points found within mask hydro {idx}. Adding to masks_without_points.") - masks_without_points = pd.concat([masks_without_points, gpd.GeoDataFrame([mask], crs=mask_hydro.crs)]) - # Save the resulting masks_without_points to a GeoJSON file - output_mask_hydro_error = os.path.join(output_dir, "mask_hydro_no_virtual_points.geojson") - masks_without_points.to_file(output_mask_hydro_error, driver="GeoJSON") - else: - # Step 1: Generates a regular 2D grid of evenly spaced points within a Mask Hydro - gdf_grid = generate_grid_from_geojson(mask_hydro, spacing) - # Calculate the length of the river - river_length = float(line.length.iloc[0]) + # Generate grid 2D inside mask hydro + grid_gdf = generate_grid_from_geojson(mask_hydro, spacing) + + # Extract the polyline from skeleton hydro (assuming there is only one geometry in line_gdf) + skeleton_line = line_gdf.iloc[0].geometry + # Extract points and their coordinates + polyline_coords = [Point(coord) for coord in skeleton_line.coords] + # Extract Z values from the 3D coordinates + z_values = np.array([coord[2] for coord in skeleton_line.coords]) # assuming the polyline is in 3D + + # Compute cumulative distances along the polyline for interpolation + cumulative_distances = [0] + [ + skeleton_line.project(Point(polyline_coords[i].x, polyline_coords[i].y)) + for i in range(1, len(polyline_coords)) + ] + + # Interpolation function for Z based on cumulative distances along the polyline + z_interpolator = interp1d(cumulative_distances, z_values, kind="linear", fill_value="extrapolate") + + # Assign Z values to grid points by projecting each point onto the polyline + z_values_for_grid = [] + for point in grid_gdf.geometry: + # Project the 2D point onto the polyline to get the distance along the polyline + distance_on_polyline = skeleton_line.project(point) + # Interpolate the Z value at this distance + z_value = z_interpolator(distance_on_polyline) + z_values_for_grid.append(z_value) - # Step 2 : Model linear regression for river's length > 150 m - if river_length > length: - model, r2 = calculate_linear_regression_line(points, line, crs) - if model == np.poly1d([0, 0]) and r2 == 0.0: - masks_without_points = gpd.GeoDataFrame(columns=mask_hydro.columns, crs=mask_hydro.crs) - for idx, mask in mask_hydro.iterrows(): - masks_without_points = pd.concat( - [masks_without_points, gpd.GeoDataFrame([mask], crs=mask_hydro.crs)] - ) - # Save the resulting masks_without_points because of linear regression is impossible to a GeoJSON file - output_mask_hydro_error = os.path.join( - output_dir, "mask_hydro_no_virtual_points_with_regression.geojson" - ) - logging.warning( - f"Save masks_without_points because linear regression is impossible: {output_mask_hydro_error}" - ) - masks_without_points.to_file(output_mask_hydro_error, driver="GeoJSON") - # Apply linear regression model at the rivers - gdf_grid_with_z = calculate_grid_z_with_model(gdf_grid, line, model) - # Step 2 bis: Model flattening for river's length < 150 m or river's length == 150 m - else: - predicted_z = flatten_little_river(points, line) - if predicted_z == 0: - masks_without_points = gpd.GeoDataFrame(columns=mask_hydro.columns, crs=mask_hydro.crs) - for idx, mask in mask_hydro.iterrows(): - masks_without_points = pd.concat( - [masks_without_points, gpd.GeoDataFrame([mask], crs=mask_hydro.crs)] - ) - # Save the resulting masks_without_points because of flattening river is impossible to a GeoJSON file - output_mask_hydro_error = os.path.join( - output_dir, "mask_hydro_no_virtual_points_for_little_river.geojson" - ) - logging.warning( - f"Save masks_without_points because of flattening river is impossible: {output_mask_hydro_error}" - ) - masks_without_points.to_file(output_mask_hydro_error, driver="GeoJSON") - # Apply flattening model at the rivers - gdf_grid_with_z = calculate_grid_z(gdf_grid, predicted_z) + # Create a new GeoDataFrame with X, Y, and interpolated Z values as 3D points + grid_with_z = gpd.GeoDataFrame( + geometry=gpd.points_from_xy(grid_gdf.geometry.x, grid_gdf.geometry.y, z_values_for_grid), crs=grid_gdf.crs + ) - return gdf_grid_with_z + return grid_with_z diff --git a/lidro/create_virtual_point/vectors/run_update_skeleton_with_z.py b/lidro/create_virtual_point/vectors/run_update_skeleton_with_z.py new file mode 100644 index 00000000..2ba3809d --- /dev/null +++ b/lidro/create_virtual_point/vectors/run_update_skeleton_with_z.py @@ -0,0 +1,96 @@ +# -*- coding: utf-8 -*- +""" Run function "update skeleton with z" +""" +import logging +import os +from typing import List + +import geopandas as gpd +import numpy as np +import pandas as pd + +from lidro.create_virtual_point.vectors.create_skeleton_3d import ( + skeleton_3d_with_flatten, + skeleton_3d_with_linear_regression, +) +from lidro.create_virtual_point.vectors.flatten_river import flatten_little_river +from lidro.create_virtual_point.vectors.linear_regression_model import ( + calculate_linear_regression_line, +) + + +def compute_skeleton_with_z( + points: gpd.GeoDataFrame, + line: gpd.GeoDataFrame, + mask_hydro: gpd.GeoDataFrame, + crs: str, + length: int, + output_dir: str, +) -> List: + """This function update skeleton with Z + + Args: + points (gpd.GeoDataFrame): A GeoDataFrame containing points along Hydro's Skeleton + line (gpd.GeoDataFrame): A GeoDataFrame containing each line from Hydro's Skeleton + mask_hydro (gpd.GeoDataFrame): A GeoDataFrame containing each mask hydro from Hydro's Skeleton + length (int, optional): Minimum length of a river to use the linear regression model. + The default value is 150 meters. + output_dir (str): Folder to output Mask Hydro without virtual points + + Returns: + List : the 3D skeleton's geometrie (Polyligne Z) Z with Z mean + """ + # Check if the points DataFrame is empty and all the values in the "points_knn" column are null + if points.empty or points["points_knn"].isnull().all(): + logging.warning("The points GeoDataFrame is empty. Saving the skeleton and mask hydro to GeoJSON.") + masks_without_points = gpd.GeoDataFrame(columns=mask_hydro.columns, crs=mask_hydro.crs) + for idx, mask in mask_hydro.iterrows(): + logging.warning(f"No points found within mask hydro {idx}. Adding to masks_without_points.") + masks_without_points = pd.concat([masks_without_points, gpd.GeoDataFrame([mask], crs=mask_hydro.crs)]) + # Save the resulting masks_without_points to a GeoJSON file + output_mask_hydro_error = os.path.join(output_dir, "mask_hydro_no_virtual_points.geojson") + masks_without_points.to_file(output_mask_hydro_error, driver="GeoJSON") + else: + # Calculate the length of the river + river_length = float(line.length.iloc[0]) + + # Step 2 : Model linear regression for river's length > 150 m + if river_length > length: + model, r2 = calculate_linear_regression_line(points, line, crs) + if model == np.poly1d([0, 0]) and r2 == 0.0: + masks_without_points = gpd.GeoDataFrame(columns=mask_hydro.columns, crs=mask_hydro.crs) + for idx, mask in mask_hydro.iterrows(): + masks_without_points = pd.concat( + [masks_without_points, gpd.GeoDataFrame([mask], crs=mask_hydro.crs)] + ) + # Save the resulting masks_without_points because of linear regression is impossible to a GeoJSON file + output_mask_hydro_error = os.path.join( + output_dir, "mask_hydro_no_virtual_points_with_regression.geojson" + ) + logging.warning( + f"Save masks_without_points because linear regression is impossible: {output_mask_hydro_error}" + ) + masks_without_points.to_file(output_mask_hydro_error, driver="GeoJSON") + # Apply model from skeleton + skeleton_hydro_3D = skeleton_3d_with_linear_regression(line, model, crs) + # Step 2 bis: Model flattening for river's length < 150 m or river's length == 150 m + else: + predicted_z = flatten_little_river(points, line) + if predicted_z == 0: + masks_without_points = gpd.GeoDataFrame(columns=mask_hydro.columns, crs=mask_hydro.crs) + for idx, mask in mask_hydro.iterrows(): + masks_without_points = pd.concat( + [masks_without_points, gpd.GeoDataFrame([mask], crs=mask_hydro.crs)] + ) + # Save the resulting masks_without_points because of flattening river is impossible to a GeoJSON file + output_mask_hydro_error = os.path.join( + output_dir, "mask_hydro_no_virtual_points_for_little_river.geojson" + ) + logging.warning( + f"Save masks_without_points because of flattening river is impossible: {output_mask_hydro_error}" + ) + masks_without_points.to_file(output_mask_hydro_error, driver="GeoJSON") + # Apply model from skeleton + skeleton_hydro_3D = skeleton_3d_with_flatten(line, predicted_z, crs) + + return skeleton_hydro_3D diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index b17861fc..d8384e27 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -17,9 +17,15 @@ from lidro.create_virtual_point.vectors.merge_skeleton_by_mask import ( merge_skeleton_by_mask, ) +from lidro.create_virtual_point.vectors.rectify_alert import ( + create_list_rectify_skeleton_with_mask, +) from lidro.create_virtual_point.vectors.run_create_virtual_points import ( compute_virtual_points_by_section, ) +from lidro.create_virtual_point.vectors.run_update_skeleton_with_z import ( + compute_skeleton_with_z, +) @hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") @@ -52,6 +58,7 @@ def main(config: DictConfig): input_mask_hydro = config.io.input_mask_hydro input_skeleton = config.io.input_skeleton input_dir_points_skeleton = config.io.dir_points_skeleton + input_bridge = config.io.input_bridge crs = CRS.from_user_input(config.io.srid) river_length = config.virtual_point.vector.river_length points_grid_spacing = config.virtual_point.pointcloud.points_grid_spacing @@ -80,22 +87,42 @@ def process_points_knn(points_knn): # Combine skeleton lines into a single polyline for each hydro entity gdf_merged = merge_skeleton_by_mask(input_skeleton, input_mask_hydro, output_dir, crs) - # Step 3 : Generate a regular grid of 3D points spaced every N meters inside each hydro entity - list_virtual_points = [ - compute_virtual_points_by_section( + # Step 3 : Apply Z from skeleton + list_skeleton_with_z = [ + compute_skeleton_with_z( points_gdf, gpd.GeoDataFrame([{"geometry": row["geometry_skeleton"]}], crs=crs), gpd.GeoDataFrame([{"geometry": row["geometry_mask"]}], crs=crs), crs, - points_grid_spacing, river_length, output_dir, ) for idx, row in gdf_merged.iterrows() ] + logging.info("Apply Z to skeleton") + # print(list_skeleton_with_z) + + # Step 4 : intersect skeletons by bridge and return info + gdf_skeleton_rectify_with_mask = create_list_rectify_skeleton_with_mask( + input_bridge, input_mask_hydro, list_skeleton_with_z, crs + ) + logging.info("intersect skeletons by bridge and return info") + print(gdf_skeleton_rectify_with_mask) + + # Step 5 : Generate a regular grid of 3D points spaced every N meters inside each hydro entity + list_virtual_points = [ + compute_virtual_points_by_section( + gpd.GeoDataFrame(geometry=row["geometry"], crs=crs), + gpd.GeoDataFrame(geometry=row["geometry_mask"], crs=crs), + crs, + points_grid_spacing, + ) + for idx, row in enumerate(gdf_skeleton_rectify_with_mask) + ] logging.info("Calculate virtuals points by mask hydro and skeleton") + # print(list_virtual_points) - # Step 4 : Save the virtual points in a file (.LAZ) + # Step 6 : Save the virtual points in a file (.LAZ) list_points_to_las(list_virtual_points, output_dir, crs, classes) else: logging.error("Error when merged all points around skeleton by lidar tile")