diff --git a/README.md b/README.md index 4df90f02..9b281420 100644 --- a/README.md +++ b/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 diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index 4f98e4b3..22ecdc28 100644 --- a/configs/configs_lidro.yaml +++ b/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: diff --git a/data b/data index 19f969cd..36a33e74 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 19f969cd0f938dcd5ca3d7019578fa74785c9934 +Subproject commit 36a33e74ac841a0c59a0f1692781fdec2eaee8eb diff --git a/images/chaine_traitement_lidro.jpg b/images/chaine_traitement_lidro.jpg new file mode 100644 index 00000000..e48c8dbe Binary files /dev/null and b/images/chaine_traitement_lidro.jpg differ diff --git a/images/process_points_virtuels.jpg b/images/process_points_virtuels.jpg new file mode 100644 index 00000000..123ca1d5 Binary files /dev/null and b/images/process_points_virtuels.jpg differ diff --git a/lidro/create_virtual_point/pointcloud/convert_list_points_to_las.py b/lidro/create_virtual_point/pointcloud/convert_list_points_to_las.py new file mode 100644 index 00000000..b647b3a5 --- /dev/null +++ b/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}") diff --git a/lidro/create_virtual_point/stats/knn_distance.py b/lidro/create_virtual_point/stats/knn_distance.py index 7362bcbe..3b47bac8 100644 --- a/lidro/create_virtual_point/stats/knn_distance.py +++ b/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 diff --git a/lidro/create_virtual_point/vectors/apply_Z_from_grid.py b/lidro/create_virtual_point/vectors/apply_Z_from_grid.py new file mode 100644 index 00000000..9e6d497b --- /dev/null +++ b/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 diff --git a/lidro/create_virtual_point/vectors/create_grid_2D_inside_maskhydro.py b/lidro/create_virtual_point/vectors/create_grid_2D_inside_maskhydro.py new file mode 100644 index 00000000..4f9b5f11 --- /dev/null +++ b/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 diff --git a/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py b/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py index a7a5cbda..7bd4b11b 100644 --- a/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py +++ b/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py @@ -1,11 +1,13 @@ # -*- coding: utf-8 -*- """ Extract points around skeleton by tile """ +import json import logging import os from typing import List import geopandas as gpd +import pandas as pd from pdaltools.las_info import las_get_xy_bounds from shapely.geometry import Point @@ -21,8 +23,10 @@ def extract_points_around_skeleton_points_one_tile( filename: str, input_dir: str, + output_dir: str, input_mask_hydro_buffer: gpd.GeoDataFrame, points_skeleton_gdf: gpd.GeoDataFrame, + crs: str | int, classes: List[int:int], k: int, ): @@ -33,13 +37,12 @@ def extract_points_around_skeleton_points_one_tile( Args: filename (str): filename to the LAS file input_dir (str): input folder + output_dir (str): ouput folder input_mask_hydro_buffer (gpd.GeoDataFrame): hydro mask with buffer - points_skeleton_gdf (gpd.GeoDataFrame): Points every 2 meters (by default) along skeleton hydro + points_skeleton_gdf (gpd.GeoDataFrame): Points every N meters along skeleton hydro + crs (str | int): Make a CRS from an EPSG code : CRS WKT string classes (List): List of classes to use for the filtering k (int): the number of nearest neighbors to find - - Returns: - points_clip (np.array) : Numpy array containing point coordinates (X, Y, Z) after filtering and croping """ # Step 1 : Crop filtered pointcloud by Mask Hydro with buffer input_dir_points = os.path.join(input_dir, "pointcloud") @@ -58,6 +61,27 @@ def extract_points_around_skeleton_points_one_tile( ] # Step 2 : Extract a Z elevation value along the hydrographic skeleton logging.info(f"\nExtract a Z elevation value along the hydrographic skeleton for tile : {tilename}") - result = filter_las_around_point(points_skeleton_with_z_clip, points_clip, k) + points_Z = filter_las_around_point(points_skeleton_with_z_clip, points_clip, k) + + if len(points_Z) > 0: + # Limit the precision of coordinates using numpy arrays or tuples + knn_points = [ + [[round(coord, 3) for coord in point] for point in item["points_knn"]] + for item in points_Z + if len(item["points_knn"]) > 0 and isinstance(item["geometry"], Point) + ] + # Create a DataFrame with lists or numpy arrays + df_points_z = pd.DataFrame({"geometry": [item["geometry"] for item in points_Z], "points_knn": knn_points}) + + # Encode the 'points_knn' lists as JSON + df_points_z["points_knn"] = df_points_z["points_knn"].apply(lambda x: json.dumps(x)) + + # Convert DataFrame to GeoDataFrame + gdf_points_z = gpd.GeoDataFrame(df_points_z, geometry="geometry") + gdf_points_z.set_crs(crs, inplace=True) + + # Save the GeoDataFrame to a GeoJSON file + output_geojson_path = os.path.join(output_dir, "_points_skeleton".join([tilename, ".geojson"])) + gdf_points_z.to_file(output_geojson_path, driver="GeoJSON") - return result + logging.info(f"Result saved to {output_geojson_path}") diff --git a/lidro/create_virtual_point/vectors/flatten_river.py b/lidro/create_virtual_point/vectors/flatten_river.py new file mode 100644 index 00000000..1792e399 --- /dev/null +++ b/lidro/create_virtual_point/vectors/flatten_river.py @@ -0,0 +1,55 @@ +# -*- coding: utf-8 -*- +""" this function flattens the Z of the river +""" +import logging + +import geopandas as gpd +import numpy as np + +from lidro.create_virtual_point.vectors.intersect_points_by_line import ( + return_points_by_line, +) + + +def flatten_little_river(points: gpd.GeoDataFrame, line: gpd.GeoDataFrame): + """This function flattens the Z of the little rivers, because of + regression model doesn't work for these type of rivers. + + Args: + points (gpd.GeoDataFrame): A GeoDataFrame containing points along Hydro's Skeleton + line (gpd.GeoDataFrame): A GeoDataFrame containing each line from Hydro's Skeleton + + Returns: + float: Z of the river + """ + # Inputs + gdf_points = return_points_by_line(points, line) + # Extract points_knn and drop NaNs + points_knn_values = gdf_points["points_knn"].values + + # Check if points_knn_values is empty + if points_knn_values.size == 0: + logging.warning("For little river: no neighbor found to calculate Z quartile: set Z to 0") + z_first_quartile = 0 + return z_first_quartile + + # Merge points and remove duplicates + try: + all_points_knn = np.vstack(points_knn_values) + except ValueError as e: + logging.error(f"Error during stacking points: {e}") + z_first_quartile = 0 + return z_first_quartile + + unique_points_knn = np.unique(all_points_knn, axis=0) + + # Check if unique_points_knn is empty before calculating the first quartile + if unique_points_knn.size == 0: + logging.warning("No unique points found after filtering. Setting Z to 0.") + z_first_quartile = 0 + else: + # Calculate the 1st quartile of all points + first_quartile = np.percentile(unique_points_knn, 25, axis=0) + z_first_quartile = first_quartile[-1] + + return z_first_quartile diff --git a/lidro/create_virtual_point/vectors/intersect_points_by_line.py b/lidro/create_virtual_point/vectors/intersect_points_by_line.py new file mode 100644 index 00000000..0a94592c --- /dev/null +++ b/lidro/create_virtual_point/vectors/intersect_points_by_line.py @@ -0,0 +1,28 @@ +# -*- coding: utf-8 -*- +""" Identifies all points that intersect each line +""" +import geopandas as gpd +from shapely.geometry import CAP_STYLE + + +def return_points_by_line(points: gpd.GeoDataFrame, line: gpd.GeoDataFrame) -> gpd.GeoDataFrame: + """This function identifies all points that intersect each Skeleton by hydro's section + + Args: + points (gpd.GeoDataFrame): A GeoDataFrame containing points along Hydro's Skeleton + line (gpd.GeoDataFrame): A GeoDataFrame containing each line from Hydro's Skeleton + + Returns: + gpd.GeoDataFrame: A GeoDataframe containing only points that intersect each hydro's section + """ + # Apply buffer (0.1 meters) from Mask Hydro + line_buffer = line.buffer(0.1, cap_style=CAP_STYLE.square) + gdf_line_buffer = gpd.GeoDataFrame(geometry=line_buffer) + + # Perform spatial join to find intersections + pts_intersect = gpd.sjoin(points, gdf_line_buffer, how="left", predicate="intersects") + + # Drop lines which contain one or more NaN values + pts_intersect = pts_intersect.dropna() + + return pts_intersect diff --git a/lidro/create_virtual_point/vectors/las_around_point.py b/lidro/create_virtual_point/vectors/las_around_point.py index ece0dec5..13cbd9e8 100644 --- a/lidro/create_virtual_point/vectors/las_around_point.py +++ b/lidro/create_virtual_point/vectors/las_around_point.py @@ -1,17 +1,16 @@ # -*- coding: utf-8 -*- -""" Extract a Z elevation value every N meters along the hydrographic skeleton +""" Extract a list of N "ground" points closest to the N points of the hydro skeleton """ from typing import List import numpy as np from shapely.geometry import Point -from lidro.create_virtual_point.stats.calculate_stat import calculate_percentile from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors def filter_las_around_point(points_skeleton: List, points_clip: np.array, k: int) -> List: - """Extract a Z elevation value every N meters along the hydrographic skeleton + """Extract a list of N "ground" points closest to the N points of the hydro skeleton Args: points_skeleton (list) : points every N meters (by default: 2) along skeleton Hydro @@ -19,7 +18,8 @@ def filter_las_around_point(points_skeleton: List, points_clip: np.array, k: int k (int): The number of nearest neighbors to find Returns: - List : Result {'geometry': Point 3D, 'z_q1': points KNN} + List : Result {'geometry': Point 2D on skeleton, + 'points_knn': List of N "ground" points closest to 2D points of the hydro skeleton} """ # Finds the K nearest neighbors of a given point from a list of 3D points points_knn_list = [ @@ -27,11 +27,4 @@ def filter_las_around_point(points_skeleton: List, points_clip: np.array, k: int for point in points_skeleton ] - # Calcule Z "Q1" for each points every N meters along skeleton hydro - results = [ - ({"geometry": p["geometry"], "z_q1": calculate_percentile(p["points_knn"], 10)}) - for p in points_knn_list - if not np.all(np.isnan(p["points_knn"])) - ] - - return results + return points_knn_list diff --git a/lidro/create_virtual_point/vectors/linear_regression_model.py b/lidro/create_virtual_point/vectors/linear_regression_model.py new file mode 100644 index 00000000..f2a7879f --- /dev/null +++ b/lidro/create_virtual_point/vectors/linear_regression_model.py @@ -0,0 +1,96 @@ +# -*- coding: utf-8 -*- +""" This function calculates a linear regression line +in order to read the Zs along the hydro skeleton to guarantee the flow +""" +import logging + +import geopandas as gpd +import numpy as np +from shapely import line_locate_point + +from lidro.create_virtual_point.vectors.intersect_points_by_line import ( + return_points_by_line, +) + + +def calculate_linear_regression_line(points: gpd.GeoDataFrame, line: gpd.GeoDataFrame, crs: str): + """This function calculates a linear regression line in order to read the Zs + along the hydro skeleton to guarantee the flow. + A river is a natural watercourse whose slope is less than 1%, and + the error tolerance specified in the specifications is 0.5 m. + So, the linear regression model must have a step of 50 m. + + Args: + points (gpd.GeoDataFrame): A GeoDataFrame containing points along Hydro's Skeleton + line (gpd.GeoDataFrame): A GeoDataFrame containing each line from Hydro's Skeleton + crs (str): a pyproj CRS object used to create the output GeoJSON file + + Returns: + np.poly1d: Regression model + numpy.array: Determination coefficient + """ + if points.empty or line.empty: + logging.warning("Input GeoDataFrames 'points' and 'line' must not be empty") + return np.poly1d([0, 0]), 0.0 + + # Retrieve points along the line + gdf_points = return_points_by_line(points, line) + + if gdf_points.empty or len(gdf_points) < 3: + logging.warning("Not enough points for regression analysis") + return np.poly1d([0, 0]), 0.0 + + # Combine all `points_knn` into a single list + all_points_knn = [] + for knn_list in gdf_points["points_knn"]: + all_points_knn.extend(knn_list) + + # Convert the list of points to a numpy array and remove duplicates + all_points_knn_array = np.array(all_points_knn) + unique_points_knn = np.unique(all_points_knn_array, axis=0) + + # Create final data structure + final_data = {"geometry": [line.iloc[0]["geometry"]], "points_knn": unique_points_knn} + + # Generate projected coordinates + points_gs = gpd.GeoSeries().from_xy( + final_data["points_knn"][:, 0], # X coordinates + final_data["points_knn"][:, 1], # Y coordinates + final_data["points_knn"][:, 2], # Z coordinates + crs=crs, # Coordinate reference system + ) + + # Locate each point along the polyline + d_line = line_locate_point(final_data["geometry"], points_gs, normalized=False) + + # Create a GeoDataFrame for regression analysis + regression_gpd = gpd.GeoDataFrame(geometry=points_gs) + # This column contains the curvilinear abscissa + regression_gpd["ac"] = d_line + bins = np.arange(regression_gpd["ac"].max(), step=50) # Create bins for curvilinear abscissa + regression_gpd["ac_bin"] = np.digitize(regression_gpd["ac"], bins) # Digitize curvilinear abscissa into bins + # This column contains the Z value of the point + regression_gpd["z"] = regression_gpd.geometry.z + + # Linear regression model using binned data + temp = regression_gpd.groupby("ac_bin").aggregate( + { + "ac": "mean", # Mean of curvilinear abscissa + "z": [ + ("quantile", lambda x: x.quantile(0.1)), # 10th percentile of Z values + ("std", lambda x: x.quantile(0.75) - x.quantile(0.25)), + ], # Interquartile range of Z values + } + ) + # Linear regression with weights + coeff, sse, *_ = np.polyfit(temp["ac"]["mean"], temp["z"]["quantile"], deg=1, full=True) + + # Calculate SST (TOTAL SQUARE SUM) + sst = np.sum((temp["z"]["quantile"] - np.mean(temp["z"]["quantile"])) ** 2) + + # Determination coefficient [0, 1] + r2 = 1 - (sse / sst) + + model = np.poly1d(coeff) + + return model, r2 diff --git a/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py b/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py index c8ab157b..e0530059 100644 --- a/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py +++ b/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py @@ -1,12 +1,12 @@ # -*- coding: utf-8 -*- -""" Extract a Z elevation value every N meters along the hydrographic skeleton +""" Create a new Mask Hydro at the edge of the bank with a buffer """ import geopandas as gpd from shapely.geometry import CAP_STYLE def import_mask_hydro_with_buffer(file_mask: str, buffer: float, crs: str | int) -> gpd.GeoDataFrame: - """Apply buffer (2 meters by default) from Mask Hydro + """Create a new Mask Hydro at the edge of the bank with a buffer (+50 cm and -N cm) Args: file_mask (str): Path from Mask Hydro @@ -19,10 +19,8 @@ def import_mask_hydro_with_buffer(file_mask: str, buffer: float, crs: str | int) # Read Mask Hydro merged gdf = gpd.read_file(file_mask, crs=crs).unary_union - # Apply buffer (2 meters by default) from Mask Hydro - gdf_buffer = gdf.buffer(buffer, cap_style=CAP_STYLE.square) - + # Apply negative's buffer + difference from Mask Hydro # Return a polygon representing the limit of the bank with a buffer of N meters - limit_bank = gdf_buffer.difference(gdf) + gdf_buffer = gdf.difference(gdf.buffer(-abs(buffer), cap_style=CAP_STYLE.square)) - return limit_bank + return gdf_buffer diff --git a/lidro/create_virtual_point/vectors/merge_skeleton_by_mask.py b/lidro/create_virtual_point/vectors/merge_skeleton_by_mask.py new file mode 100644 index 00000000..fd7f770b --- /dev/null +++ b/lidro/create_virtual_point/vectors/merge_skeleton_by_mask.py @@ -0,0 +1,114 @@ +# -*- coding: utf-8 -*- +""" Merge skeleton by Mask Hydro +""" +import logging +import os + +import geopandas as gpd +import pandas as pd +from shapely.geometry import LineString, MultiLineString + + +def explode_multipart(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame: + """Explode multi-part geometries into single-part geometries + + Args: + gdf (gpd.GeoDataFrame): A GeoDataFrame + + Returns: + gpd.GeoDataframe : a GeoDataFrame exploded + """ + exploded_geom = gdf.explode(index_parts=True) + exploded_geom = exploded_geom.reset_index(drop=True) + + return exploded_geom + + +def combine_and_connect_lines(geometries: list) -> LineString: + """Combines and connects a list of LineString and MultiLineString geometries into a single LineString. + + Args: + geometries (list): A list of LineString and/or MultiLineString objects. + + Returns: + LineString: The merged and connected geometry. + """ + # Convert all geometries into lists of lines + lines = [geom.geoms if isinstance(geom, MultiLineString) else [geom] for geom in geometries] + # Extract and combine all coordinates from the nested list using a list comprehension + all_coords = [ + coord for sublist in lines for line in sublist if isinstance(line, LineString) for coord in line.coords + ] + # Create a new connected line + connected_line = LineString(all_coords) + + return connected_line + + +def merge_skeleton_by_mask( + input_skeleton: gpd.GeoDataFrame, input_mask_hydro: gpd.GeoDataFrame, output_dir: str, crs: str +) -> pd.DataFrame: + """Combine skeleton lines into a single polyline for each hydro entity. + Save a invalid mask to GeoDataframe + + Args: + input_skeleton (gpd.GeoDataFrame): A GeoDataFrame containing each line from Hydro's Skeleton + input_mask_hydro (gpd.GeoDataFrame): A GeoDataFrame containing Mask Hydro + output_dir (str): output folder + crs (str): a pyproj CRS object used to create the out + + Returns: + Dataframe : Dataframe with single polyline "skeleton" by hydro entity + """ # Inputs + gdf_skeleton = gpd.read_file(input_skeleton) + gdf_skeleton = explode_multipart(gdf_skeleton) + gdf_mask_hydro = gpd.read_file(input_mask_hydro) + gdf_mask_hydro = explode_multipart(gdf_mask_hydro) + + # # Perform a spatial join to find skeletons within each mask_hydro + gdf_joined = gpd.sjoin(gdf_skeleton, gdf_mask_hydro, how="inner", predicate="intersects") + # geometry intersections with funtion "overlay" + gdf_intersections = gpd.overlay(gdf_joined, gdf_mask_hydro, how="intersection", keep_geom_type=True) + + # Combine skeleton lines into a single polyline for each hydro entity + combined_skeletons = gdf_intersections.groupby("index_right")["geometry"].apply(combine_and_connect_lines) + gdf_combined_skeletons = gpd.GeoDataFrame(combined_skeletons, columns=["geometry"], crs=crs).reset_index() + # # Re-join with mask_hydro to keep only the combined skeletons within masks + gdf_joined_combined = gpd.sjoin( + gdf_combined_skeletons, gdf_mask_hydro, how="inner", predicate="intersects", lsuffix="combined", rsuffix="mask" + ) + + # # Count the number of skeletons per mask + skeleton_counts = gdf_joined_combined["index_mask"].value_counts() + + valid_masks = skeleton_counts[skeleton_counts == 1].index + invalid_masks = skeleton_counts[skeleton_counts != 1].index + + # Filter to keep only masks with a single skeleton + if not valid_masks.empty: + # Filter the joined GeoDataFrame to keep only valid masks + gdf_valid_joined = gdf_joined_combined[gdf_joined_combined["index_mask"].isin(valid_masks)] + # Merge the geometries of masks and skeletons into the resulting GeoDataFrame + df_result = gdf_valid_joined.merge( + gdf_mask_hydro, left_on="index_mask", right_index=True, suffixes=("_skeleton", "_mask") + ) + # Keep only necessary columns + df_result = df_result[["geometry_skeleton", "geometry_mask"]] + + if not invalid_masks.empty: + # Filter the joined GeoDataFrame to keep invalid masks + gdf_invalid_joined = gdf_joined_combined[gdf_joined_combined["index_mask"].isin(valid_masks)] + # Merge the geometries of masks and skeletons into the resulting GeoDataFrames + df_exclusion = gdf_invalid_joined.merge( + gdf_mask_hydro, left_on="index_mask", right_index=True, suffixes=("_skeleton", "_mask") + ) + # Keep only necessary columns + df_exclusion = df_exclusion[["geometry_mask"]] + # Save the results and exclusions to separate GeoJSON files + # Convert DataFrame to GeoDataFrame + gdf_exclusion = gpd.GeoDataFrame(df_exclusion, geometry="geometry_mask") + exclusions_outputs = os.path.join(output_dir, "mask_skeletons_exclusions.geojson") + gdf_exclusion.to_file(exclusions_outputs, driver="GeoJSON") + logging.info(f"Error: Save the exclusions rivers in GeoJSON : {exclusions_outputs}") + + return df_result diff --git a/lidro/create_virtual_point/vectors/run_create_virtual_points.py b/lidro/create_virtual_point/vectors/run_create_virtual_points.py new file mode 100644 index 00000000..d30ad444 --- /dev/null +++ b/lidro/create_virtual_point/vectors/run_create_virtual_points.py @@ -0,0 +1,105 @@ +# -*- coding: utf-8 -*- +""" Run function "virtual points" +""" +import logging +import os + +import geopandas as gpd +import numpy as np +import pandas as pd + +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 launch_virtual_points_by_section( + points: gpd.GeoDataFrame, + line: gpd.GeoDataFrame, + mask_hydro: gpd.GeoDataFrame, + crs: str, + spacing: float, + length: int, + 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 + 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]) + + # 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) + + return gdf_grid_with_z diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 2ead2845..3e9ef544 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -1,10 +1,10 @@ -""" Main script for calculate Mask HYDRO 1 +""" Main script for create virtuals points """ +import ast import logging import os import sys - import geopandas as gpd import hydra import pandas as pd @@ -13,20 +13,20 @@ sys.path.append('../lidro') -from lidro.create_virtual_point.vectors.extract_points_around_skeleton import ( - extract_points_around_skeleton_points_one_tile, +from lidro.create_virtual_point.pointcloud.convert_list_points_to_las import ( + list_points_to_las, ) -from lidro.create_virtual_point.vectors.mask_hydro_with_buffer import ( - import_mask_hydro_with_buffer, +from lidro.create_virtual_point.vectors.merge_skeleton_by_mask import ( + merge_skeleton_by_mask, ) -from lidro.create_virtual_point.vectors.points_along_skeleton import ( - generate_points_along_skeleton, +from lidro.create_virtual_point.vectors.run_create_virtual_points import ( + launch_virtual_points_by_section, ) @hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") def main(config: DictConfig): - """Create a virtual point along hydro surfaces from the points classification of + """Create a virtual point inside hydro surfaces (3D grid) from the points classification of the input LAS/LAZ file and the Hyro Skeleton (GeoJSON) and save it as LAS file. It can run either on a single file, or on each file of a folder @@ -36,7 +36,6 @@ def main(config: DictConfig): 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: @@ -51,55 +50,56 @@ def main(config: DictConfig): os.makedirs(output_dir, exist_ok=True) - # If input filename is not provided, lidro runs on the whole input_dir directory - initial_las_filename = config.io.input_filename - + # Parameters for creating virtual point input_mask_hydro = config.io.input_mask_hydro input_skeleton = config.io.input_skeleton - - # Parameters for creating virtual point - distance_meters = config.virtual_point.vector.distance_meters - buffer = config.virtual_point.vector.buffer + input_dir_points_skeleton = config.io.dir_points_skeleton crs = CRS.from_user_input(config.io.srid) - classes = config.virtual_point.filter.keep_neighbors_classes - k = config.virtual_point.vector.k - - # Step 1 : Import Mask Hydro, then apply a buffer - # Return GeoDataframe - input_mask_hydro_buffer = import_mask_hydro_with_buffer(input_mask_hydro, buffer, crs).wkt - - # Step 2 : Create several points every 2 meters (by default) along skeleton Hydro - # Return GeoDataframe - points_skeleton_gdf = generate_points_along_skeleton(input_skeleton, distance_meters, crs) - - # Step 3 : Extract points around skeleton by tile - if initial_las_filename: - # Lauch croping filtered pointcloud by Mask Hydro with buffer by one tile: - points_clip = extract_points_around_skeleton_points_one_tile( - initial_las_filename, input_dir, input_mask_hydro_buffer, points_skeleton_gdf, classes, k - ) - - else: - # Lauch croping filtered pointcloud by Mask Hydro with buffer tile by tile - input_dir_points = os.path.join(input_dir, "pointcloud") - points_clip_list = [ - extract_points_around_skeleton_points_one_tile( - file, input_dir, input_mask_hydro_buffer, points_skeleton_gdf, classes, k + river_length = config.virtual_point.vector.river_length + points_grid_spacing = config.virtual_point.pointcloud.points_grid_spacing + classes = config.virtual_point.pointcloud.virtual_points_classes + + # Step 1 : Merged all "points around skeleton" by lidar tile + def process_points_knn(points_knn): + # Check if points_knn is a string and convert it to a list if necessary + if isinstance(points_knn, str): + points_knn = ast.literal_eval(points_knn) # Convert the string to a list of lists + return [[round(coord, 3) for coord in point] for point in points_knn] + + points_clip_list = [ + {"geometry": row["geometry"], "points_knn": process_points_knn(row["points_knn"])} + for filename in os.listdir(input_dir_points_skeleton) + if filename.endswith(".geojson") + for _, row in gpd.read_file(os.path.join(input_dir_points_skeleton, filename)).iterrows() + ] + # List match Z elevation values every N meters along the hydrographic skeleton + df = pd.DataFrame(points_clip_list) + + # Step 2: Combine skeleton lines into a single polyline for each hydro entity + if not df.empty and "points_knn" in df.columns and "geometry" in df.columns: + points_gdf = gpd.GeoDataFrame(df, geometry="geometry") + points_gdf.set_crs(crs, inplace=True) + # 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 = [ + launch_virtual_points_by_section( + 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 file in os.listdir(input_dir_points) + for idx, row in gdf_merged.iterrows() ] - # Flatten the list of lists into a single list of dictionaries - points_clip = [item for sublist in points_clip_list for item in sublist] - - # Create a pandas DataFrame from the flattened list - df = pd.DataFrame(points_clip) - # Create a GeoDataFrame from the pandas DataFrame - result_gdf = gpd.GeoDataFrame(df, geometry="geometry") - result_gdf.set_crs(crs, inplace=True) - # path to the Points along Skeleton Hydro file - output_file = os.path.join(output_dir, "Points_Skeleton.GeoJSON") - # Save to GeoJSON - result_gdf.to_file(output_file, driver="GeoJSON") + logging.info("Calculate virtuals points by mask hydro and skeleton") + # Step 4 : 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") if __name__ == "__main__": diff --git a/lidro/main_extract_points_from_skeleton.py b/lidro/main_extract_points_from_skeleton.py new file mode 100644 index 00000000..07df6f63 --- /dev/null +++ b/lidro/main_extract_points_from_skeleton.py @@ -0,0 +1,84 @@ +""" Main script for creating N points along Skeleton +""" + +import logging +import os + +import hydra +from omegaconf import DictConfig +from pyproj import CRS + +from lidro.create_virtual_point.vectors.extract_points_around_skeleton import ( + extract_points_around_skeleton_points_one_tile, +) +from lidro.create_virtual_point.vectors.mask_hydro_with_buffer import ( + import_mask_hydro_with_buffer, +) +from lidro.create_virtual_point.vectors.points_along_skeleton import ( + generate_points_along_skeleton, +) + + +@hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") +def main(config: DictConfig): + """Create N points along Skeleton from the points classification of + the input LAS/LAZ file and the Hydro Skeleton (GeoJSON). Save a result by LIDAR tiles. + + It can run either on a single file, or 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) + + # If input filename is not provided, lidro runs on the whole input_dir directory + initial_las_filename = config.io.input_filename + input_mask_hydro = config.io.input_mask_hydro + input_skeleton = config.io.input_skeleton + + # Parameters for creating virtual point + distance_meters = config.virtual_point.vector.distance_meters + buffer = config.virtual_point.vector.buffer + crs = CRS.from_user_input(config.io.srid) + classes = config.virtual_point.filter.keep_neighbors_classes + k = config.virtual_point.vector.k + + # Step 1 : Import Mask Hydro, then apply a buffer + # Return GeoDataframe + input_mask_hydro_buffer = import_mask_hydro_with_buffer(input_mask_hydro, buffer, crs).wkt + + # Step 2 : Create several points every 2 meters (by default) along skeleton Hydro + # Return GeoDataframe + points_skeleton_gdf = generate_points_along_skeleton(input_skeleton, distance_meters, crs) + + # Step 3 : Extract points around skeleton by tile + if initial_las_filename: + # Lauch croping filtered pointcloud by Mask Hydro with buffer by one tile: + extract_points_around_skeleton_points_one_tile( + initial_las_filename, input_dir, output_dir, input_mask_hydro_buffer, points_skeleton_gdf, crs, classes, k + ) + else: + # Lauch croping filtered pointcloud by Mask Hydro with buffer tile by tile + input_dir_points = os.path.join(input_dir, "pointcloud") + for file in os.listdir(input_dir_points): + extract_points_around_skeleton_points_one_tile( + file, input_dir, output_dir, input_mask_hydro_buffer, points_skeleton_gdf, crs, classes, k + ) + + +if __name__ == "__main__": + main() diff --git a/scripts/example_create_mask_default.sh b/scripts/example_create_mask_default.sh index 8818c5b4..1287c98c 100644 --- a/scripts/example_create_mask_default.sh +++ b/scripts/example_create_mask_default.sh @@ -2,6 +2,7 @@ python -m lidro.main_create_mask \ io.input_dir=./data/pointcloud/ \ io.output_dir=./tmp/ \ - +io.mask_generation.raster.dilatation_size=3 \ +io.mask_generation.filter.keep_classes="[0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67]" \ diff --git a/scripts/example_create_virtual_point_by_tile.sh b/scripts/example_create_virtual_point_by_tile.sh index 4a1a8ccc..4c48facd 100644 --- a/scripts/example_create_virtual_point_by_tile.sh +++ b/scripts/example_create_virtual_point_by_tile.sh @@ -5,6 +5,6 @@ io.input_filename=Semis_2021_0830_6291_LA93_IGN69.laz \ io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geojson \ io.input_skeleton=./data/skeleton_hydro/Skeleton_Hydro.geojson \ io.output_dir=./tmp/ \ - +io.virtual_point.pointcloud.points_grid_spacing=0.5 \ diff --git a/scripts/example_create_virtual_point_default.sh b/scripts/example_create_virtual_point_default.sh index 3622ffbf..b37bab10 100644 --- a/scripts/example_create_virtual_point_default.sh +++ b/scripts/example_create_virtual_point_default.sh @@ -3,7 +3,9 @@ python -m lidro.main_create_virtual_point \ io.input_dir=./data/ \ io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geojson \ io.input_skeleton=./data/skeleton_hydro/Skeleton_Hydro.geojson \ +io.dir_points_skeleton=./tmp/point_skeleton/ \ io.output_dir=./tmp/ \ +io.virtual_point.pointcloud.points_grid_spacing=0.5 \ diff --git a/scripts/example_extract_points_around_skeleton_default.sh b/scripts/example_extract_points_around_skeleton_default.sh new file mode 100644 index 00000000..11913cd8 --- /dev/null +++ b/scripts/example_extract_points_around_skeleton_default.sh @@ -0,0 +1,13 @@ +# Launch hydro mask merging +python -m lidro.main_extract_points_from_skeleton \ +io.input_dir=./data/ \ +io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geosjon \ +io.input_skeleton=./data/skeleton_hydro/Skeleton_Hydro.geojson \ +io.output_dir=./tmp/points_skeleton/ \ +io.virtual_point.vector.distance_meters=5 \ +io.virtual_point.vector.buffer=3 \ +io.virtual_point.vector.k=3 \ +io.virtual_point.filter.keep_neighbors_classes="[2, 9]" \ + + + diff --git a/scripts/example_merge_mask_default.sh b/scripts/example_merge_mask_default.sh index 7e162fc7..0d0d1165 100644 --- a/scripts/example_merge_mask_default.sh +++ b/scripts/example_merge_mask_default.sh @@ -2,6 +2,10 @@ python -m lidro.main_merge_mask \ io.input_dir=./data/mask_hydro/ \ io.output_dir=./tmp/merge_mask_hydro/ \ +io.mask_generation.vector.min_water_area=150 \ +io.mask_generation.vector.buffer_positive=1 \ +io.mask_generation.vector.buffer_negative=-1.5 \ +io.mask_generation.vector.tolerance=1 \ diff --git a/test/pointcloud/test_convert_list_points_to_las.py b/test/pointcloud/test_convert_list_points_to_las.py new file mode 100644 index 00000000..50d2d9bc --- /dev/null +++ b/test/pointcloud/test_convert_list_points_to_las.py @@ -0,0 +1,68 @@ +import os +import shutil +from pathlib import Path + +import geopandas as gpd +import laspy +import numpy as np +import pyproj +from shapely.geometry import Point + +# Import the function from your module +from lidro.create_virtual_point.pointcloud.convert_list_points_to_las import ( + list_points_to_las, +) + +TMP_PATH = Path("./tmp/virtual_points/pointcloud/convert_geodataframe_to_las") + + +def setup_module(module): + if TMP_PATH.is_dir(): + shutil.rmtree(TMP_PATH) + os.makedirs(TMP_PATH) + + +def test_list_points_to_las_default(): + # Create a sample GeoDataFrame with virtual points (EPSG:2154) + points = gpd.GeoDataFrame( + { + "geometry": [ + Point(700000, 6600000, 10), + Point(700001, 6600001, 15), + Point(700002, 6600002, 20), + Point(700010, 6600010, 25), + ] + }, + crs="EPSG:2154", + ) + + # Define the CRS + crs = pyproj.CRS("EPSG:2154") + + # Call the function to test + list_points_to_las([points], TMP_PATH, crs, 68) + + # Verify the output LAS file + output_las = os.path.join(TMP_PATH, "virtual_points.laz") + assert os.path.exists(output_las), "The output LAS file should exist" + + # Read the LAS file + las = laspy.read(output_las) + + # Check that the coordinates match the input points + expected_coords = np.array([[p.x, p.y, p.z] for p in points.geometry]) + actual_coords = np.vstack((las.x, las.y, las.z)).T + np.testing.assert_array_almost_equal( + actual_coords, + expected_coords, + decimal=5, + err_msg="The coordinates in the LAS file do not match the input points", + ) + + # Check that the classification is 68 for all points + expected_classification = np.full(len(points), 68, dtype=np.uint8) + np.testing.assert_array_equal( + las.classification, + expected_classification, + err_msg="The classification in the LAS file should be 68 for all points", + ) diff --git a/test/test_main_create_virtual_point.py b/test/test_main_create_virtual_point.py index 27401864..add4a2ce 100644 --- a/test/test_main_create_virtual_point.py +++ b/test/test_main_create_virtual_point.py @@ -2,7 +2,6 @@ import subprocess as sp from pathlib import Path -import pytest from hydra import compose, initialize from lidro.main_create_virtual_point import main @@ -21,51 +20,23 @@ def test_main_run_okay(): io.input_dir="{repo_dir}/lidro/data/"\ io.input_filename=Semis_2021_0830_6291_LA93_IGN69.laz \ io.input_mask_hydro="{repo_dir}/lidro/data/merge_mask_hydro/MaskHydro_merge.geojson"\ - io.input_skeleton="{repo_dir}/lidro/data/skeleton_hydro/Skeleton_Hydro.geojson"\ + io.input_skeleton="{repo_dir}/lidro/data/skeleton_hydro/Skeleton_Hydro_small.geojson"\ + io.dir_points_skeleton="{repo_dir}/lidro/data/point_virtual/"\ io.output_dir="{repo_dir}/lidro/tmp/create_virtual_point/main/" """ sp.run(cmd, shell=True, check=True) -def test_main_lidro_default(): - input_dir = INPUT_DIR - output_dir = OUTPUT_DIR / "main_lidro_default" - input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" - input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson" - distances_meters = 10 - buffer = 10 - srid = 2154 - k = 20 - - 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.output_dir={output_dir}", - f"io.input_mask_hydro={input_mask_hydro}", - f"io.input_skeleton={input_skeleton}", - f"virtual_point.vector.distance_meters={distances_meters}", - f"virtual_point.vector.buffer={buffer}", - f"virtual_point.vector.k={k}", - f"io.srid={srid}", - ], - ) - main(cfg) - assert (Path(output_dir) / "Points_Skeleton.GeoJSON").is_file() - - def test_main_lidro_input_file(): input_dir = INPUT_DIR output_dir = OUTPUT_DIR / "main_lidro_input_file" input_filename = "Semis_2021_0830_6291_LA93_IGN69.laz" input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" - input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson" - distances_meters = 5 - buffer = 2 + input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro_small.geojson" + dir_points_skeleton = INPUT_DIR / "point_virtual/" srid = 2154 - k = 10 + points_grid_spacing = 1 + with initialize(version_base="1.2", config_path="../configs"): # config is relative to a module cfg = compose( @@ -76,142 +47,10 @@ def test_main_lidro_input_file(): f"io.output_dir={output_dir}", f"io.input_mask_hydro={input_mask_hydro}", f"io.input_skeleton={input_skeleton}", - f"virtual_point.vector.distance_meters={distances_meters}", - f"virtual_point.vector.buffer={buffer}", - f"virtual_point.vector.k={k}", + f"io.dir_points_skeleton={dir_points_skeleton}", f"io.srid={srid}", + f"virtual_point.pointcloud.points_grid_spacing={points_grid_spacing}", ], ) main(cfg) - assert (Path(output_dir) / "Points_Skeleton.GeoJSON").is_file() - - -def test_main_lidro_fail_no_input_dir(): - output_dir = OUTPUT_DIR / "main_no_input_dir" - input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" - input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson" - distances_meters = 5 - buffer = 2 - srid = 2154 - k = 10 - 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_mask_hydro={input_mask_hydro}", - f"io.input_skeleton={input_skeleton}", - f"io.output_dir={output_dir}", - f"virtual_point.vector.distance_meters={distances_meters}", - f"virtual_point.vector.buffer={buffer}", - f"virtual_point.vector.k={k}", - 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" - input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" - input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson" - distances_meters = 5 - buffer = 2 - srid = 2154 - k = 10 - 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.input_mask_hydro={input_mask_hydro}", - f"io.input_skeleton={input_skeleton}", - f"io.output_dir={output_dir}", - f"virtual_point.vector.distance_meters={distances_meters}", - f"virtual_point.vector.buffer={buffer}", - f"virtual_point.vector.k={k}", - f"io.srid={srid}", - ], - ) - with pytest.raises(FileNotFoundError): - main(cfg) - - -def test_main_lidro_fail_no_output_dir(): - input_dir = INPUT_DIR - input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" - input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson" - distances_meters = 5 - buffer = 2 - srid = 2154 - k = 10 - 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.input_mask_hydro={input_mask_hydro}", - f"io.input_skeleton={input_skeleton}", - f"virtual_point.vector.distance_meters={distances_meters}", - f"virtual_point.vector.buffer={buffer}", - f"virtual_point.vector.k={k}", - f"io.srid={srid}", - ], - ) - with pytest.raises(ValueError): - main(cfg) - - -def test_main_lidro_fail_no_input_mask_hydro(): - input_dir = INPUT_DIR - output_dir = OUTPUT_DIR / "main_no_input_dir" - input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro.geojson" - distances_meters = 5 - buffer = 2 - srid = 2154 - k = 10 - with initialize(version_base="1.2", config_path="../configs"): - cfg = compose( - config_name="configs_lidro", - overrides=[ - f"io.input_dir={input_dir}", - f"io.output_dir={output_dir}", - "io.input_mask_hydro=some_random_input_mask_hydro", - f"io.input_skeleton={input_skeleton}", - f"virtual_point.vector.distance_meters={distances_meters}", - f"virtual_point.vector.buffer={buffer}", - f"virtual_point.vector.k={k}", - f"io.srid={srid}", - ], - ) - with pytest.raises(ValueError): - main(cfg) - - -def test_main_lidro_fail_no_input_skeleton(): - input_dir = INPUT_DIR - output_dir = OUTPUT_DIR / "main_no_input_dir" - input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" - distances_meters = 5 - buffer = 2 - srid = 2154 - k = 10 - with initialize(version_base="1.2", config_path="../configs"): - cfg = compose( - config_name="configs_lidro", - overrides=[ - f"io.input_dir={input_dir}", - f"io.output_dir={output_dir}", - f"io.input_mask_hydro={input_mask_hydro}", - "io.input_skeleton=some_random_input_skeleton", - f"virtual_point.vector.distance_meters={distances_meters}", - f"virtual_point.vector.buffer={buffer}", - f"virtual_point.vector.k={k}", - f"io.srid={srid}", - ], - ) - with pytest.raises(ValueError): - main(cfg) + assert (Path(output_dir) / "virtual_points.laz").is_file() diff --git a/test/test_main_extract_points_from_skeleton.py b/test/test_main_extract_points_from_skeleton.py new file mode 100644 index 00000000..dff01e6d --- /dev/null +++ b/test/test_main_extract_points_from_skeleton.py @@ -0,0 +1,188 @@ +import os +import subprocess as sp +from pathlib import Path + +import pytest +from hydra import compose, initialize + +from lidro.main_extract_points_from_skeleton import main + +INPUT_DIR = Path("data") +OUTPUT_DIR = Path("tmp") / "extract_points_around_skeleton/main" + + +def setup_module(module): + os.makedirs("tmp/extract_points_around_skeleton/main", exist_ok=True) + + +def test_main_run_okay(): + repo_dir = Path.cwd().parent + cmd = f"""python -m lidro.main_extract_points_from_skeleton \ + io.input_dir="{repo_dir}/lidro/data/"\ + io.input_filename=Semis_2021_0830_6291_LA93_IGN69.laz \ + io.input_mask_hydro="{repo_dir}/lidro/data/merge_mask_hydro/MaskHydro_merge.geojson"\ + io.input_skeleton="{repo_dir}/lidro/data/skeleton_hydro/Skeleton_Hydro_small.geojson"\ + io.output_dir="{repo_dir}/lidro/tmp/extract_points_around_skeleton/main/" + """ + sp.run(cmd, shell=True, check=True) + + +def test_main_extract_points_skeleton_input_file(): + input_dir = INPUT_DIR + output_dir = OUTPUT_DIR / "main_extract_points_skeleton_input_file" + input_filename = "Semis_2021_0830_6291_LA93_IGN69.laz" + input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" + input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro_small.geojson" + distances_meters = 5 + buffer = 2 + srid = 2154 + k = 10 + 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_filename={input_filename}", + f"io.input_dir={input_dir}", + f"io.output_dir={output_dir}", + f"io.input_mask_hydro={input_mask_hydro}", + f"io.input_skeleton={input_skeleton}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", + f"io.srid={srid}", + ], + ) + main(cfg) + assert (Path(output_dir) / "Semis_2021_0830_6291_LA93_IGN69_points_skeleton.geojson").is_file() + + +def test_main_extract_points_skeleton_fail_no_input_dir(): + output_dir = OUTPUT_DIR / "main_no_input_dir" + input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" + input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro_small.geojson" + distances_meters = 5 + buffer = 2 + srid = 2154 + k = 10 + 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_mask_hydro={input_mask_hydro}", + f"io.input_skeleton={input_skeleton}", + f"io.output_dir={output_dir}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", + f"io.srid={srid}", + ], + ) + with pytest.raises(ValueError): + main(cfg) + + +def test_main_extract_points_skeleton_fail_wrong_input_dir(): + output_dir = OUTPUT_DIR / "main_wrong_input_dir" + input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" + input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro_small.geojson" + distances_meters = 5 + buffer = 2 + srid = 2154 + k = 10 + 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.input_mask_hydro={input_mask_hydro}", + f"io.input_skeleton={input_skeleton}", + f"io.output_dir={output_dir}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", + f"io.srid={srid}", + ], + ) + with pytest.raises(FileNotFoundError): + main(cfg) + + +def test_main_extract_points_skeleton_fail_no_output_dir(): + input_dir = INPUT_DIR + input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" + input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro_small.geojson" + distances_meters = 5 + buffer = 2 + srid = 2154 + k = 10 + 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.input_mask_hydro={input_mask_hydro}", + f"io.input_skeleton={input_skeleton}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", + f"io.srid={srid}", + ], + ) + with pytest.raises(ValueError): + main(cfg) + + +def test_main_extract_points_skeleton_fail_no_input_mask_hydro(): + input_dir = INPUT_DIR + output_dir = OUTPUT_DIR / "main_no_input_dir" + input_skeleton = INPUT_DIR / "skeleton_hydro/Skeleton_Hydro_small.geojson" + distances_meters = 5 + buffer = 2 + srid = 2154 + k = 10 + with initialize(version_base="1.2", config_path="../configs"): + cfg = compose( + config_name="configs_lidro", + overrides=[ + f"io.input_dir={input_dir}", + f"io.output_dir={output_dir}", + "io.input_mask_hydro=some_random_input_mask_hydro", + f"io.input_skeleton={input_skeleton}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", + f"io.srid={srid}", + ], + ) + with pytest.raises(ValueError): + main(cfg) + + +def test_main_extract_points_skeleton_fail_no_input_skeleton(): + input_dir = INPUT_DIR + output_dir = OUTPUT_DIR / "main_no_input_dir" + input_mask_hydro = INPUT_DIR / "merge_mask_hydro/MaskHydro_merge.geojson" + distances_meters = 5 + buffer = 2 + srid = 2154 + k = 10 + with initialize(version_base="1.2", config_path="../configs"): + cfg = compose( + config_name="configs_lidro", + overrides=[ + f"io.input_dir={input_dir}", + f"io.output_dir={output_dir}", + f"io.input_mask_hydro={input_mask_hydro}", + "io.input_skeleton=some_random_input_skeleton", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", + f"io.srid={srid}", + ], + ) + with pytest.raises(ValueError): + main(cfg) diff --git a/test/vectors/test_apply_z_from_grid.py b/test/vectors/test_apply_z_from_grid.py new file mode 100644 index 00000000..c5e14a82 --- /dev/null +++ b/test/vectors/test_apply_z_from_grid.py @@ -0,0 +1,60 @@ +import geopandas as gpd +import numpy as np +from shapely.geometry import LineString, Point + +from lidro.create_virtual_point.vectors.apply_Z_from_grid import ( + calculate_grid_z, + calculate_grid_z_with_model, +) + + +def test_calculate_grid_z_with_model(): + # Create a sample GeoDataFrame of points + points = gpd.GeoDataFrame({"geometry": [Point(0, 0), Point(1, 1), Point(2, 2)]}, crs="EPSG:4326") + + # Create a sample GeoDataFrame of line + line = gpd.GeoDataFrame({"geometry": [LineString([(0, 0), (2, 2)])]}, crs="EPSG:4326") + + # Sample model function + def model(x): + return np.array(x) * 2 # Simple linear model for testing + + # Call the function to test + result = calculate_grid_z_with_model(points, line, model) + + # Check that the result is a GeoDataFrame + assert isinstance(result, gpd.GeoDataFrame), "The result should be a GeoDataFrame" + + # Check that the result has the same number of points + assert len(result) == len(points), "The number of points should be the same" + + # Check the Z values + curvilinear_abs = [0, np.sqrt(2), np.sqrt(8)] + expected_z = model(curvilinear_abs) + result_z = result.geometry.apply(lambda geom: geom.z) + + # Use assert_array_almost_equal to check the values + np.testing.assert_array_almost_equal( + result_z, expected_z, decimal=5, err_msg="Z values do not match the expected values" + ) + + +def test_calculate_grid_z_for_flattening(): + # Create a sample GeoDataFrame of points + points = gpd.GeoDataFrame({"geometry": [Point(0, 0), Point(1, 1), Point(2, 2)]}, crs="EPSG:4326") + + # Predicted Z for flattening + predicted_z = 10.0 + + # Call the function to test + result = calculate_grid_z(points, predicted_z) + + # Check that the result is a GeoDataFrame + assert isinstance(result, gpd.GeoDataFrame), "The result should be a GeoDataFrame" + + # Check that the result has the same number of points + assert len(result) == len(points), "The number of points should be the same" + + # Check the Z values + result_z = result.geometry.apply(lambda geom: geom.z) + assert (result_z == predicted_z).all(), "All Z values should be equal to predicted_z" diff --git a/test/vectors/test_create_grid_2D_inside_maskhydro.py b/test/vectors/test_create_grid_2D_inside_maskhydro.py new file mode 100644 index 00000000..a134bc21 --- /dev/null +++ b/test/vectors/test_create_grid_2D_inside_maskhydro.py @@ -0,0 +1,44 @@ +import geopandas as gpd +from shapely.geometry import Polygon + +from lidro.create_virtual_point.vectors.create_grid_2D_inside_maskhydro import ( + generate_grid_from_geojson, +) + + +def test_generate_grid_from_geojson_default(): + # Create a sample GeoDataFrame with a smaller polygon (200m x 200m) in EPSG:2154 + polygon = Polygon([(700000, 6600000), (700000, 6600200), (700200, 6600200), (700200, 6600000), (700000, 6600000)]) + mask_hydro = gpd.GeoDataFrame({"geometry": [polygon]}, crs="EPSG:2154") + + # Define the spacing and margin + spacing = 0.5 # 0.5 meter spacing + margin = 0 # no marging + + # Call the function to test + result = generate_grid_from_geojson(mask_hydro, spacing, margin) + + # Check that the result is a GeoDataFrame + assert isinstance(result, gpd.GeoDataFrame), "The result should be a GeoDataFrame" + + # Check that the points are within the polygon bounds + assert result.within(polygon).all(), "All points should be within the polygon" + + +def test_generate_grid_from_geojson(): + # Create a sample GeoDataFrame with a smaller polygon (200m x 200m) in EPSG:2154 + polygon = Polygon([(700000, 6600000), (700000, 6600200), (700200, 6600200), (700200, 6600000), (700000, 6600000)]) + mask_hydro = gpd.GeoDataFrame({"geometry": [polygon]}, crs="EPSG:2154") + + # Define the spacing and margin + spacing = 2 # 2 meter spacing + margin = 1 # 1 meter marging + + # Call the function to test + result = generate_grid_from_geojson(mask_hydro, spacing, margin) + + # Check that the result is a GeoDataFrame + assert isinstance(result, gpd.GeoDataFrame), "The result should be a GeoDataFrame" + + # Check that the points are within the polygon bounds + assert result.within(polygon).all(), "All points should be within the polygon" diff --git a/test/vectors/test_extract_points_around_skeleton.py b/test/vectors/test_extract_points_around_skeleton.py new file mode 100644 index 00000000..5baa06b3 --- /dev/null +++ b/test/vectors/test_extract_points_around_skeleton.py @@ -0,0 +1,68 @@ +import os +import shutil +from pathlib import Path + +import geopandas as gpd +from pyproj import CRS +from shapely.geometry import Point + +from lidro.create_virtual_point.vectors.extract_points_around_skeleton import ( + extract_points_around_skeleton_points_one_tile, +) + +TMP_PATH = Path("./tmp/create_virtual_point/vectors/extract_points_around_skeleton") +INPUT_DIR = Path("./data/") +MASK_HYDRO = "./data/merge_mask_hydro/MaskHydro_merge.geojson" +POINTS_SKELETON = "./data/skeleton_hydro/points_along_skeleton_small.geojson" +OUTPUT_GEOJSON = TMP_PATH / "Semis_2021_0830_6291_LA93_IGN69_points_skeleton.geojson" + + +def setup_module(module): + if TMP_PATH.is_dir(): + shutil.rmtree(TMP_PATH) + os.makedirs(TMP_PATH) + + +def test_extract_points_around_skeleton_default(): + # Parameters + crs = CRS.from_epsg(2154) + classes = [2] # Example classes + k = 3 # Example k value + + # Mask Hydro with buffer + mask_hydro_buffered = "MULTIPOLYGON (((829969.0856167737 6292032.553442742,\ + 830452.9506643447 6292032.553442742, \ + 830498.1716968281 6291675.307286125, \ + 830624.7905877812 6291073.867554097, \ + 830783.0642014727 6290675.922468244, \ + 830972.9925379024 6290400.074170096, \ + 831045.3461898757 6290228.23424666, \ + 831036.3019833791 6289952.385948512, \ + 830783.0642014727 6289947.863845264, \ + 830371.5528058748 6290599.046713023, \ + 830055.005578492 6290947.248663144, \ + 829919.3424810421 6291399.458987976, \ + 829969.0856167737 6292032.553442742)))" + + # Execute the function + extract_points_around_skeleton_points_one_tile( + "Semis_2021_0830_6291_LA93_IGN69.laz", + INPUT_DIR, + TMP_PATH, + mask_hydro_buffered, + gpd.read_file(POINTS_SKELETON), + crs, + classes, + k, + ) + + # Verify that the output file exists + assert OUTPUT_GEOJSON.exists() + + # Load the output GeoJSON + gdf = gpd.read_file(OUTPUT_GEOJSON) + + # Assertions to check the integrity of the output + assert not gdf.empty, "GeoDataFrame shouldn't be empty" + assert gdf.crs.to_string() == crs.to_string(), "CRS should match the specified CRS" + assert all(isinstance(geom, Point) for geom in gdf.geometry), "All geometries should be Points" diff --git a/test/vectors/test_flatten_river.py b/test/vectors/test_flatten_river.py new file mode 100644 index 00000000..02a59040 --- /dev/null +++ b/test/vectors/test_flatten_river.py @@ -0,0 +1,52 @@ +import geopandas as gpd +import numpy as np +from shapely.geometry import LineString, Point + +from lidro.create_virtual_point.vectors.flatten_river import flatten_little_river + + +def test_flatten_little_river_points_knn_not_empty(): + # Create example points GeoDataFrame + points = gpd.GeoDataFrame( + { + "geometry": [Point(0, 0, 10), Point(1, 1, 20), Point(2, 2, 30)], + "points_knn": [np.array([[0, 0, 10], [1, 1, 20], [2, 2, 30]]) for _ in range(3)], + } + ) + + # Create example line GeoDataFrame + line = gpd.GeoDataFrame({"geometry": [LineString([(0, 0), (2, 2)])]}) + + # Test when points_knn is not empty + expected_result = 15.0 # Adjusted to the correct 1st quartile value + result = flatten_little_river(points, line) + assert result == expected_result, f"Expected {expected_result}, but got {result}" + + +def test_flatten_little_river_points_knn_empty(): + # Test when points_knn is empty + empty_points = gpd.GeoDataFrame({"geometry": [Point(0, 0)], "points_knn": [np.array([])]}) + + # Create example line GeoDataFrame + line = gpd.GeoDataFrame({"geometry": [LineString([(0, 0), (2, 2)])]}) + + expected_result = 0 + result = flatten_little_river(empty_points, line) + assert result == expected_result, f"Expected {expected_result}, but got {result}" + + +def test_flatten_little_river_duplicate_points(): + # Test with duplicate points + duplicate_points = gpd.GeoDataFrame( + { + "geometry": [Point(0, 0, 10), Point(1, 1, 10), Point(2, 2, 10)], + "points_knn": [np.array([[0, 0, 10], [1, 1, 10], [2, 2, 10]]) for _ in range(3)], + } + ) + + # Create example line GeoDataFrame + line = gpd.GeoDataFrame({"geometry": [LineString([(0, 0), (2, 2)])]}) + + expected_result = 10 + result = flatten_little_river(duplicate_points, line) + assert result == expected_result, f"Expected {expected_result}, but got {result}" diff --git a/test/vectors/test_intersection_points_by_line.py b/test/vectors/test_intersection_points_by_line.py new file mode 100644 index 00000000..00582009 --- /dev/null +++ b/test/vectors/test_intersection_points_by_line.py @@ -0,0 +1,34 @@ +import geopandas as gpd +from shapely.geometry import LineString, Point + +# Import the functions from your module +from lidro.create_virtual_point.vectors.intersect_points_by_line import ( + return_points_by_line, +) + + +def test_return_points_by_line_default(): + # Create a sample GeoDataFrame with points along a river (EPSG:2154) + points = gpd.GeoDataFrame( + {"geometry": [Point(700000, 6600000), Point(700001, 6600001), Point(700002, 6600002), Point(700010, 6600010)]}, + crs="EPSG:2154", + ) + + # Create a sample GeoDataFrame with a line representing the river (EPSG:2154) + line = gpd.GeoDataFrame({"geometry": [LineString([(700000, 6600000), (700002, 6600002)])]}, crs="EPSG:2154") + + # Call the function to test + result = return_points_by_line(points, line) + result_geom = result["geometry"] + + # Expected points are those within the buffer of the line + expected_points = gpd.GeoDataFrame( + {"geometry": [Point(700000, 6600000), Point(700001, 6600001), Point(700002, 6600002)]}, crs="EPSG:2154" + ) + expected_points_geom = expected_points["geometry"] + + # Check that the result is a GeoDataFrame + assert isinstance(result, gpd.GeoDataFrame), "The result should be a GeoDataFrame" + + # Check that the result contains the expected points + assert result_geom.equals(expected_points_geom), f"The result does not match the expected points. Got: {result}" diff --git a/test/vectors/test_las_around_point.py b/test/vectors/test_las_around_point.py index 5cdfb756..214bf239 100644 --- a/test/vectors/test_las_around_point.py +++ b/test/vectors/test_las_around_point.py @@ -4,6 +4,7 @@ import geopandas as gpd import numpy as np +import pandas as pd from pyproj import CRS from shapely.geometry import Point @@ -40,14 +41,12 @@ def test_las_around_point_default(): result = filter_las_around_point(points_along_skeleton, points_clip, k) - # Convert results to GeoDataFrame - result_gdf = gpd.GeoDataFrame(result) - result_gdf.set_crs(crs, inplace=True) - # Save to GeoJSON - result_gdf.to_file(output, driver="GeoJSON") + # Create a pandas DataFrame from the flattened list + df = pd.DataFrame(result) + # Create a GeoDataFrame from the pandas DataFrame + points_gdf = gpd.GeoDataFrame(df, geometry="geometry") + points_gdf.set_crs(crs, inplace=True) - gdf = gpd.read_file(output) - - assert not gdf.empty # GeoDataFrame shouldn't empty - assert gdf.crs.to_string() == crs # CRS is identical - assert all(isinstance(geom, Point) for geom in gdf.geometry) # All geometry should Points + assert not points_gdf.empty # GeoDataFrame shouldn't empty + assert points_gdf.crs.to_string() == crs # CRS is identical + assert all(isinstance(geom, Point) for geom in points_gdf.geometry) # All geometry should Points diff --git a/test/vectors/test_linear_regression_model.py b/test/vectors/test_linear_regression_model.py new file mode 100644 index 00000000..804340f5 --- /dev/null +++ b/test/vectors/test_linear_regression_model.py @@ -0,0 +1,101 @@ +import geopandas as gpd +import numpy as np +from shapely.geometry import LineString, Point +from sklearn.neighbors import NearestNeighbors + +from lidro.create_virtual_point.vectors.linear_regression_model import ( + calculate_linear_regression_line, +) + + +def create_knn_column(points): + if points.empty: + return points + + # Extracting the 3D points + coords = np.array([[p.x, p.y, p.z] for p in points.geometry]) + + # Using Nearest Neighbors to find the 3 nearest neighbors + nbrs = NearestNeighbors(n_neighbors=min(3, len(points)), algorithm="ball_tree").fit(coords) + distances, indices = nbrs.kneighbors(coords) + + # Creating the 'points_knn' column + points_knn = [coords[indices[i]] for i in range(len(points))] + points["points_knn"] = points_knn + return points + + +def test_calculate_linear_regression_line_default(): + # Create a sample GeoDataFrame with points designed to produce an R² of 0.25 + np.random.seed(0) + x_coords = np.linspace(0, 100, num=10) + y_coords = np.linspace(0, 100, num=10) + z_coords = x_coords * 0.1 + np.random.normal(0, 5, 10) # Linear relation with noise + + points = gpd.GeoDataFrame( + {"geometry": [Point(x, y, z) for x, y, z in zip(x_coords, y_coords, z_coords)]}, crs="EPSG:2154" + ) + + # Add the points_knn column + points = create_knn_column(points) + + # Create a sample GeoDataFrame with a line representing the river (EPSG:2154) + line = gpd.GeoDataFrame({"geometry": [LineString([(0, 0), (100, 100)])]}, crs="EPSG:2154") + + # Define the CRS + crs = "EPSG:2154" + + # Call the function to test + model, r2 = calculate_linear_regression_line(points, line, crs) + + # Check that the r2 is within the range [0, 1] + assert 0 <= r2 <= 1, "The determination coefficient should be between 0 and 1" + + # Check the type of the model + assert isinstance(model, np.poly1d), "The regression model should be of type np.poly1d" + + # Check the length of the model's coefficients + assert len(model.coefficients) == 2, "The regression model should be linear with 2 coefficients" + + # Check the model's coefficients are within expected ranges + expected_slope_range = (0.01, 0.15) # Adjusted range + expected_intercept_range = (-10, 10) + slope, intercept = model.coefficients + assert ( + expected_slope_range[0] <= slope <= expected_slope_range[1] + ), f"Slope {slope} is not within the expected range {expected_slope_range}" + assert ( + expected_intercept_range[0] <= intercept <= expected_intercept_range[1] + ), f"Intercept {intercept} is not within the expected range {expected_intercept_range}" + print("All tests passed successfully!") + + +def test_calculate_linear_regression_line_empty_inputs(): + # Test with empty GeoDataFrames + points = gpd.GeoDataFrame(columns=["geometry"], crs="EPSG:2154") + line = gpd.GeoDataFrame(columns=["geometry"], crs="EPSG:2154") + crs = "EPSG:2154" + + model, r2 = calculate_linear_regression_line(points, line, crs) + assert isinstance(model, np.poly1d), "The regression model should be of type np.poly1d" + assert r2 == 0.0, "The determination coefficient should be 0 for empty inputs" + print("Empty input test passed successfully!") + + +def test_calculate_linear_regression_line_single_point(): + # Test with a single point + points = gpd.GeoDataFrame({"geometry": [Point(0, 0, 0)]}, crs="EPSG:2154") + + # Add the points_knn column + points = create_knn_column(points) + + # Create a sample GeoDataFrame with a line representing the river (EPSG:2154) + line = gpd.GeoDataFrame({"geometry": [LineString([(0, 0), (100, 100)])]}, crs="EPSG:2154") + + # Define the CRS + crs = "EPSG:2154" + + model, r2 = calculate_linear_regression_line(points, line, crs) + assert isinstance(model, np.poly1d), "The regression model should be of type np.poly1d" + assert r2 == 0.0, "The determination coefficient should be 0 for a single point input" + print("Single point input test passed successfully!") diff --git a/test/vectors/test_merge_skeleton_by_mask.py b/test/vectors/test_merge_skeleton_by_mask.py new file mode 100644 index 00000000..43640e8a --- /dev/null +++ b/test/vectors/test_merge_skeleton_by_mask.py @@ -0,0 +1,36 @@ +import os +import shutil +from pathlib import Path + +import pandas as pd +from pyproj import CRS + +from lidro.create_virtual_point.vectors.merge_skeleton_by_mask import ( + merge_skeleton_by_mask, +) + +TMP_PATH = Path("./tmp/virtual_points/vectors/merge_skeleton_by_mask") + +mask = Path("./data/merge_mask_hydro/MaskHydro_merge.geojson") +skeleton = Path("./data/skeleton_hydro/Skeleton_Hydro.geojson") +output = Path("./tmp/virtual_points/vectors/merge_skeleton_by_mask/merge_skeleton_by_mask.geojson") + + +def setup_module(module): + if TMP_PATH.is_dir(): + shutil.rmtree(TMP_PATH) + os.makedirs(TMP_PATH) + + +def test_merge_skeleton_by_mask_default(): + # Parameters + crs = CRS.from_epsg(2154) + + result = merge_skeleton_by_mask(skeleton, mask, TMP_PATH, crs) + + # Check if the result is a DataFrame + assert isinstance(result, pd.DataFrame), "Result is not a DataFrame" + + # Check if the necessary columns are present + assert "geometry_skeleton" in result.columns, "Missing 'geometry_skeleton' column" + assert "geometry_mask" in result.columns, "Missing 'geometry_mask' column" diff --git a/test/vectors/test_run_create_virtual_points.py b/test/vectors/test_run_create_virtual_points.py new file mode 100644 index 00000000..0f9b61ee --- /dev/null +++ b/test/vectors/test_run_create_virtual_points.py @@ -0,0 +1,284 @@ +import os +import shutil +from pathlib import Path + +import geopandas as gpd +import numpy as np +from pyproj import CRS +from shapely.geometry import LineString, Point, Polygon + +from lidro.create_virtual_point.vectors.run_create_virtual_points import ( + launch_virtual_points_by_section, +) + +TMP_PATH = Path("./tmp/create_virtual_point/vectors/run_create_virtual_points/") + + +def setup_module(module): + if TMP_PATH.is_dir(): + shutil.rmtree(TMP_PATH) + os.makedirs(TMP_PATH) + + +def create_test_data_no_points(): + # Create empty GeoDataFrame with the required structure + points = gpd.GeoDataFrame({"geometry": [], "points_knn": []}, crs="EPSG:2154") + + # Create example lines and mask_hydro GeoDataFrames with two entries each + lines = gpd.GeoDataFrame( + { + "geometry": [ + LineString([(700000, 6600000), (700100, 6600100)]), + LineString([(700200, 6600200), (700300, 6600300)]), + ] + }, + crs="EPSG:2154", + ) + + mask_hydro = gpd.GeoDataFrame( + { + "geometry": [ + Polygon( + [(700000, 6600000), (700100, 6600000), (700100, 6600100), (700000, 6600100), (700000, 6600000)] + ), + Polygon( + [(700200, 6600200), (700300, 6600200), (700300, 6600300), (700200, 6600300), (700200, 6600200)] + ), + ] + }, + crs="EPSG:2154", + ) + + return points, lines, mask_hydro + + +def create_test_data_with_points(): + # Create GeoDataFrame with points and points_knn columns + points = gpd.GeoDataFrame( + { + "geometry": [Point(700050, 6600050)], + "points_knn": [ + [ + (700000, 6600000, 10), + (700050, 6600050, 15), + (700100, 6600100, 20), + (700200, 6600200, 30), + (700300, 6600300, 40), + ] + ], + }, + crs="EPSG:2154", + ) + + # Example lines and mask_hydro GeoDataFrames + lines = gpd.GeoDataFrame( + { + "geometry": [ + LineString([(700000, 6600000), (700300, 6600300)]), + LineString([(700400, 6600400), (700500, 6600500)]), + ] + }, + crs="EPSG:2154", + ) + + mask_hydro = gpd.GeoDataFrame( + { + "geometry": [ + Polygon( + [(699900, 6599900), (700400, 6599900), (700400, 6600400), (699900, 6600400), (699900, 6599900)] + ), + Polygon( + [(700400, 6600400), (700500, 6600400), (700500, 6600500), (700400, 6600500), (700400, 6600400)] + ), + ] + }, + crs="EPSG:2154", + ) + + return points, lines, mask_hydro + + +def create_test_data_with_geometry_only(): + # Create GeoDataFrame with points column but without points_knn + points = gpd.GeoDataFrame( + {"geometry": [Point(700050, 6600050), Point(700250, 6600250)], "points_knn": [None, None]}, crs="EPSG:2154" + ) + + # Example line and mask_hydro GeoDataFrames + lines = gpd.GeoDataFrame( + { + "geometry": [ + LineString([(700000, 6600000), (700100, 6600100)]), + LineString([(700200, 6600200), (700300, 6600300)]), + ] + }, + crs="EPSG:2154", + ) + mask_hydro = gpd.GeoDataFrame( + { + "geometry": [ + Polygon( + [(700000, 6600000), (700100, 6600000), (700100, 6600100), (700000, 6600100), (700000, 6600000)] + ), + Polygon( + [(700200, 6600200), (700300, 6600200), (700300, 6600300), (700200, 6600300), (700200, 6600200)] + ), + ] + }, + crs="EPSG:2154", + ) + return points, lines, mask_hydro + + +def create_test_data_with_regression_failure(): + # Create GeoDataFrame with points and points_knn columns that will lead to regression failure + points = gpd.GeoDataFrame( + { + "geometry": [Point(700000, 6600000), Point(700300, 6600300)], + "points_knn": [ + [ + (700000, 6600000, 10), + (700100, 6600100, 10), # Same Z to force failure + (700200, 6600200, 10), # Same Z to force failure + (700300, 6600300, 10), + ], + [ + (700000, 6600000, 10), + (700100, 6600100, 10), # Same Z to force failure + (700200, 6600200, 10), # Same Z to force failure + (700300, 6600300, 10), + ], + ], + }, + crs="EPSG:2154", + ) + + # Create a longer LineString (>150m) to represent a larger river section + lines = gpd.GeoDataFrame( + { + "geometry": [ + LineString([(700000, 6600000), (700200, 6600200)]), # Length is approximately 282.84 meters + ] + }, + crs="EPSG:2154", + ) + + # Create mask_hydro polygons to resemble river segments (long and narrow) + mask_hydro = gpd.GeoDataFrame( + { + "geometry": [ + Polygon( + [(699950, 6599950), (700350, 6599950), (700350, 6600350), (699950, 6600350), (699950, 6599950)] + ), # A polygon resembling a river segment + ] + }, + crs="EPSG:2154", + ) + + return points, lines, mask_hydro + + +def create_test_data_with_flattening_failure(): + # Create GeoDataFrame with points and points_knn columns : the points_knn is None + points = gpd.GeoDataFrame( + {"geometry": [Point(700000, 6600000)], "points_knn": [np.array([])]}, + crs="EPSG:2154", + ) + # Create a short LineString to represent a small river section (<150 meters) + lines = gpd.GeoDataFrame( + { + "geometry": [ + LineString([(700000, 6600000), (700050, 6600050)]), # Length is 70.71 meters + ] + }, + crs="EPSG:2154", + ) + + # Create mask_hydro polygons to resemble river segments (long and narrow) + mask_hydro = gpd.GeoDataFrame( + { + "geometry": [ + Polygon( + [ + (699980, 6599980), + (700020, 6599980), + (700030, 6600000), + (700020, 6600020), + (699980, 6600020), + (699970, 6600000), + ] + ), # A polygon resembling a river bend + ] + }, + crs="EPSG:2154", + ) + + return points, lines, mask_hydro + + +def test_launch_virtual_points_by_section_with_geometry_only(): + points, lines, mask_hydro = create_test_data_with_geometry_only() + crs = CRS.from_epsg(2154) + spacing = 1.0 + river_length = 150 + output_filename = os.path.join(TMP_PATH, "mask_hydro_no_virtual_points.geojson") + + launch_virtual_points_by_section(points, lines, mask_hydro, crs, spacing, river_length, TMP_PATH) + + assert (Path(TMP_PATH) / "mask_hydro_no_virtual_points.geojson").is_file() + masks_without_points = gpd.read_file(output_filename) + assert len(masks_without_points) == len(mask_hydro) + + +def test_launch_virtual_points_by_section_no_points(): + points, line, mask_hydro = create_test_data_no_points() + crs = CRS.from_epsg(2154) + river_length = 150 + spacing = 1.0 + output_filename = os.path.join(TMP_PATH, "mask_hydro_no_virtual_points.geojson") + + launch_virtual_points_by_section(points, line, mask_hydro, crs, spacing, river_length, TMP_PATH) + + assert (Path(TMP_PATH) / "mask_hydro_no_virtual_points.geojson").is_file() + masks_without_points = gpd.read_file(output_filename) + assert len(masks_without_points) == len(mask_hydro) + + +def test_launch_virtual_points_by_section_with_points(): + points, lines, mask_hydro = create_test_data_with_points() + crs = CRS.from_epsg(2154) + spacing = 1.0 + river_length = 150 + + grid_with_z = launch_virtual_points_by_section(points, lines, mask_hydro, crs, spacing, river_length, TMP_PATH) + + assert all(isinstance(geom, Point) for geom in grid_with_z.geometry) + assert all(geom.has_z for geom in grid_with_z.geometry) # Check that all points have a Z coordinate + + +def test_launch_virtual_points_by_section_regression_failure(): + points, lines, mask_hydro = create_test_data_with_regression_failure() + crs = CRS.from_epsg(2154) + spacing = 1.0 + river_length = 150 + output_filename = os.path.join(TMP_PATH, "mask_hydro_no_virtual_points_with_regression.geojson") + + launch_virtual_points_by_section(points, lines, mask_hydro, crs, spacing, river_length, TMP_PATH) + + assert (Path(TMP_PATH) / "mask_hydro_no_virtual_points_with_regression.geojson").is_file() + masks_without_points = gpd.read_file(output_filename) + assert len(masks_without_points) == len(mask_hydro) + + +def test_launch_virtual_points_by_section_flattening_failure(): + points, lines, mask_hydro = create_test_data_with_flattening_failure() + crs = CRS.from_epsg(2154) + spacing = 1.0 + river_length = 150 + output_filename = os.path.join(TMP_PATH, "mask_hydro_no_virtual_points_for_little_river.geojson") + + launch_virtual_points_by_section(points, lines, mask_hydro, crs, spacing, river_length, TMP_PATH) + + assert (Path(TMP_PATH) / "mask_hydro_no_virtual_points_for_little_river.geojson").is_file() + masks_without_points = gpd.read_file(output_filename) + assert len(masks_without_points) == len(mask_hydro)