From 3341e6470afc5f01f682204f9aef0e985015903e Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 30 May 2024 15:37:02 +0200 Subject: [PATCH 01/64] add function : crop lidar pointcloud --- .../pointcloud/crop_las.py | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 lidro/create_virtual_point/pointcloud/crop_las.py diff --git a/lidro/create_virtual_point/pointcloud/crop_las.py b/lidro/create_virtual_point/pointcloud/crop_las.py new file mode 100644 index 00000000..0a1578e1 --- /dev/null +++ b/lidro/create_virtual_point/pointcloud/crop_las.py @@ -0,0 +1,33 @@ +# -*- coding: utf-8 -*- +""" Crop filtered pointcloud """ +from typing import List +import numpy as np +import pdal +from shapely.geometry import MultiPolygon + +def crop_pointcloud_by_points(input_points: str, geom : MultiPolygon, classes: List[int:int])-> np.array: + + """Crop filtered pointcloud : + 1. Filter pointcloud for keeping only classes "1" and "2" + 2. Crop filtered pointcloud by "Mask HYDRO + buffer" + + Args: + input_points (str): Path to the input LAS/LAZ file + geom (MultiPolygon): An array of WKT or GeoJSON 2D MultiPolygon (Mask Hydro with buffer) + classes (list): List of classes to use for the filtering + + Returns: + filtered_points (np.ndarray) : Numpy array containing point coordinates (X, Y, Z) after filtering and croping + """ + # Crop pointcloud by point + pipeline = pdal.Reader.las(filename=input_points, nosrs=True) | pdal.Filter.range( + limits=f"Classification{classes}", + ) | pdal.Filter.crop( + polygon=geom, + ) + pipeline.execute() + # extract points + cropped_points = pipeline.arrays[0] + + return np.array([cropped_points['X'], cropped_points['Y'], cropped_points['Z']]).T + From 0bc5b801b0e9f9fc9278b9b6e5236628f25feae9 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 30 May 2024 15:41:58 +0200 Subject: [PATCH 02/64] add test for function : cropping las --- test/pointcloud/test_crop_las.py | 38 ++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 test/pointcloud/test_crop_las.py diff --git a/test/pointcloud/test_crop_las.py b/test/pointcloud/test_crop_las.py new file mode 100644 index 00000000..043a353c --- /dev/null +++ b/test/pointcloud/test_crop_las.py @@ -0,0 +1,38 @@ +import os +import shutil +from pathlib import Path + +import numpy as np + +from lidro.create_virtual_point.pointcloud.crop_las import crop_pointcloud_by_points + +TMP_PATH = Path("./tmp/create_virtual_point/pointcloud/crop_las") +LAS_FILE = "./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz" +file_mask = "./data/merge_mask_hydro/mask_test2.geojson" + + +def setup_module(module): + if TMP_PATH.is_dir(): + shutil.rmtree(TMP_PATH) + os.makedirs(TMP_PATH) + + +def test_crop_pointcloud_default(): + classes = "[1:2]" + + geom = str( + "MULTIPOLYGON (((830873.1249999998 6290475.625000002, \ + 830610.0480769228 6290335.43269231, \ + 830546.0096153843 6290405.528846156, \ + 830803.0288461535 6290559.567307694, \ + 830873.1249999998 6290475.625000002)), \ + ((830210.5288461539 6290861.298076924, \ + 830205.0480769231 6290760.336538463, \ + 830101.2019230769 6290759.471153848, \ + 830103.2211538461 6290861.586538463, \ + 830210.5288461539 6290861.298076924)))" + ) + output = crop_pointcloud_by_points(LAS_FILE, geom, classes) + + assert isinstance(output, np.ndarray) is True + assert (output == 0).all() is False From 2b4fa04e474d27ae23ab2024af02fb1737f5aaf4 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 30 May 2024 16:06:43 +0200 Subject: [PATCH 03/64] add function who calculates quartile and median from np.array --- .../stats/calculate_stat.py | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 lidro/create_virtual_point/stats/calculate_stat.py diff --git a/lidro/create_virtual_point/stats/calculate_stat.py b/lidro/create_virtual_point/stats/calculate_stat.py new file mode 100644 index 00000000..e9aa7357 --- /dev/null +++ b/lidro/create_virtual_point/stats/calculate_stat.py @@ -0,0 +1,39 @@ +# -*- coding: utf-8 -*- +""" Calculates statistics from np.array (ex. pointcloud) """ +import numpy as np + + +def calculate_quartile(points: np.array, q: int) -> float: + """ + Calculates the quartile of the value (ex. altitude from z-coordinate) of points. + + Parameters: + - points (numpy.ndarray): An array of 3D points [x, y, z] + - q (int): Percentage or sequence of percentages for the percentiles to compute. + Values must be between 0 and 100 inclusive. + + + Returns: + - float: The quartile of the value (ex. z-coordinates). + """ + altitudes = points[:, 2] # Extract the altitude column + n_quartile = round(np.percentile(altitudes, q), 2) + + return n_quartile + + +def calculate_median(points: np.array) -> float: + """ + Calculates the median of the value (ex. altitude from z-coordinate) of points. + + Parameters: + - points (numpy.ndarray): An array of 3D points [x, y, z] + + + Returns: + - float: The median of the value (ex. z-coordinates). + """ + altitudes = points[:, 2] # Extract the altitude column + n_median = round(np.median(altitudes), 2) + + return n_median From 4e0ae4f26e721f1d16a3b6524b1b22289a61c831 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 30 May 2024 16:10:32 +0200 Subject: [PATCH 04/64] add test from calculate_stats.py --- test/stats/test_calculate_stats.py | 51 ++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 test/stats/test_calculate_stats.py diff --git a/test/stats/test_calculate_stats.py b/test/stats/test_calculate_stats.py new file mode 100644 index 00000000..0e0c4c08 --- /dev/null +++ b/test/stats/test_calculate_stats.py @@ -0,0 +1,51 @@ +import numpy as np + +from lidro.create_virtual_point.stats.calculate_stat import ( + calculate_median, + calculate_quartile, +) + + +def test_calculate_quartile_25(): + points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) + assert calculate_quartile(points, 25) == 6.0 + + +def test_calculate_quartile_50(): + points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) + assert calculate_quartile(points, 50) == 9.0 + + +def test_calculate_quartile_75(): + points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) + assert calculate_quartile(points, 75) == 12.0 + + +def test_calculate_quartile_100(): + points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) + assert calculate_quartile(points, 100) == 15.0 + + +def test_calculate_quartile_invalid_q(): + points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) + try: + calculate_quartile(points, 110) + except ValueError: + pass + else: + assert False, "ValueError non levée pour q=110" + + +def test_calculate_quartile_negative_q(): + points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) + try: + calculate_quartile(points, -10) + except ValueError: + pass + else: + assert False, "ValueError non levée pour q=-10" + + +def test_calculate_median(): + points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) + assert calculate_median(points) == 9.0 From 4888491eb2ce08e647b7bec75327b9538d32e1ee Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 30 May 2024 17:28:20 +0200 Subject: [PATCH 05/64] delete paramaters geom --- test/pointcloud/test_crop_las.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/pointcloud/test_crop_las.py b/test/pointcloud/test_crop_las.py index 043a353c..c2db782d 100644 --- a/test/pointcloud/test_crop_las.py +++ b/test/pointcloud/test_crop_las.py @@ -8,7 +8,6 @@ TMP_PATH = Path("./tmp/create_virtual_point/pointcloud/crop_las") LAS_FILE = "./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz" -file_mask = "./data/merge_mask_hydro/mask_test2.geojson" def setup_module(module): From 4f87a7f39ec09e6b70042aac008bb280d81efff7 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 30 May 2024 17:29:14 +0200 Subject: [PATCH 06/64] add function knn from np.array --- .../stats/knn_distance.py | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 lidro/create_virtual_point/stats/knn_distance.py diff --git a/lidro/create_virtual_point/stats/knn_distance.py b/lidro/create_virtual_point/stats/knn_distance.py new file mode 100644 index 00000000..7191167e --- /dev/null +++ b/lidro/create_virtual_point/stats/knn_distance.py @@ -0,0 +1,44 @@ +# -*- coding: utf-8 -*- +""" KNN """ +import numpy as np +from shapely.geometry import Point +from sklearn.neighbors import NearestNeighbors + + +def point_to_numpy(point: Point) -> np.array: + """Converts a 3D shapely.geometry Point to a numpy array + + Args: + point (shapely.geometry.Point): A 3D point with x, and y coordinates + + Returns: + numpy.ndarray: An array of the K nearest 3D point [x, y, z]. + By default [z] is '0' meters + """ + return np.array([point.x, point.y, 0]) + + +def find_k_nearest_neighbors(point: Point, points_array: np.array, k: int) -> np.array: + """Finds the K nearest neighbors of a given point from a list of 3D points + + Args: + point (Point): The reference 2D point + points_array (np.array): An array of 3D points [x, y, z] + k (int): The number of nearest neighbors to find + + Returns: + numpy.ndarray: An array of the K nearest 3D points [x, y, z] + """ + # Convert the point to a numpy array + point_coords = point_to_numpy(point).reshape(1, -1) + + # Initialize and fit the NearestNeighbors model + nbrs = NearestNeighbors(n_neighbors=k).fit(points_array) + + # Find the distances and indices of the nearest neighbors + distances, indices = nbrs.kneighbors(point_coords) + + # Retrieve the points corresponding to the indices + k_nearest_points = points_array[indices[0]] + + return k_nearest_points From faeecc43f5814ddbb4988a0ede833c84e7be5fef Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 30 May 2024 17:29:47 +0200 Subject: [PATCH 07/64] add test for knn --- test/stats/test_knn_distance.py | 34 +++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 test/stats/test_knn_distance.py diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py new file mode 100644 index 00000000..b5e665aa --- /dev/null +++ b/test/stats/test_knn_distance.py @@ -0,0 +1,34 @@ +import numpy as np +from shapely.geometry import Point + +from lidro.create_virtual_point.stats.knn_distance import ( + find_k_nearest_neighbors, + point_to_numpy, +) + + +def test_point_to_numpy_default(): + point = Point(830574.89, 6290648.53, 0) + + result = point_to_numpy(point) + expected = np.array([830574.89, 6290648.53, 0]) + + assert isinstance(result, np.ndarray) is True + assert np.array_equal(result, expected) + + +def test_find_k_nearest_neighbors_default(): + points_array = np.array( + [ + [830438.91, 6290854.32, 2.56], + [830721.84, 6290447.79, 2.23], + [830861.04, 6290242.06, 2.78], + [830867.61, 6290202.62, 2.89], + ] + ) + point = Point(830574.89, 6290648.53) + k = 3 + + result = find_k_nearest_neighbors(point, points_array, k) + + assert isinstance(result, np.ndarray) is True From 3d760fe0ab16cad8a43a707a40202d7b30b81c8d Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 31 May 2024 10:20:13 +0200 Subject: [PATCH 08/64] modify classe : keep only ground --- .../pointcloud/crop_las.py | 22 +++++++++++-------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/lidro/create_virtual_point/pointcloud/crop_las.py b/lidro/create_virtual_point/pointcloud/crop_las.py index 0a1578e1..dfb394b0 100644 --- a/lidro/create_virtual_point/pointcloud/crop_las.py +++ b/lidro/create_virtual_point/pointcloud/crop_las.py @@ -1,14 +1,15 @@ # -*- coding: utf-8 -*- """ Crop filtered pointcloud """ from typing import List + import numpy as np import pdal from shapely.geometry import MultiPolygon -def crop_pointcloud_by_points(input_points: str, geom : MultiPolygon, classes: List[int:int])-> np.array: - + +def crop_pointcloud_by_points(input_points: str, geom: MultiPolygon, classes: List[int:int]) -> np.array: """Crop filtered pointcloud : - 1. Filter pointcloud for keeping only classes "1" and "2" + 1. Filter pointcloud for keeping only classe "GROUND" 2. Crop filtered pointcloud by "Mask HYDRO + buffer" Args: @@ -20,14 +21,17 @@ def crop_pointcloud_by_points(input_points: str, geom : MultiPolygon, classes: filtered_points (np.ndarray) : Numpy array containing point coordinates (X, Y, Z) after filtering and croping """ # Crop pointcloud by point - pipeline = pdal.Reader.las(filename=input_points, nosrs=True) | pdal.Filter.range( - limits=f"Classification{classes}", - ) | pdal.Filter.crop( - polygon=geom, + pipeline = ( + pdal.Reader.las(filename=input_points, nosrs=True) + | pdal.Filter.range( + limits=f"Classification{classes}", + ) + | pdal.Filter.crop( + polygon=geom, ) + ) pipeline.execute() # extract points cropped_points = pipeline.arrays[0] - - return np.array([cropped_points['X'], cropped_points['Y'], cropped_points['Z']]).T + return np.array([cropped_points["X"], cropped_points["Y"], cropped_points["Z"]]).T From 2865602df57b518502e9f8c8450c1263594642c9 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 3 Jun 2024 15:48:55 +0200 Subject: [PATCH 09/64] refacto function calculate_stat.py --- lidro/create_virtual_point/stats/calculate_stat.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/lidro/create_virtual_point/stats/calculate_stat.py b/lidro/create_virtual_point/stats/calculate_stat.py index e9aa7357..6207f73e 100644 --- a/lidro/create_virtual_point/stats/calculate_stat.py +++ b/lidro/create_virtual_point/stats/calculate_stat.py @@ -12,7 +12,6 @@ def calculate_quartile(points: np.array, q: int) -> float: - q (int): Percentage or sequence of percentages for the percentiles to compute. Values must be between 0 and 100 inclusive. - Returns: - float: The quartile of the value (ex. z-coordinates). """ @@ -29,7 +28,6 @@ def calculate_median(points: np.array) -> float: Parameters: - points (numpy.ndarray): An array of 3D points [x, y, z] - Returns: - float: The median of the value (ex. z-coordinates). """ From 6860779d087a158e67d813052ae725367a4efd81 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 3 Jun 2024 15:52:37 +0200 Subject: [PATCH 10/64] refacto function knn and test --- .../stats/knn_distance.py | 32 ++++++------------- test/stats/test_knn_distance.py | 17 ++-------- 2 files changed, 12 insertions(+), 37 deletions(-) diff --git a/lidro/create_virtual_point/stats/knn_distance.py b/lidro/create_virtual_point/stats/knn_distance.py index 7191167e..fe6b3a46 100644 --- a/lidro/create_virtual_point/stats/knn_distance.py +++ b/lidro/create_virtual_point/stats/knn_distance.py @@ -1,44 +1,32 @@ # -*- coding: utf-8 -*- """ KNN """ import numpy as np -from shapely.geometry import Point from sklearn.neighbors import NearestNeighbors -def point_to_numpy(point: Point) -> np.array: - """Converts a 3D shapely.geometry Point to a numpy array - - Args: - point (shapely.geometry.Point): A 3D point with x, and y coordinates - - Returns: - numpy.ndarray: An array of the K nearest 3D point [x, y, z]. - By default [z] is '0' meters - """ - return np.array([point.x, point.y, 0]) - - -def find_k_nearest_neighbors(point: Point, points_array: np.array, k: int) -> np.array: +def find_k_nearest_neighbors(point_2d: np.array, points_3d: np.array, k: int) -> np.array: """Finds the K nearest neighbors of a given point from a list of 3D points Args: - point (Point): The reference 2D point - points_array (np.array): An array of 3D points [x, y, z] + point_2d (np.array): An array of 3D point [x, y, 0] + points_3D (np.array): An array of 3D points [x, y, z] k (int): The number of nearest neighbors to find Returns: numpy.ndarray: An array of the K nearest 3D points [x, y, z] """ - # Convert the point to a numpy array - point_coords = point_to_numpy(point).reshape(1, -1) + # Convert point_2d to nump.arrat + point_2d_array = np.array(point_2d).reshape(1, -1) # Initialize and fit the NearestNeighbors model - nbrs = NearestNeighbors(n_neighbors=k).fit(points_array) + nbrs = NearestNeighbors(n_neighbors=k, algorithm="auto").fit(points_3d) # Find the distances and indices of the nearest neighbors - distances, indices = nbrs.kneighbors(point_coords) + indices = nbrs.kneighbors(point_2d_array, return_distance=False) # Retrieve the points corresponding to the indices - k_nearest_points = points_array[indices[0]] + # Use indices to retrieve the closest 3D points + k_nearest_points = points_3d[indices[0]] + # = points_3d[indices.flatten()] return k_nearest_points diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py index b5e665aa..37a36e00 100644 --- a/test/stats/test_knn_distance.py +++ b/test/stats/test_knn_distance.py @@ -1,20 +1,7 @@ import numpy as np from shapely.geometry import Point -from lidro.create_virtual_point.stats.knn_distance import ( - find_k_nearest_neighbors, - point_to_numpy, -) - - -def test_point_to_numpy_default(): - point = Point(830574.89, 6290648.53, 0) - - result = point_to_numpy(point) - expected = np.array([830574.89, 6290648.53, 0]) - - assert isinstance(result, np.ndarray) is True - assert np.array_equal(result, expected) +from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors def test_find_k_nearest_neighbors_default(): @@ -26,7 +13,7 @@ def test_find_k_nearest_neighbors_default(): [830867.61, 6290202.62, 2.89], ] ) - point = Point(830574.89, 6290648.53) + point = Point(830574.89, 6290648.53, 0) k = 3 result = find_k_nearest_neighbors(point, points_array, k) From b90d069859c92696f8aa49b7964f6f7d0a6f2db1 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 3 Jun 2024 15:57:34 +0200 Subject: [PATCH 11/64] create function and test: calculate Z value along skeleton --- .../vectors/las_around_point.py | 86 +++++++++++++++++++ test/vectors/test_las_around_point.py | 44 ++++++++++ 2 files changed, 130 insertions(+) create mode 100644 lidro/create_virtual_point/vectors/las_around_point.py create mode 100644 test/vectors/test_las_around_point.py diff --git a/lidro/create_virtual_point/vectors/las_around_point.py b/lidro/create_virtual_point/vectors/las_around_point.py new file mode 100644 index 00000000..e343b05b --- /dev/null +++ b/lidro/create_virtual_point/vectors/las_around_point.py @@ -0,0 +1,86 @@ +# -*- coding: utf-8 -*- +""" Extract a Z elevation value every N meters along the hydrographic skeleton +""" +from typing import List + +import geopandas as gpd +import numpy as np +from shapely.geometry import CAP_STYLE, Point + +from lidro.create_virtual_point.pointcloud.crop_las import crop_pointcloud_by_points +from lidro.create_virtual_point.stats.calculate_stat import calculate_quartile +from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors +from lidro.create_virtual_point.vectors.points_along_skeleton import ( + generate_points_along_skeleton, +) + + +def apply_buffers_to_geometry(hydro_mask, buffer: float): + """Buffer geometry + Objective: create a HYDRO mask largest + + Args: + hydro_mask (gpd.GeoDataFrame): geopandas dataframe with input geometry + buffer (float): buffer to apply to the input geometry + Returns: + GeoJSON: updated geometry + """ + # Buffer + geom = hydro_mask.buffer(buffer, cap_style=CAP_STYLE.square) + return geom + + +def filter_las_around_point( + file_skeleton: str, + file_mask: str, + file_lidar: str, + distance_meters: float, + k: int, + classes: List[int], + crs: str | int, +) -> List: + """Extract a Z elevation value every 2 meters along the hydrographic skeleton + + Args: + file_skeleton (str): Path from Skeleton + file_mask (str): Path from Mask Hydro + file_lidar (str): Path to the input LAS/LAZ file + distance_meters (float): distance in meters betwen each point + k (int): The number of nearest neighbors to find + classes (List[int]): List of classes to use for the binarisation (points with other + classification values are ignored) + crs (str | int): Make a CRS from an EPSG code : CRS WKT string + + Returns: + List : Result {'geometry': Point 3D, 'z_q1': points KNN} + """ + # Read Mask Hydro merged + gdf = gpd.read_file(file_mask, crs=crs).unary_union + + # Apply buffer (2 meters) from Mask Hydro + gdf = apply_buffers_to_geometry(gdf, 2) + + # Create severals points every 2 meters (by default) along skeleton Hydro + # Return GeoDataframe + points_gdf = generate_points_along_skeleton(file_skeleton, distance_meters, crs) + + points_list = [([geom.x, geom.y, 0]) for geom in points_gdf.geometry if isinstance(geom, Point)] + + # Filter pointcloud by classes: keep only points from selected classe(s) ("2" : DEFAULT) + # Crop filtered pointcloud by Mask HYDRO who have a buffer + points_clip = crop_pointcloud_by_points(file_lidar, str(gdf.wkt), classes) + + # Finds the K nearest neighbors of a given point from a list of 3D points + points_knn_list = [ + ({"geometry": Point(point), "points_knn": find_k_nearest_neighbors(point, points_clip, k)}) + for point in points_list + ] + + # Calcule Z "Q1" for each points every 2 meters (by default) along skeleton hydro + results = [ + ({"geometry": p["geometry"], "z_q1": calculate_quartile(p["points_knn"], 10)}) + for p in points_knn_list + if not np.all(np.isnan(p["points_knn"])) + ] + + return results diff --git a/test/vectors/test_las_around_point.py b/test/vectors/test_las_around_point.py new file mode 100644 index 00000000..4d66ebdc --- /dev/null +++ b/test/vectors/test_las_around_point.py @@ -0,0 +1,44 @@ +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.las_around_point import filter_las_around_point + +TMP_PATH = Path("./tmp/create_virtual_point/vectors/las_around_point") + +file_skeleton = "./data/skeleton_hydro/Skeleton_Hydro.geojson" +file_mask = "./data/merge_mask_hydro/MaskHydro_merge.geojson" +file_lidar = "./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz" +output = Path("./tmp/create_virtual_point/vectors/las_around_point/Points.geojson") + + +def setup_module(module): + if TMP_PATH.is_dir(): + shutil.rmtree(TMP_PATH) + os.makedirs(TMP_PATH) + + +def test_las_around_point_default(): + # Parameters + distance_meters = 10 + k = 6 + classes = "[2:2]" + crs = CRS.from_epsg(2154) + + result = filter_las_around_point(file_skeleton, file_mask, file_lidar, distance_meters, k, classes, crs) + + # 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") + + 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 From 365db60f5089e550a5bc2ae4c7de97dbb3c5592d Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 3 Jun 2024 16:53:42 +0200 Subject: [PATCH 12/64] rectify test: (output == 0).all() == false --- test/pointcloud/test_crop_las.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/pointcloud/test_crop_las.py b/test/pointcloud/test_crop_las.py index c2db782d..b3269f81 100644 --- a/test/pointcloud/test_crop_las.py +++ b/test/pointcloud/test_crop_las.py @@ -34,4 +34,4 @@ def test_crop_pointcloud_default(): output = crop_pointcloud_by_points(LAS_FILE, geom, classes) assert isinstance(output, np.ndarray) is True - assert (output == 0).all() is False + assert isinstance((output == 0).all(), bool) is False From 52c554186e4f5e5a0403aa452c4905949fa11161 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 3 Jun 2024 16:54:18 +0200 Subject: [PATCH 13/64] rectify test knn distance --- test/stats/test_knn_distance.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py index 37a36e00..8f943eb6 100644 --- a/test/stats/test_knn_distance.py +++ b/test/stats/test_knn_distance.py @@ -1,5 +1,4 @@ import numpy as np -from shapely.geometry import Point from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors @@ -13,7 +12,7 @@ def test_find_k_nearest_neighbors_default(): [830867.61, 6290202.62, 2.89], ] ) - point = Point(830574.89, 6290648.53, 0) + point = np.array([[830574.89, 6290648.53, 0]]) k = 3 result = find_k_nearest_neighbors(point, points_array, k) From ab9ad1c3c3098662f0e417739e0457ebf0f17e09 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 10:39:43 +0200 Subject: [PATCH 14/64] refacto configs with function create virtual point --- configs/configs_lidro.yaml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index 117e4114..8ea6d30a 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -8,6 +8,8 @@ work_dir: ${hydra:runtime.cwd} io: input_filename: null + input_mask_hydro: null + input_skeleton: null input_dir: null output_dir: null srid: 2154 @@ -31,10 +33,15 @@ vector: tolerance: 1 # 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 filter: # Classes to be considered as "non-water" keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67] # All classes + keep_classe_ground: "[2:2]" hydra: output_subdir: null From cc5c2cdc134426e6c93565b984065adf3e77af44 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 10:40:36 +0200 Subject: [PATCH 15/64] refacto function crop_las --- lidro/create_virtual_point/pointcloud/crop_las.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lidro/create_virtual_point/pointcloud/crop_las.py b/lidro/create_virtual_point/pointcloud/crop_las.py index dfb394b0..ba1b017a 100644 --- a/lidro/create_virtual_point/pointcloud/crop_las.py +++ b/lidro/create_virtual_point/pointcloud/crop_las.py @@ -18,7 +18,7 @@ def crop_pointcloud_by_points(input_points: str, geom: MultiPolygon, classes: Li classes (list): List of classes to use for the filtering Returns: - filtered_points (np.ndarray) : Numpy array containing point coordinates (X, Y, Z) after filtering and croping + np.array : Numpy array containing point coordinates (X, Y, Z) after filtering and croping """ # Crop pointcloud by point pipeline = ( From dbadb2d8ef6f8fcdce49694d9d8bd5ff6a92a917 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 10:52:20 +0200 Subject: [PATCH 16/64] refacto script for lauching steps --- example_create_virtual_point_by_tile.sh | 10 ++++++++++ example_create_virtual_point_default.sh | 9 +++++++++ example_merge_mask_default.sh | 2 +- 3 files changed, 20 insertions(+), 1 deletion(-) create mode 100644 example_create_virtual_point_by_tile.sh create mode 100644 example_create_virtual_point_default.sh diff --git a/example_create_virtual_point_by_tile.sh b/example_create_virtual_point_by_tile.sh new file mode 100644 index 00000000..6a50c135 --- /dev/null +++ b/example_create_virtual_point_by_tile.sh @@ -0,0 +1,10 @@ +# For lauching create virtual points +python -m lidro.main_create_virtual_point \ +io.input_dir=./data/ \ +io.input_filename=Semis_2021_0830_6291_LA93_IGN69.laz \ +io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geosjon \ +io.input_skeleton=./data/skeleton_hydro/skeleton_hydro.geojson \ +io.output_dir=./tmp/ \ + + + diff --git a/example_create_virtual_point_default.sh b/example_create_virtual_point_default.sh new file mode 100644 index 00000000..c184bcf7 --- /dev/null +++ b/example_create_virtual_point_default.sh @@ -0,0 +1,9 @@ +# For lauching Create virtual point +python -m lidro.main_create_virtual_point \ +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/ \ + + + diff --git a/example_merge_mask_default.sh b/example_merge_mask_default.sh index a2d6e65f..7e162fc7 100644 --- a/example_merge_mask_default.sh +++ b/example_merge_mask_default.sh @@ -1,4 +1,4 @@ -# For lauching Mask Hydro +# For lauching Mask Hydro merged python -m lidro.main_merge_mask \ io.input_dir=./data/mask_hydro/ \ io.output_dir=./tmp/merge_mask_hydro/ \ From 8819752a3ccb658039a46a34e46146f13de933ef Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 10:53:47 +0200 Subject: [PATCH 17/64] refacto script las_around_point for lauching --- .../vectors/las_around_point.py | 43 +++---------------- 1 file changed, 5 insertions(+), 38 deletions(-) diff --git a/lidro/create_virtual_point/vectors/las_around_point.py b/lidro/create_virtual_point/vectors/las_around_point.py index e343b05b..3eb6a22c 100644 --- a/lidro/create_virtual_point/vectors/las_around_point.py +++ b/lidro/create_virtual_point/vectors/las_around_point.py @@ -7,15 +7,11 @@ import numpy as np from shapely.geometry import CAP_STYLE, Point -from lidro.create_virtual_point.pointcloud.crop_las import crop_pointcloud_by_points from lidro.create_virtual_point.stats.calculate_stat import calculate_quartile from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors -from lidro.create_virtual_point.vectors.points_along_skeleton import ( - generate_points_along_skeleton, -) -def apply_buffers_to_geometry(hydro_mask, buffer: float): +def apply_buffers_to_geometry(hydro_mask: gpd.GeoDataFrame, buffer: float): """Buffer geometry Objective: create a HYDRO mask largest @@ -30,50 +26,21 @@ def apply_buffers_to_geometry(hydro_mask, buffer: float): return geom -def filter_las_around_point( - file_skeleton: str, - file_mask: str, - file_lidar: str, - distance_meters: float, - k: int, - classes: List[int], - crs: str | int, -) -> List: +def filter_las_around_point(points_skeleton: List, points_clip: np.array, k: int) -> List: """Extract a Z elevation value every 2 meters along the hydrographic skeleton Args: - file_skeleton (str): Path from Skeleton - file_mask (str): Path from Mask Hydro - file_lidar (str): Path to the input LAS/LAZ file - distance_meters (float): distance in meters betwen each point + points_skeleton (list) : points every 2 meters (by default) along skeleton Hydro + points_clip (np.array): Numpy array containing point coordinates (X, Y, Z) after filtering and croping k (int): The number of nearest neighbors to find - classes (List[int]): List of classes to use for the binarisation (points with other - classification values are ignored) - crs (str | int): Make a CRS from an EPSG code : CRS WKT string Returns: List : Result {'geometry': Point 3D, 'z_q1': points KNN} """ - # Read Mask Hydro merged - gdf = gpd.read_file(file_mask, crs=crs).unary_union - - # Apply buffer (2 meters) from Mask Hydro - gdf = apply_buffers_to_geometry(gdf, 2) - - # Create severals points every 2 meters (by default) along skeleton Hydro - # Return GeoDataframe - points_gdf = generate_points_along_skeleton(file_skeleton, distance_meters, crs) - - points_list = [([geom.x, geom.y, 0]) for geom in points_gdf.geometry if isinstance(geom, Point)] - - # Filter pointcloud by classes: keep only points from selected classe(s) ("2" : DEFAULT) - # Crop filtered pointcloud by Mask HYDRO who have a buffer - points_clip = crop_pointcloud_by_points(file_lidar, str(gdf.wkt), classes) - # Finds the K nearest neighbors of a given point from a list of 3D points points_knn_list = [ ({"geometry": Point(point), "points_knn": find_k_nearest_neighbors(point, points_clip, k)}) - for point in points_list + for point in points_skeleton ] # Calcule Z "Q1" for each points every 2 meters (by default) along skeleton hydro From 06f8e6bf6f4a7da778892927e3df793c1f24b291 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 10:54:40 +0200 Subject: [PATCH 18/64] refacto test for las around point --- test/vectors/test_las_around_point.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/test/vectors/test_las_around_point.py b/test/vectors/test_las_around_point.py index 4d66ebdc..083b0ea2 100644 --- a/test/vectors/test_las_around_point.py +++ b/test/vectors/test_las_around_point.py @@ -3,6 +3,7 @@ from pathlib import Path import geopandas as gpd +import numpy as np from pyproj import CRS from shapely.geometry import Point @@ -10,9 +11,6 @@ TMP_PATH = Path("./tmp/create_virtual_point/vectors/las_around_point") -file_skeleton = "./data/skeleton_hydro/Skeleton_Hydro.geojson" -file_mask = "./data/merge_mask_hydro/MaskHydro_merge.geojson" -file_lidar = "./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz" output = Path("./tmp/create_virtual_point/vectors/las_around_point/Points.geojson") @@ -24,12 +22,23 @@ def setup_module(module): def test_las_around_point_default(): # Parameters - distance_meters = 10 + points_along_skeleton = [[830864.5181373736, 6290217.943739296, 0], [830866.5957826116, 6290208.162525126, 0]] + + points_clip = np.array( + [ + [8.30822700e05, 6.29000133e06, 2.59000000e00], + [8.30836950e05, 6.29000254e06, 2.47000000e00], + [8.30837730e05, 6.29000379e06, 2.58000000e00], + [8.30042740e05, 6.29041476e06, 3.09000000e00], + [8.30033610e05, 6.29041739e06, 3.87000000e00], + [8.30042180e05, 6.29041447e06, 3.10000000e00], + ] + ) + k = 6 - classes = "[2:2]" crs = CRS.from_epsg(2154) - result = filter_las_around_point(file_skeleton, file_mask, file_lidar, distance_meters, k, classes, crs) + result = filter_las_around_point(points_along_skeleton, points_clip, k) # Convert results to GeoDataFrame result_gdf = gpd.GeoDataFrame(result) From 0be63c9164290c9580abb43f3218c6e48991f486 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 10:55:57 +0200 Subject: [PATCH 19/64] add new function create mask hydro buffer --- .../vectors/create_mask_hydro_buffer.py | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 lidro/create_virtual_point/vectors/create_mask_hydro_buffer.py diff --git a/lidro/create_virtual_point/vectors/create_mask_hydro_buffer.py b/lidro/create_virtual_point/vectors/create_mask_hydro_buffer.py new file mode 100644 index 00000000..fb7089e6 --- /dev/null +++ b/lidro/create_virtual_point/vectors/create_mask_hydro_buffer.py @@ -0,0 +1,41 @@ +# -*- coding: utf-8 -*- +""" Extract a Z elevation value every N meters along the hydrographic skeleton +""" +import geopandas as gpd +from shapely.geometry import CAP_STYLE, MultiPolygon + + +def apply_buffers_to_geometry(hydro_mask: gpd.GeoDataFrame, buffer: float) -> MultiPolygon: + """Buffer geometry + Objective: create a HYDRO mask largest + + Args: + hydro_mask (gpd.GeoDataFrame): geopandas dataframe with input geometry + buffer (float): buffer to apply to the input geometry + + Returns: + MutliPolygon: Mask Hydro with buffer + """ + # Buffer + geom = hydro_mask.buffer(buffer, cap_style=CAP_STYLE.square) + return geom + + +def create_mask_hydro_buffer(file_mask: str, buffer: float, crs: str | int) -> gpd.GeoDataFrame: + """Apply buffer (2 meters by default) from Mask Hydro + + Args: + file_mask (str): Path from Mask Hydro + buffer (float): buffer to apply to the input geometry + crs (str | int): Make a CRS from an EPSG code : CRS WKT string + + Returns: + gpd.GeoDataFrame : geometry columns are encoded to WKT + """ + # Read Mask Hydro merged + gdf = gpd.read_file(file_mask, crs=crs).unary_union + + # Apply buffer (2 meters by default) from Mask Hydro + gdf = apply_buffers_to_geometry(gdf, buffer) + + return gdf From b01963d9422e96c208005e512ec313afb601be25 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 10:57:59 +0200 Subject: [PATCH 20/64] add function for lauching create virtual point without rectify Z --- lidro/main_create_virtual_point.py | 114 +++++++++++++++++++++++++++++ 1 file changed, 114 insertions(+) create mode 100644 lidro/main_create_virtual_point.py diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py new file mode 100644 index 00000000..cd3702c5 --- /dev/null +++ b/lidro/main_create_virtual_point.py @@ -0,0 +1,114 @@ +""" Main script for calculate Mask HYDRO 1 +""" + +import logging +import os + +import geopandas as gpd +import hydra +import numpy as np +from omegaconf import DictConfig +from pyproj import CRS +from shapely.geometry import Point + +from lidro.create_virtual_point.pointcloud.crop_las import crop_pointcloud_by_points +from lidro.create_virtual_point.vectors.create_mask_hydro_buffer import ( + create_mask_hydro_buffer, +) +from lidro.create_virtual_point.vectors.las_around_point import filter_las_around_point +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 a virtual point along hydro surfaces 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 + + 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.vector.distance_meters + buffer = config.vector.buffer + crs = CRS.from_user_input(config.io.srid) + classes = config.filter.keep_classe_ground + k = config.vector.k + + # Step 1 : Create Mask Hydro with buffer + # Return GeoDataframe + input_mask_hydro_buffer = create_mask_hydro_buffer(input_mask_hydro, buffer, crs).wkt + + # Step 2 : Create several points every 2 meters (by default) along skeleton Hydro + # Return GeoDataframe + points_gdf = generate_points_along_skeleton(input_skeleton, distance_meters, crs) + + points_list_skeleton = [([geom.x, geom.y, 0]) for geom in points_gdf.geometry if isinstance(geom, Point)] + + # Step 3 : Crope filtered pointcloud by Mask Hydro with buffer + def main_on_one_tile(filename): + """Lauch main.py on one tile + + Args: + filename (str): filename to the LAS file + + Returns: + points_clip (np.array) : Numpy array containing point coordinates (X, Y, Z) after filtering and croping + """ + input_dir_points = os.path.join(input_dir, "pointcloud") + tilename = os.path.splitext(filename)[0] # filename to the LAS file + input_pointcloud = os.path.join(input_dir_points, filename) # path to the LAS file + logging.info(f"\nCroping filtered pointcloud by Mask Hydro with buffer for tile : {tilename}") + points_clip = crop_pointcloud_by_points(input_pointcloud, str(input_mask_hydro_buffer), classes) + return points_clip + + if initial_las_filename: + # Lauch croping filtered pointcloud by Mask Hydro with buffer by one tile: + points_clip = main_on_one_tile(initial_las_filename) + + 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 = [main_on_one_tile(file) for file in os.listdir(input_dir_points)] + points_clip = np.vstack(points_clip_list) # merge ONE numpy.array + + # Step 4 : Extract a Z elevation value every 2 meters along the hydrographic skeleton + result = filter_las_around_point(points_list_skeleton, points_clip, k) + + # Convert results to GeoDataFrame + result_gdf = gpd.GeoDataFrame(result) + 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") + + +if __name__ == "__main__": + main() From 12652d2b6955ef9f536c16b46571ddafc047449b Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 10:58:39 +0200 Subject: [PATCH 21/64] add test for point along skeleton --- test/vectors/test_points_along_skeleton.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/vectors/test_points_along_skeleton.py b/test/vectors/test_points_along_skeleton.py index 97c3329c..19e09be7 100644 --- a/test/vectors/test_points_along_skeleton.py +++ b/test/vectors/test_points_along_skeleton.py @@ -11,10 +11,10 @@ generate_points_along_skeleton, ) -TMP_PATH = Path("./tmp/create_point_virtual/vectors/points_along_skeleton") +TMP_PATH = Path("./tmp/create_virtual_point/vectors/points_along_skeleton") file = "./data/skeleton_hydro/Skeleton_Hydro.geojson" -output = Path("./tmp/create_point_virtual/vectors/points_along_skeleton/Points.geojson") +output = Path("./tmp/create_virtual_point/vectors/points_along_skeleton/Points.geojson") def setup_module(module): From f8d0880a3e8fc8646bd7c5e8f82d6b4dbac13967 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 10:59:18 +0200 Subject: [PATCH 22/64] add test for lauching virtual point without rectify Z --- test/test_main_create_virtual_point.py | 217 +++++++++++++++++++++++++ 1 file changed, 217 insertions(+) create mode 100644 test/test_main_create_virtual_point.py diff --git a/test/test_main_create_virtual_point.py b/test/test_main_create_virtual_point.py new file mode 100644 index 00000000..51415c61 --- /dev/null +++ b/test/test_main_create_virtual_point.py @@ -0,0 +1,217 @@ +import os +import subprocess as sp +from pathlib import Path + +import pytest +from hydra import compose, initialize + +from lidro.main_create_virtual_point import main + +INPUT_DIR = Path("data") +OUTPUT_DIR = Path("tmp") / "create_virtual_point/main" + + +def setup_module(module): + os.makedirs("tmp/create_virtual_point/main", exist_ok=True) + + +def test_main_run_okay(): + repo_dir = Path.cwd().parent + cmd = f"""python -m lidro.main_create_virtual_point \ + 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.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"vector.distance_meters={distances_meters}", + f"vector.buffer={buffer}", + f"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 + 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"vector.distance_meters={distances_meters}", + f"vector.buffer={buffer}", + f"vector.k={k}", + f"io.srid={srid}", + ], + ) + 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"vector.distance_meters={distances_meters}", + f"vector.buffer={buffer}", + f"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"vector.distance_meters={distances_meters}", + f"vector.buffer={buffer}", + f"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"vector.distance_meters={distances_meters}", + f"vector.buffer={buffer}", + f"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"vector.distance_meters={distances_meters}", + f"vector.buffer={buffer}", + f"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"vector.distance_meters={distances_meters}", + f"vector.buffer={buffer}", + f"vector.k={k}", + f"io.srid={srid}", + ], + ) + with pytest.raises(ValueError): + main(cfg) From b79e75e520859d1c67ac7cd8acb41c377e6702f2 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 6 Jun 2024 11:00:28 +0200 Subject: [PATCH 23/64] add test for creating mask hydro buffer --- test/vectors/test_create_mask_hydro_buffer.py | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) create mode 100644 test/vectors/test_create_mask_hydro_buffer.py diff --git a/test/vectors/test_create_mask_hydro_buffer.py b/test/vectors/test_create_mask_hydro_buffer.py new file mode 100644 index 00000000..d8f9b57d --- /dev/null +++ b/test/vectors/test_create_mask_hydro_buffer.py @@ -0,0 +1,44 @@ +import os +import shutil +from pathlib import Path + +import geopandas as gpd +from pyproj import CRS +from shapely.geometry import MultiPolygon + +from lidro.create_virtual_point.vectors.create_mask_hydro_buffer import ( + create_mask_hydro_buffer, +) + +TMP_PATH = Path("./tmp/create_virtual_point/vectors/create_mask_hydro_buffer/") + +file_mask = "./data/merge_mask_hydro/MaskHydro_merge.geojson" +output = Path("./tmp/create_virtual_point/vectors/create_mask_hydro_buffer/MaskHydro_merge_buffer.geojson") + + +def setup_module(module): + if TMP_PATH.is_dir(): + shutil.rmtree(TMP_PATH) + os.makedirs(TMP_PATH) + + +def test_create_mask_hydro_buffer_default(): + # Parameters + buffer = 2 + crs = CRS.from_epsg(2154) + + result = create_mask_hydro_buffer(file_mask, buffer, crs) + + # Save the result to GeoJSON + d = {"nom": "mask_hydro_merge_buffer", "geometry": [result]} + gdf = gpd.GeoDataFrame(d, crs=crs) + gdf.to_file(output, driver="GeoJSON", crs=crs) + + assert Path(output).exists() + + 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, MultiPolygon) for geom in gdf.geometry) # All geometry should MultiPolygone From d65f4c9db4747f9e321d963e0ec30bf4f34ac63f Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 7 Jun 2024 15:58:12 +0200 Subject: [PATCH 24/64] add function identify upstream/downstream --- .../vectors/identify_upstream_downstream.py | 53 +++++++++++++++++++ .../test_identify_upstream_downstream.py | 43 +++++++++++++++ 2 files changed, 96 insertions(+) create mode 100644 lidro/create_virtual_point/vectors/identify_upstream_downstream.py create mode 100644 test/vectors/test_identify_upstream_downstream.py diff --git a/lidro/create_virtual_point/vectors/identify_upstream_downstream.py b/lidro/create_virtual_point/vectors/identify_upstream_downstream.py new file mode 100644 index 00000000..d97d33d4 --- /dev/null +++ b/lidro/create_virtual_point/vectors/identify_upstream_downstream.py @@ -0,0 +1,53 @@ +# -*- coding: utf-8 -*- +""" Identify upstream and downstream by skeleton's section +""" +from typing import Dict, List, Tuple + + +def calculate_z_q1_averages(points: List) -> Tuple: + """Calculate the average 'Z_q1' values for the first 10 percent and the last 10 percent of a list of points + + Parameters: + points (list): A list of dictionaries where each dictionary contains a 'geometry' key + with a Point object and a 'z_q1' key with a numeric value. + + Returns: + tuple: A pair of values representing the average 'Z_q1' values + for the first 10 percent and the last 10 percent of the list, respectively. + """ + total_points = len(points) + percent_10 = max(1, total_points // 10) # Ensure at least 1 point to avoid division by zero + + # Extract the first 10 percent and the last 10 percent of points + first_10_percent_points = points[:percent_10] + last_10_percent_points = points[-percent_10:] + + # Calculate the average 'Z_q1' values for each group + average_z_q1_first_10_percent = sum(point["z_q1"] for point in first_10_percent_points) / percent_10 + average_z_q1_last_10_percent = sum(point["z_q1"] for point in last_10_percent_points) / percent_10 + + return average_z_q1_first_10_percent, average_z_q1_last_10_percent + + +def identify_upstream_downstream(points: List) -> Dict: + """Determinate point for upstream and downstream by skeleton's section + + Args: + points (List): A list of dictionaries where each dictionary contains a 'geometry' key + with a Point object and a 'z_q1' key with a numeric value. + + Returns: + Dict : A dictionnary representing geometry (Point 2D) and value of "z_q1" + for upstream and downstream by skeleton's section + """ + # Calculate the average 'Z_q1' values for the first 10 percent and the last 10 percent of a list of points + average_z_q1_first_10_percent, average_z_q1_last_10_percent = calculate_z_q1_averages(points) + + if average_z_q1_first_10_percent < average_z_q1_last_10_percent: + upstream_point = points[0] + downstream_point = points[-1] + else: + upstream_point = points[-1] + downstream_point = points[0] + + return upstream_point, downstream_point diff --git a/test/vectors/test_identify_upstream_downstream.py b/test/vectors/test_identify_upstream_downstream.py new file mode 100644 index 00000000..c6216efd --- /dev/null +++ b/test/vectors/test_identify_upstream_downstream.py @@ -0,0 +1,43 @@ +from typing import Dict + +from shapely.geometry import Point + +from lidro.create_virtual_point.vectors.identify_upstream_downstream import ( + identify_upstream_downstream, +) + + +def test_identify_upstream_downstream_default(): + # Example list of points + points = [ + {"geometry": Point(830211.190425319829956, 6292000.43919328507036), "z_q1": 2.95}, + {"geometry": Point(830211.950744876405224, 6291980.453650656156242), "z_q1": 3.05}, + {"geometry": Point(830212.711064432864077, 6291960.468108027242124), "z_q1": 3.20}, + {"geometry": Point(830213.471383989439346, 6291940.482565398328006), "z_q1": 3.15}, + {"geometry": Point(830214.231703546014614, 6291920.497022769413888), "z_q1": 3.22}, + {"geometry": Point(830214.992023102473468, 6291900.511480140499771), "z_q1": 4.12}, + {"geometry": Point(830218.033301328658126, 6291820.569309624843299), "z_q1": 4.12}, + {"geometry": Point(830218.793620885233395, 6291800.583766995929182), "z_q1": 4.55}, + {"geometry": Point(830219.553940441692248, 6291780.598224367015064), "z_q1": 4.54}, + {"geometry": Point(830224.663624010980129, 6291680.817003472708166), "z_q1": 5.00}, + {"geometry": Point(830228.343047266826034, 6291661.158370648510754), "z_q1": 5.05}, + {"geometry": Point(830232.022470522555523, 6291641.499737825244665), "z_q1": 5.05}, + {"geometry": Point(830246.740163545822725, 6291562.865206529386342), "z_q1": 5.10}, + {"geometry": Point(830250.419586801668629, 6291543.20657370518893), "z_q1": 5.15}, + {"geometry": Point(830254.099010057398118, 6291523.547940881922841), "z_q1": 2.95}, + {"geometry": Point(830257.778433313244022, 6291503.88930805772543), "z_q1": 2.96}, + {"geometry": Point(830261.457856569089927, 6291484.230675233528018), "z_q1": 6.00}, + {"geometry": Point(830265.137279824819416, 6291464.572042410261929), "z_q1": 6.10}, + {"geometry": Point(830268.81670308066532, 6291444.913409586064517), "z_q1": 6.15}, + {"geometry": Point(830272.496126336511225, 6291425.254776761867106), "z_q1": 6.10}, + {"geometry": Point(830287.621702362317592, 6291346.699770421721041), "z_q1": 7.00}, + ] + + upstream, downstream = identify_upstream_downstream(points) + + assert isinstance(upstream, Dict) + assert isinstance(downstream, Dict) + assert upstream["geometry"] == Point(830211.1904253198, 6292000.439193285) + assert upstream["z_q1"] == 2.95 + assert downstream["geometry"] == Point(830287.6217023623, 6291346.699770422) + assert downstream["z_q1"] == 7.0 From eb89910601664000b47e414688e740d682ea061d Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 7 Jun 2024 16:00:43 +0200 Subject: [PATCH 25/64] add function sorted geospatial points by skeleton --- .../vectors/sort_geospatial_points.py | 34 ++++++++++ test/vectors/test_sort_geospatial_points.py | 64 +++++++++++++++++++ 2 files changed, 98 insertions(+) create mode 100644 lidro/create_virtual_point/vectors/sort_geospatial_points.py create mode 100644 test/vectors/test_sort_geospatial_points.py diff --git a/lidro/create_virtual_point/vectors/sort_geospatial_points.py b/lidro/create_virtual_point/vectors/sort_geospatial_points.py new file mode 100644 index 00000000..6f97cf92 --- /dev/null +++ b/lidro/create_virtual_point/vectors/sort_geospatial_points.py @@ -0,0 +1,34 @@ +# -*- coding: utf-8 -*- +""" Sorts a list of 2D points based on their distance from the first point in the list +""" +from typing import List + +from shapely.geometry import Point + + +def sort_geospatial_points(points: List, reference_point: Point) -> List: + """Sorts a list of 2D points in EPSG:2154 (metric) based on their distance + from the first point in the list (= upstream and downstream by skeleton's section) + + Parameters: + points (list of dictionnary): A list of tuples where each tuple represents + a point (x, y) in metric coordinates. + reference_point (Point): upstream by skeleton's section + + Returns: + list of tuple: The sorted list of points based on their increasing distance from the first point. + """ + + def euclidean_distance(point1, point2): + """Calculates the Euclidean distance between two points.""" + x1, y1 = point1.x, point1.y + x2, y2 = point2.x, point2.y + return ((x2 - x1) ** 2 + (y2 - y1) ** 2) ** 0.5 + + if not points: + return [] + + # Sort the points by distance to the reference point + sorted_points = sorted(points, key=lambda point: euclidean_distance(reference_point, point["geometry"])) + + return sorted_points diff --git a/test/vectors/test_sort_geospatial_points.py b/test/vectors/test_sort_geospatial_points.py new file mode 100644 index 00000000..60f9225d --- /dev/null +++ b/test/vectors/test_sort_geospatial_points.py @@ -0,0 +1,64 @@ +from typing import List + +from shapely.geometry import Point + +from lidro.create_virtual_point.vectors.sort_geospatial_points import ( + sort_geospatial_points, +) + + +def test_sort_geospatial_points_default(): + # Example list of points + points = [ + {"geometry": Point(830211.190425319829956, 6292000.43919328507036), "z_q1": 2.95}, + {"geometry": Point(830211.950744876405224, 6291980.453650656156242), "z_q1": 3.05}, + {"geometry": Point(830214.992023102473468, 6291900.511480140499771), "z_q1": 4.12}, + {"geometry": Point(830212.711064432864077, 6291960.468108027242124), "z_q1": 3.20}, + {"geometry": Point(830213.471383989439346, 6291940.482565398328006), "z_q1": 3.15}, + {"geometry": Point(830214.231703546014614, 6291920.497022769413888), "z_q1": 3.22}, + {"geometry": Point(830218.033301328658126, 6291820.569309624843299), "z_q1": 4.12}, + {"geometry": Point(830218.793620885233395, 6291800.583766995929182), "z_q1": 4.55}, + {"geometry": Point(830219.553940441692248, 6291780.598224367015064), "z_q1": 4.54}, + {"geometry": Point(830224.663624010980129, 6291680.817003472708166), "z_q1": 5.00}, + {"geometry": Point(830228.343047266826034, 6291661.158370648510754), "z_q1": 5.05}, + {"geometry": Point(830257.778433313244022, 6291503.88930805772543), "z_q1": 2.96}, + {"geometry": Point(830232.022470522555523, 6291641.499737825244665), "z_q1": 5.05}, + {"geometry": Point(830246.740163545822725, 6291562.865206529386342), "z_q1": 5.10}, + {"geometry": Point(830250.419586801668629, 6291543.20657370518893), "z_q1": 5.15}, + {"geometry": Point(830254.099010057398118, 6291523.547940881922841), "z_q1": 2.95}, + {"geometry": Point(830261.457856569089927, 6291484.230675233528018), "z_q1": 6.00}, + {"geometry": Point(830265.137279824819416, 6291464.572042410261929), "z_q1": 6.10}, + {"geometry": Point(830268.81670308066532, 6291444.913409586064517), "z_q1": 6.15}, + {"geometry": Point(830272.496126336511225, 6291425.254776761867106), "z_q1": 6.10}, + {"geometry": Point(830287.621702362317592, 6291346.699770421721041), "z_q1": 7.00}, + ] + + reference_point = Point(830211.190425319829956, 6292000.43919328507036) + + result = sort_geospatial_points(points, reference_point) + + assert isinstance(result, List) + assert points[0]["geometry"] == Point(830211.190425319829956, 6292000.43919328507036) + assert result == [ + {"geometry": Point(830211.190425319829956, 6292000.43919328507036), "z_q1": 2.95}, + {"geometry": Point(830211.950744876405224, 6291980.453650656156242), "z_q1": 3.05}, + {"geometry": Point(830212.711064432864077, 6291960.468108027242124), "z_q1": 3.20}, + {"geometry": Point(830213.471383989439346, 6291940.482565398328006), "z_q1": 3.15}, + {"geometry": Point(830214.231703546014614, 6291920.497022769413888), "z_q1": 3.22}, + {"geometry": Point(830214.992023102473468, 6291900.511480140499771), "z_q1": 4.12}, + {"geometry": Point(830218.033301328658126, 6291820.569309624843299), "z_q1": 4.12}, + {"geometry": Point(830218.793620885233395, 6291800.583766995929182), "z_q1": 4.55}, + {"geometry": Point(830219.553940441692248, 6291780.598224367015064), "z_q1": 4.54}, + {"geometry": Point(830224.663624010980129, 6291680.817003472708166), "z_q1": 5.00}, + {"geometry": Point(830228.343047266826034, 6291661.158370648510754), "z_q1": 5.05}, + {"geometry": Point(830232.022470522555523, 6291641.499737825244665), "z_q1": 5.05}, + {"geometry": Point(830246.740163545822725, 6291562.865206529386342), "z_q1": 5.10}, + {"geometry": Point(830250.419586801668629, 6291543.20657370518893), "z_q1": 5.15}, + {"geometry": Point(830254.099010057398118, 6291523.547940881922841), "z_q1": 2.95}, + {"geometry": Point(830257.778433313244022, 6291503.88930805772543), "z_q1": 2.96}, + {"geometry": Point(830261.457856569089927, 6291484.230675233528018), "z_q1": 6.00}, + {"geometry": Point(830265.137279824819416, 6291464.572042410261929), "z_q1": 6.10}, + {"geometry": Point(830268.81670308066532, 6291444.913409586064517), "z_q1": 6.15}, + {"geometry": Point(830272.496126336511225, 6291425.254776761867106), "z_q1": 6.10}, + {"geometry": Point(830287.621702362317592, 6291346.699770421721041), "z_q1": 7.00}, + ] From a78d9d301f36cd2902e4d200466f43294e23185d Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 13 Jun 2024 14:58:10 +0200 Subject: [PATCH 26/64] refacto config.yaml with under-categories --- configs/configs_lidro.yaml | 48 ++++++++++++++------------ lidro/main_create_mask.py | 4 +-- lidro/main_create_virtual_point.py | 8 ++--- lidro/main_merge_mask.py | 8 ++--- test/test_main_create_virtual_point.py | 42 +++++++++++----------- test/test_main_merge_mask.py | 24 ++++++------- 6 files changed, 69 insertions(+), 65 deletions(-) diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index 8ea6d30a..13660ac2 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -19,29 +19,33 @@ io: no_data_value: -9999 tile_size: 1000 -raster: - # size for dilatation - dilation_size: 3 +mask_generation: + raster: + # size for dilatation + dilation_size: 3 + filter: + # Classes to be considered as "non-water" + keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67] # All classes + vector: + # Filter water's area (m²) + min_water_area: 150 + # Parameters for buffer + buffer_positive: 1 + buffer_negative: -1.5 # negative buffer should be bigger than positive buffer to prevent protruding over the banks + # Tolerance from Douglas-Peucker + tolerance: 1 -vector: - # Filter water's area (m²) - min_water_area: 150 - # Parameters for buffer - buffer_positive: 1 - buffer_negative: -1.5 # negative buffer should be bigger than positive buffer to prevent protruding over the banks - # Tolerance from Douglas-Peucker - tolerance: 1 - # 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 - -filter: - # Classes to be considered as "non-water" - keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67] # All classes - keep_classe_ground: "[2:2]" +virtual_point: + filter: + # Keep ground pointcloud between Hydro Mask and Hydro Mask buffer + keep_classe_ground: "[2:2]" + 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 hydra: output_subdir: null diff --git a/lidro/main_create_mask.py b/lidro/main_create_mask.py index 3e01fb1c..6c4927f5 100644 --- a/lidro/main_create_mask.py +++ b/lidro/main_create_mask.py @@ -45,8 +45,8 @@ def main(config: DictConfig): pixel_size = config.io.pixel_size tile_size = config.io.tile_size crs = CRS.from_user_input(config.io.srid) - classe = config.filter.keep_classes - dilation_size = config.raster.dilation_size + classe = config.mask_generation.filter.keep_classes + dilation_size = config.mask_generation.raster.dilation_size def main_on_one_tile(filename): """Lauch main.py on one tile diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index cd3702c5..2f1b6baa 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -55,11 +55,11 @@ def main(config: DictConfig): input_skeleton = config.io.input_skeleton # Parameters for creating virtual point - distance_meters = config.vector.distance_meters - buffer = config.vector.buffer + distance_meters = config.virtual_point.vector.distance_meters + buffer = config.virtual_point.vector.buffer crs = CRS.from_user_input(config.io.srid) - classes = config.filter.keep_classe_ground - k = config.vector.k + classes = config.virtual_point.filter.keep_classe_ground + k = config.virtual_point.vector.k # Step 1 : Create Mask Hydro with buffer # Return GeoDataframe diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py index 350625de..49854413 100644 --- a/lidro/main_merge_mask.py +++ b/lidro/main_merge_mask.py @@ -43,10 +43,10 @@ def main(config: DictConfig): input_dir, output_dir, CRS.from_user_input(config.io.srid), - config.vector.min_water_area, - config.vector.buffer_positive, - config.vector.buffer_negative, - config.vector.tolerance) + config.mask_generation.vector.min_water_area, + config.mask_generation.vector.buffer_positive, + config.mask_generation.vector.buffer_negative, + config.mask_generation.vector.tolerance) if __name__ == "__main__": diff --git a/test/test_main_create_virtual_point.py b/test/test_main_create_virtual_point.py index 51415c61..27401864 100644 --- a/test/test_main_create_virtual_point.py +++ b/test/test_main_create_virtual_point.py @@ -46,9 +46,9 @@ def test_main_lidro_default(): f"io.output_dir={output_dir}", f"io.input_mask_hydro={input_mask_hydro}", f"io.input_skeleton={input_skeleton}", - f"vector.distance_meters={distances_meters}", - f"vector.buffer={buffer}", - f"vector.k={k}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", f"io.srid={srid}", ], ) @@ -76,9 +76,9 @@ 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"vector.distance_meters={distances_meters}", - f"vector.buffer={buffer}", - f"vector.k={k}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", f"io.srid={srid}", ], ) @@ -102,9 +102,9 @@ def test_main_lidro_fail_no_input_dir(): f"io.input_mask_hydro={input_mask_hydro}", f"io.input_skeleton={input_skeleton}", f"io.output_dir={output_dir}", - f"vector.distance_meters={distances_meters}", - f"vector.buffer={buffer}", - f"vector.k={k}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", f"io.srid={srid}", ], ) @@ -129,9 +129,9 @@ def test_main_lidro_fail_wrong_input_dir(): f"io.input_mask_hydro={input_mask_hydro}", f"io.input_skeleton={input_skeleton}", f"io.output_dir={output_dir}", - f"vector.distance_meters={distances_meters}", - f"vector.buffer={buffer}", - f"vector.k={k}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", f"io.srid={srid}", ], ) @@ -155,9 +155,9 @@ def test_main_lidro_fail_no_output_dir(): f"io.input_dir={input_dir}", f"io.input_mask_hydro={input_mask_hydro}", f"io.input_skeleton={input_skeleton}", - f"vector.distance_meters={distances_meters}", - f"vector.buffer={buffer}", - f"vector.k={k}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", f"io.srid={srid}", ], ) @@ -181,9 +181,9 @@ def test_main_lidro_fail_no_input_mask_hydro(): f"io.output_dir={output_dir}", "io.input_mask_hydro=some_random_input_mask_hydro", f"io.input_skeleton={input_skeleton}", - f"vector.distance_meters={distances_meters}", - f"vector.buffer={buffer}", - f"vector.k={k}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", f"io.srid={srid}", ], ) @@ -207,9 +207,9 @@ def test_main_lidro_fail_no_input_skeleton(): f"io.output_dir={output_dir}", f"io.input_mask_hydro={input_mask_hydro}", "io.input_skeleton=some_random_input_skeleton", - f"vector.distance_meters={distances_meters}", - f"vector.buffer={buffer}", - f"vector.k={k}", + f"virtual_point.vector.distance_meters={distances_meters}", + f"virtual_point.vector.buffer={buffer}", + f"virtual_point.vector.k={k}", f"io.srid={srid}", ], ) diff --git a/test/test_main_merge_mask.py b/test/test_main_merge_mask.py index a4d8beb3..e7750a09 100644 --- a/test/test_main_merge_mask.py +++ b/test/test_main_merge_mask.py @@ -39,10 +39,10 @@ def test_main_lidro_fail_no_input_dir(): overrides=[ f"io.output_dir={output_dir}", f"io.pixel_size={pixel_size}", - f"vector.min_water_area={min_water_area}", - f"vector.buffer_positive={buffer_positive}", - f"vector.buffer_negative={buffer_negative}", - f"vector.tolerance={tolerance}", + f"mask_generation.vector.min_water_area={min_water_area}", + f"mask_generation.vector.buffer_positive={buffer_positive}", + f"mask_generation.vector.buffer_negative={buffer_negative}", + f"mask_generation.vector.tolerance={tolerance}", f"io.srid={srid}", ], ) @@ -66,10 +66,10 @@ def test_main_lidro_fail_wrong_input_dir(): "io.input_dir=some_random_input_dir", f"io.output_dir={output_dir}", f"io.pixel_size={pixel_size}", - f"vector.min_water_area={min_water_area}", - f"vector.buffer_positive={buffer_positive}", - f"vector.buffer_negative={buffer_negative}", - f"vector.tolerance={tolerance}", + f"mask_generation.vector.min_water_area={min_water_area}", + f"mask_generation.vector.buffer_positive={buffer_positive}", + f"mask_generation.vector.buffer_negative={buffer_negative}", + f"mask_generation.vector.tolerance={tolerance}", f"io.srid={srid}", ], ) @@ -92,10 +92,10 @@ def test_main_lidro_fail_no_output_dir(): overrides=[ f"io.input_dir={input_dir}", f"io.pixel_size={pixel_size}", - f"vector.min_water_area={min_water_area}", - f"vector.buffer_positive={buffer_positive}", - f"vector.buffer_negative={buffer_negative}", - f"vector.tolerance={tolerance}", + f"mask_generation.vector.min_water_area={min_water_area}", + f"mask_generation.vector.buffer_positive={buffer_positive}", + f"mask_generation.vector.buffer_negative={buffer_negative}", + f"mask_generation.vector.tolerance={tolerance}", f"io.srid={srid}", ], ) From 47e9d25271e49a09b3bbbd0de91e2bcfa9f5888c Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 13 Jun 2024 15:00:37 +0200 Subject: [PATCH 27/64] refacto gitignore: delete all logs --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 08b71877..d6262489 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,3 @@ **/__pycache__ tmp -main.log \ No newline at end of file +**.log \ No newline at end of file From 592eefc26944fb926f06327c486754f317edc14e Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 13 Jun 2024 15:08:57 +0200 Subject: [PATCH 28/64] refacto function step 1: import mask hydro and apply buffer --- .../vectors/create_mask_hydro_buffer.py | 41 ------------------- .../vectors/mask_hydro_with_buffer.py | 25 +++++++++++ lidro/main_create_virtual_point.py | 10 ++--- test/vectors/test_create_mask_hydro_buffer.py | 6 +-- 4 files changed, 33 insertions(+), 49 deletions(-) delete mode 100644 lidro/create_virtual_point/vectors/create_mask_hydro_buffer.py create mode 100644 lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py diff --git a/lidro/create_virtual_point/vectors/create_mask_hydro_buffer.py b/lidro/create_virtual_point/vectors/create_mask_hydro_buffer.py deleted file mode 100644 index fb7089e6..00000000 --- a/lidro/create_virtual_point/vectors/create_mask_hydro_buffer.py +++ /dev/null @@ -1,41 +0,0 @@ -# -*- coding: utf-8 -*- -""" Extract a Z elevation value every N meters along the hydrographic skeleton -""" -import geopandas as gpd -from shapely.geometry import CAP_STYLE, MultiPolygon - - -def apply_buffers_to_geometry(hydro_mask: gpd.GeoDataFrame, buffer: float) -> MultiPolygon: - """Buffer geometry - Objective: create a HYDRO mask largest - - Args: - hydro_mask (gpd.GeoDataFrame): geopandas dataframe with input geometry - buffer (float): buffer to apply to the input geometry - - Returns: - MutliPolygon: Mask Hydro with buffer - """ - # Buffer - geom = hydro_mask.buffer(buffer, cap_style=CAP_STYLE.square) - return geom - - -def create_mask_hydro_buffer(file_mask: str, buffer: float, crs: str | int) -> gpd.GeoDataFrame: - """Apply buffer (2 meters by default) from Mask Hydro - - Args: - file_mask (str): Path from Mask Hydro - buffer (float): buffer to apply to the input geometry - crs (str | int): Make a CRS from an EPSG code : CRS WKT string - - Returns: - gpd.GeoDataFrame : geometry columns are encoded to WKT - """ - # Read Mask Hydro merged - gdf = gpd.read_file(file_mask, crs=crs).unary_union - - # Apply buffer (2 meters by default) from Mask Hydro - gdf = apply_buffers_to_geometry(gdf, buffer) - - return gdf diff --git a/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py b/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py new file mode 100644 index 00000000..531baa8c --- /dev/null +++ b/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py @@ -0,0 +1,25 @@ +# -*- coding: utf-8 -*- +""" Extract a Z elevation value every N meters along the hydrographic skeleton +""" +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 + + Args: + file_mask (str): Path from Mask Hydro + buffer (float): buffer to apply to the input geometry + crs (str | int): Make a CRS from an EPSG code : CRS WKT string + + Returns: + gpd.GeoDataFrame : geometry columns are encoded to WKT + """ + # 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) + + return gdf_buffer diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 2f1b6baa..2fada611 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -12,10 +12,10 @@ from shapely.geometry import Point from lidro.create_virtual_point.pointcloud.crop_las import crop_pointcloud_by_points -from lidro.create_virtual_point.vectors.create_mask_hydro_buffer import ( - create_mask_hydro_buffer, -) from lidro.create_virtual_point.vectors.las_around_point import filter_las_around_point +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, ) @@ -61,9 +61,9 @@ def main(config: DictConfig): classes = config.virtual_point.filter.keep_classe_ground k = config.virtual_point.vector.k - # Step 1 : Create Mask Hydro with buffer + # Step 1 : Import Mask Hydro, then apply a buffer # Return GeoDataframe - input_mask_hydro_buffer = create_mask_hydro_buffer(input_mask_hydro, buffer, crs).wkt + 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 diff --git a/test/vectors/test_create_mask_hydro_buffer.py b/test/vectors/test_create_mask_hydro_buffer.py index d8f9b57d..6591425a 100644 --- a/test/vectors/test_create_mask_hydro_buffer.py +++ b/test/vectors/test_create_mask_hydro_buffer.py @@ -6,8 +6,8 @@ from pyproj import CRS from shapely.geometry import MultiPolygon -from lidro.create_virtual_point.vectors.create_mask_hydro_buffer import ( - create_mask_hydro_buffer, +from lidro.create_virtual_point.vectors.mask_hydro_with_buffer import ( + import_mask_hydro_with_buffer, ) TMP_PATH = Path("./tmp/create_virtual_point/vectors/create_mask_hydro_buffer/") @@ -27,7 +27,7 @@ def test_create_mask_hydro_buffer_default(): buffer = 2 crs = CRS.from_epsg(2154) - result = create_mask_hydro_buffer(file_mask, buffer, crs) + result = import_mask_hydro_with_buffer(file_mask, buffer, crs) # Save the result to GeoJSON d = {"nom": "mask_hydro_merge_buffer", "geometry": [result]} From 98c0f21aae9cd6472645566bb813e69e4f50b4a9 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 13 Jun 2024 15:18:13 +0200 Subject: [PATCH 29/64] delete fonction apply_buffer: not only use --- .../vectors/las_around_point.py | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/lidro/create_virtual_point/vectors/las_around_point.py b/lidro/create_virtual_point/vectors/las_around_point.py index 3eb6a22c..09f8af52 100644 --- a/lidro/create_virtual_point/vectors/las_around_point.py +++ b/lidro/create_virtual_point/vectors/las_around_point.py @@ -3,29 +3,13 @@ """ from typing import List -import geopandas as gpd import numpy as np -from shapely.geometry import CAP_STYLE, Point +from shapely.geometry import Point from lidro.create_virtual_point.stats.calculate_stat import calculate_quartile from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors -def apply_buffers_to_geometry(hydro_mask: gpd.GeoDataFrame, buffer: float): - """Buffer geometry - Objective: create a HYDRO mask largest - - Args: - hydro_mask (gpd.GeoDataFrame): geopandas dataframe with input geometry - buffer (float): buffer to apply to the input geometry - Returns: - GeoJSON: updated geometry - """ - # Buffer - geom = hydro_mask.buffer(buffer, cap_style=CAP_STYLE.square) - return geom - - def filter_las_around_point(points_skeleton: List, points_clip: np.array, k: int) -> List: """Extract a Z elevation value every 2 meters along the hydrographic skeleton From 915c4e4ac4fc107504411922e1ea4b6bb0164506 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 13 Jun 2024 15:22:07 +0200 Subject: [PATCH 30/64] delete fucntion sort_geospatial_points --- .../vectors/sort_geospatial_points.py | 34 ---------- test/vectors/test_sort_geospatial_points.py | 64 ------------------- 2 files changed, 98 deletions(-) delete mode 100644 lidro/create_virtual_point/vectors/sort_geospatial_points.py delete mode 100644 test/vectors/test_sort_geospatial_points.py diff --git a/lidro/create_virtual_point/vectors/sort_geospatial_points.py b/lidro/create_virtual_point/vectors/sort_geospatial_points.py deleted file mode 100644 index 6f97cf92..00000000 --- a/lidro/create_virtual_point/vectors/sort_geospatial_points.py +++ /dev/null @@ -1,34 +0,0 @@ -# -*- coding: utf-8 -*- -""" Sorts a list of 2D points based on their distance from the first point in the list -""" -from typing import List - -from shapely.geometry import Point - - -def sort_geospatial_points(points: List, reference_point: Point) -> List: - """Sorts a list of 2D points in EPSG:2154 (metric) based on their distance - from the first point in the list (= upstream and downstream by skeleton's section) - - Parameters: - points (list of dictionnary): A list of tuples where each tuple represents - a point (x, y) in metric coordinates. - reference_point (Point): upstream by skeleton's section - - Returns: - list of tuple: The sorted list of points based on their increasing distance from the first point. - """ - - def euclidean_distance(point1, point2): - """Calculates the Euclidean distance between two points.""" - x1, y1 = point1.x, point1.y - x2, y2 = point2.x, point2.y - return ((x2 - x1) ** 2 + (y2 - y1) ** 2) ** 0.5 - - if not points: - return [] - - # Sort the points by distance to the reference point - sorted_points = sorted(points, key=lambda point: euclidean_distance(reference_point, point["geometry"])) - - return sorted_points diff --git a/test/vectors/test_sort_geospatial_points.py b/test/vectors/test_sort_geospatial_points.py deleted file mode 100644 index 60f9225d..00000000 --- a/test/vectors/test_sort_geospatial_points.py +++ /dev/null @@ -1,64 +0,0 @@ -from typing import List - -from shapely.geometry import Point - -from lidro.create_virtual_point.vectors.sort_geospatial_points import ( - sort_geospatial_points, -) - - -def test_sort_geospatial_points_default(): - # Example list of points - points = [ - {"geometry": Point(830211.190425319829956, 6292000.43919328507036), "z_q1": 2.95}, - {"geometry": Point(830211.950744876405224, 6291980.453650656156242), "z_q1": 3.05}, - {"geometry": Point(830214.992023102473468, 6291900.511480140499771), "z_q1": 4.12}, - {"geometry": Point(830212.711064432864077, 6291960.468108027242124), "z_q1": 3.20}, - {"geometry": Point(830213.471383989439346, 6291940.482565398328006), "z_q1": 3.15}, - {"geometry": Point(830214.231703546014614, 6291920.497022769413888), "z_q1": 3.22}, - {"geometry": Point(830218.033301328658126, 6291820.569309624843299), "z_q1": 4.12}, - {"geometry": Point(830218.793620885233395, 6291800.583766995929182), "z_q1": 4.55}, - {"geometry": Point(830219.553940441692248, 6291780.598224367015064), "z_q1": 4.54}, - {"geometry": Point(830224.663624010980129, 6291680.817003472708166), "z_q1": 5.00}, - {"geometry": Point(830228.343047266826034, 6291661.158370648510754), "z_q1": 5.05}, - {"geometry": Point(830257.778433313244022, 6291503.88930805772543), "z_q1": 2.96}, - {"geometry": Point(830232.022470522555523, 6291641.499737825244665), "z_q1": 5.05}, - {"geometry": Point(830246.740163545822725, 6291562.865206529386342), "z_q1": 5.10}, - {"geometry": Point(830250.419586801668629, 6291543.20657370518893), "z_q1": 5.15}, - {"geometry": Point(830254.099010057398118, 6291523.547940881922841), "z_q1": 2.95}, - {"geometry": Point(830261.457856569089927, 6291484.230675233528018), "z_q1": 6.00}, - {"geometry": Point(830265.137279824819416, 6291464.572042410261929), "z_q1": 6.10}, - {"geometry": Point(830268.81670308066532, 6291444.913409586064517), "z_q1": 6.15}, - {"geometry": Point(830272.496126336511225, 6291425.254776761867106), "z_q1": 6.10}, - {"geometry": Point(830287.621702362317592, 6291346.699770421721041), "z_q1": 7.00}, - ] - - reference_point = Point(830211.190425319829956, 6292000.43919328507036) - - result = sort_geospatial_points(points, reference_point) - - assert isinstance(result, List) - assert points[0]["geometry"] == Point(830211.190425319829956, 6292000.43919328507036) - assert result == [ - {"geometry": Point(830211.190425319829956, 6292000.43919328507036), "z_q1": 2.95}, - {"geometry": Point(830211.950744876405224, 6291980.453650656156242), "z_q1": 3.05}, - {"geometry": Point(830212.711064432864077, 6291960.468108027242124), "z_q1": 3.20}, - {"geometry": Point(830213.471383989439346, 6291940.482565398328006), "z_q1": 3.15}, - {"geometry": Point(830214.231703546014614, 6291920.497022769413888), "z_q1": 3.22}, - {"geometry": Point(830214.992023102473468, 6291900.511480140499771), "z_q1": 4.12}, - {"geometry": Point(830218.033301328658126, 6291820.569309624843299), "z_q1": 4.12}, - {"geometry": Point(830218.793620885233395, 6291800.583766995929182), "z_q1": 4.55}, - {"geometry": Point(830219.553940441692248, 6291780.598224367015064), "z_q1": 4.54}, - {"geometry": Point(830224.663624010980129, 6291680.817003472708166), "z_q1": 5.00}, - {"geometry": Point(830228.343047266826034, 6291661.158370648510754), "z_q1": 5.05}, - {"geometry": Point(830232.022470522555523, 6291641.499737825244665), "z_q1": 5.05}, - {"geometry": Point(830246.740163545822725, 6291562.865206529386342), "z_q1": 5.10}, - {"geometry": Point(830250.419586801668629, 6291543.20657370518893), "z_q1": 5.15}, - {"geometry": Point(830254.099010057398118, 6291523.547940881922841), "z_q1": 2.95}, - {"geometry": Point(830257.778433313244022, 6291503.88930805772543), "z_q1": 2.96}, - {"geometry": Point(830261.457856569089927, 6291484.230675233528018), "z_q1": 6.00}, - {"geometry": Point(830265.137279824819416, 6291464.572042410261929), "z_q1": 6.10}, - {"geometry": Point(830268.81670308066532, 6291444.913409586064517), "z_q1": 6.15}, - {"geometry": Point(830272.496126336511225, 6291425.254776761867106), "z_q1": 6.10}, - {"geometry": Point(830287.621702362317592, 6291346.699770421721041), "z_q1": 7.00}, - ] From 82402f249e2075ac6cf7e2461d24e58f02c7e957 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 13 Jun 2024 16:01:26 +0200 Subject: [PATCH 31/64] rename function main_one_tile to extract_points_around_skeleton_points_one_tile --- lidro/main_create_virtual_point.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 2fada611..6740252c 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -72,7 +72,7 @@ def main(config: DictConfig): points_list_skeleton = [([geom.x, geom.y, 0]) for geom in points_gdf.geometry if isinstance(geom, Point)] # Step 3 : Crope filtered pointcloud by Mask Hydro with buffer - def main_on_one_tile(filename): + def extract_points_around_skeleton_points_one_tile(filename): """Lauch main.py on one tile Args: @@ -90,12 +90,14 @@ def main_on_one_tile(filename): if initial_las_filename: # Lauch croping filtered pointcloud by Mask Hydro with buffer by one tile: - points_clip = main_on_one_tile(initial_las_filename) + points_clip = extract_points_around_skeleton_points_one_tile(initial_las_filename) 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 = [main_on_one_tile(file) for file in os.listdir(input_dir_points)] + points_clip_list = [ + extract_points_around_skeleton_points_one_tile(file) for file in os.listdir(input_dir_points) + ] points_clip = np.vstack(points_clip_list) # merge ONE numpy.array # Step 4 : Extract a Z elevation value every 2 meters along the hydrographic skeleton From 574ca4906cebbfe5553f4b022d3383989288df16 Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 14 Jun 2024 15:53:40 +0200 Subject: [PATCH 32/64] add function for extract bounding box with pdal --- .../pointcloud/utils_pdal.py | 35 +++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 lidro/create_virtual_point/pointcloud/utils_pdal.py diff --git a/lidro/create_virtual_point/pointcloud/utils_pdal.py b/lidro/create_virtual_point/pointcloud/utils_pdal.py new file mode 100644 index 00000000..670ee636 --- /dev/null +++ b/lidro/create_virtual_point/pointcloud/utils_pdal.py @@ -0,0 +1,35 @@ +# -*- coding: utf-8 -*- +""" Utils PDAL +""" +import pdal + + +def read_las_file(input_las: str): + """Read a las file and put it in an array""" + pipeline = pdal.Pipeline() | pdal.Reader.las(filename=input_las, nosrs=True) + pipeline.execute() + return pipeline.arrays[0] + + +def get_info_from_las(points): + """Get info from a las to put it in an array""" + pipeline = pdal.Filter.stats().pipeline(points) + pipeline.execute() + return pipeline.metadata + + +def get_bounds_from_las(in_las: str): + """Get bounds=([minx,maxx],[miny,maxy]) from las file + + Args: + in_las (str): Path to the input LAS/LAZ file + + Returns: + tuple : Bounds = ([minx,maxx],[miny,maxy]) from las file + """ + metadata = get_info_from_las(read_las_file(in_las)) + xmin = metadata["metadata"]["filters.stats"]["bbox"]["native"]["bbox"]["minx"] + xmax = metadata["metadata"]["filters.stats"]["bbox"]["native"]["bbox"]["maxx"] + ymin = metadata["metadata"]["filters.stats"]["bbox"]["native"]["bbox"]["miny"] + ymax = metadata["metadata"]["filters.stats"]["bbox"]["native"]["bbox"]["maxy"] + return ([xmin, xmax], [ymin, ymax]) From 2f922d4dfabb6aa512da92b4ef70c767ddb3c7dc Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 14 Jun 2024 15:54:11 +0200 Subject: [PATCH 33/64] add test for utils pdal --- test/pointcloud/test_utils_pdal.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 test/pointcloud/test_utils_pdal.py diff --git a/test/pointcloud/test_utils_pdal.py b/test/pointcloud/test_utils_pdal.py new file mode 100644 index 00000000..7fb8bce8 --- /dev/null +++ b/test/pointcloud/test_utils_pdal.py @@ -0,0 +1,23 @@ +import numpy as np + +from lidro.create_virtual_point.pointcloud.utils_pdal import ( + get_bounds_from_las, + read_las_file, +) + +LAS_FILE = "./data/pointcloud/LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST.las" + +# Expected : ([minx,maxx],[miny,maxy]) +BoundsExpected_FILE = ([706044.94, 706292.25], [6626159, 6626324]) + + +def test_read_las_file(): + pts = read_las_file(LAS_FILE) + assert isinstance(pts, np.ndarray) # type is array + + +def test_get_bounds_from_las(): + bounds = get_bounds_from_las(LAS_FILE) + + assert isinstance(bounds, tuple) + assert bounds == BoundsExpected_FILE From 930d731d0ad7d14a3ace9cb70f0b0b0b9ad2cbf5 Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 14 Jun 2024 15:54:57 +0200 Subject: [PATCH 34/64] refacto fucntion las_around_point: N meter default --- lidro/create_virtual_point/vectors/las_around_point.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lidro/create_virtual_point/vectors/las_around_point.py b/lidro/create_virtual_point/vectors/las_around_point.py index 09f8af52..a0be698a 100644 --- a/lidro/create_virtual_point/vectors/las_around_point.py +++ b/lidro/create_virtual_point/vectors/las_around_point.py @@ -11,10 +11,10 @@ def filter_las_around_point(points_skeleton: List, points_clip: np.array, k: int) -> List: - """Extract a Z elevation value every 2 meters along the hydrographic skeleton + """Extract a Z elevation value every N meters along the hydrographic skeleton Args: - points_skeleton (list) : points every 2 meters (by default) along skeleton Hydro + points_skeleton (list) : points every N meters (by default: 2) along skeleton Hydro points_clip (np.array): Numpy array containing point coordinates (X, Y, Z) after filtering and croping k (int): The number of nearest neighbors to find @@ -27,7 +27,7 @@ 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 2 meters (by default) along skeleton hydro + # Calcule Z "Q1" for each points every N meters along skeleton hydro results = [ ({"geometry": p["geometry"], "z_q1": calculate_quartile(p["points_knn"], 10)}) for p in points_knn_list From d1a0cfa21b56d04918c856b9b575b093409b5059 Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 14 Jun 2024 15:59:04 +0200 Subject: [PATCH 35/64] add fucntion: clip points with bounding box --- .../vectors/clip_points_with_bounding_box.py | 34 +++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py diff --git a/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py b/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py new file mode 100644 index 00000000..37e4cd10 --- /dev/null +++ b/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py @@ -0,0 +1,34 @@ +# -*- coding: utf-8 -*- +""" Clip Skeleton points by tile (pointcloud) +""" +from typing import List + +import geopandas as gpd +from shapely.geometry import Point, box + + +def clip_points_with_box(points: List, bbox: tuple) -> gpd.GeoDataFrame: + """Clip skeleton points by tile (bounding box) + + Args: + points (List): Points every 2 meters (by default) along skeleton hydro + bbox (tuple): bounding box from tile (pointcloud) + + Returns: + gpd.GeoDataframe : Points every 2 meters (by default) along skeleton hydro by tile + """ + # Extract the bounding box limits + xmin = bbox[0][0] + xmax = bbox[0][1] + ymin = bbox[1][0] + ymax = bbox[1][1] + + # Create a GeoDataFrame from the points + gdf_points = gpd.GeoDataFrame(geometry=[Point(point) for point in points]) + # Create a polygon representing the bounding box + bounding_box = box(xmin, ymin, xmax, ymax) + + # Clip points using the bounding box + clipped_points = gdf_points[gdf_points.within(bounding_box)] + + return clipped_points From 6cec9c841a4936820e58800e7d65d79cd165c08d Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 14 Jun 2024 16:27:23 +0200 Subject: [PATCH 36/64] add function: clip point with bounding box --- .../test_clip_points_with_bounding_box.py | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 test/vectors/test_clip_points_with_bounding_box.py diff --git a/test/vectors/test_clip_points_with_bounding_box.py b/test/vectors/test_clip_points_with_bounding_box.py new file mode 100644 index 00000000..fa916787 --- /dev/null +++ b/test/vectors/test_clip_points_with_bounding_box.py @@ -0,0 +1,38 @@ +import geopandas as gpd +from shapely.geometry import Point + +from lidro.create_virtual_point.vectors.clip_points_with_bounding_box import ( + clip_points_with_box, +) + + +def test_create_hydro_vector_mask_default(): + # Define test points (both inside and outside the bounding box) + points = [ + (830378.719, 6290999.42), # Inside + (830382.382, 6290990.115), # Inside + (830386.045, 6290980.81), # Inside + (830389.709, 6290971.505), # Inside + (830216.893, 6291850.548), # Outside + (830217.273, 6291840.555), # Outside + (830375.000, 6291001.000), # Outside + (830400.000, 6291000.000), # Outside + ] + + # Define a bounding box that should clip out some points + bbox = ((830375, 830390), (6290970, 6291000)) + + # Expected points within the bounding box + expected_points = [ + (830378.719, 6290999.42), + (830382.382, 6290990.115), + (830386.045, 6290980.81), + (830389.709, 6290971.505), + ] + + # Convert expected points to GeoDataFrame for comparison + expected_gdf = gpd.GeoDataFrame(geometry=[Point(point) for point in expected_points]) + + result_gdf = clip_points_with_box(points, bbox) + + assert (result_gdf.equals(expected_gdf)) is True From 0b4e98df21e0c8fe5c659d897489fa125bc7aec8 Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 14 Jun 2024 16:28:28 +0200 Subject: [PATCH 37/64] refacto function: main create virtual point --- lidro/main_create_virtual_point.py | 34 ++++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 6740252c..4e2fb48e 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -6,12 +6,16 @@ import geopandas as gpd import hydra -import numpy as np +import pandas as pd from omegaconf import DictConfig from pyproj import CRS from shapely.geometry import Point from lidro.create_virtual_point.pointcloud.crop_las import crop_pointcloud_by_points +from lidro.create_virtual_point.pointcloud.utils_pdal import get_bounds_from_las +from lidro.create_virtual_point.vectors.clip_points_with_bounding_box import ( + clip_points_with_box, +) from lidro.create_virtual_point.vectors.las_around_point import filter_las_around_point from lidro.create_virtual_point.vectors.mask_hydro_with_buffer import ( import_mask_hydro_with_buffer, @@ -68,7 +72,6 @@ def main(config: DictConfig): # Step 2 : Create several points every 2 meters (by default) along skeleton Hydro # Return GeoDataframe points_gdf = generate_points_along_skeleton(input_skeleton, distance_meters, crs) - points_list_skeleton = [([geom.x, geom.y, 0]) for geom in points_gdf.geometry if isinstance(geom, Point)] # Step 3 : Crope filtered pointcloud by Mask Hydro with buffer @@ -85,8 +88,21 @@ def extract_points_around_skeleton_points_one_tile(filename): tilename = os.path.splitext(filename)[0] # filename to the LAS file input_pointcloud = os.path.join(input_dir_points, filename) # path to the LAS file logging.info(f"\nCroping filtered pointcloud by Mask Hydro with buffer for tile : {tilename}") + # Croping filtered pointcloud by Mask Hydro with buffer for tile points_clip = crop_pointcloud_by_points(input_pointcloud, str(input_mask_hydro_buffer), classes) - return points_clip + logging.info(f"\nCropping skeleton points for tile: {tilename}") + # Extract bounding box for clipping points by tile + bbox = get_bounds_from_las(input_pointcloud) + points_skeleton_clip = clip_points_with_box(points_list_skeleton, bbox) + # Create list with Skeleton's points with Z for step 4 + points_skeleton_with_z_clip = [ + ([geom.x, geom.y, 0]) for geom in points_skeleton_clip.geometry if isinstance(geom, Point) + ] + # Step 4 : 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) + + return result if initial_las_filename: # Lauch croping filtered pointcloud by Mask Hydro with buffer by one tile: @@ -98,13 +114,13 @@ def extract_points_around_skeleton_points_one_tile(filename): points_clip_list = [ extract_points_around_skeleton_points_one_tile(file) for file in os.listdir(input_dir_points) ] - points_clip = np.vstack(points_clip_list) # merge ONE numpy.array - - # Step 4 : Extract a Z elevation value every 2 meters along the hydrographic skeleton - result = filter_las_around_point(points_list_skeleton, points_clip, k) + # Flatten the list of lists into a single list of dictionaries + points_clip = [item for sublist in points_clip_list for item in sublist] - # Convert results to GeoDataFrame - result_gdf = gpd.GeoDataFrame(result) + # 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") From 1abcbaf3051e2e59479d49272dd321e60437182f Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 10 Jul 2024 10:54:00 +0200 Subject: [PATCH 38/64] update --- configs/configs_lidro.yaml | 2 +- example_create_virtual_point_default.sh | 2 +- lidro/create_virtual_point/pointcloud/crop_las.py | 6 ++---- lidro/main_create_virtual_point.py | 2 +- 4 files changed, 5 insertions(+), 7 deletions(-) diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index 13660ac2..24f49e98 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -38,7 +38,7 @@ mask_generation: virtual_point: filter: # Keep ground pointcloud between Hydro Mask and Hydro Mask buffer - keep_classe_ground: "[2:2]" + keep_classes_ground: "[2:2]" vector: # distance in meters between 2 consecutive points distance_meters: 2 diff --git a/example_create_virtual_point_default.sh b/example_create_virtual_point_default.sh index c184bcf7..1ce57db3 100644 --- a/example_create_virtual_point_default.sh +++ b/example_create_virtual_point_default.sh @@ -1,4 +1,4 @@ -# For lauching Create virtual point +# Launch hydro mask merging python -m lidro.main_create_virtual_point \ io.input_dir=./data/ \ io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geosjon \ diff --git a/lidro/create_virtual_point/pointcloud/crop_las.py b/lidro/create_virtual_point/pointcloud/crop_las.py index ba1b017a..e5b0dddc 100644 --- a/lidro/create_virtual_point/pointcloud/crop_las.py +++ b/lidro/create_virtual_point/pointcloud/crop_las.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -""" Crop filtered pointcloud """ +""" Filter pointcloud to keep only the ground class """ from typing import List import numpy as np @@ -8,9 +8,7 @@ def crop_pointcloud_by_points(input_points: str, geom: MultiPolygon, classes: List[int:int]) -> np.array: - """Crop filtered pointcloud : - 1. Filter pointcloud for keeping only classe "GROUND" - 2. Crop filtered pointcloud by "Mask HYDRO + buffer" + """Filter pointcloud to keep only the ground class Args: input_points (str): Path to the input LAS/LAZ file diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 4e2fb48e..4f9abc48 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -62,7 +62,7 @@ def main(config: DictConfig): 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_classe_ground + classes = config.virtual_point.filter.keep_classes_ground k = config.virtual_point.vector.k # Step 1 : Import Mask Hydro, then apply a buffer From ac94d3c4a61d7fc20577ccf61fba1c5eaafd7c57 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 10 Jul 2024 11:02:47 +0200 Subject: [PATCH 39/64] refacto environnement --- environment.yml | 4 +++- requirements.txt | 2 ++ 2 files changed, 5 insertions(+), 1 deletion(-) create mode 100644 requirements.txt diff --git a/environment.yml b/environment.yml index 3d43c141..99769816 100644 --- a/environment.yml +++ b/environment.yml @@ -2,7 +2,7 @@ name: lidro channels: - conda-forge dependencies: - - python==3.11.* + - python==3.10.* - pip - pytest - gdal @@ -24,4 +24,6 @@ dependencies: - flake8 - isort - pre-commit + - pip: + - -r requirements.txt diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 00000000..691ef763 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,2 @@ +ign-pdal-tools==1.5.* + From cf9b8f8f0b27fa6948fa28d24168cd564d892b51 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 10 Jul 2024 11:12:48 +0200 Subject: [PATCH 40/64] refacto environnement with python 3.11.8 --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 99769816..f054933a 100644 --- a/environment.yml +++ b/environment.yml @@ -2,7 +2,7 @@ name: lidro channels: - conda-forge dependencies: - - python==3.10.* + - python==3.11.8 - pip - pytest - gdal From a44a1d7aca1e74b3e49bf1c0375decdfaeabe88c Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 10 Jul 2024 11:21:44 +0200 Subject: [PATCH 41/64] refacto environnement with python 3.11 --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index f054933a..77eb2bde 100644 --- a/environment.yml +++ b/environment.yml @@ -2,7 +2,7 @@ name: lidro channels: - conda-forge dependencies: - - python==3.11.8 + - python==3.11 - pip - pytest - gdal From ea985b0459de5151d04c1c549d4f9f568371b206 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 11 Jul 2024 10:08:23 +0200 Subject: [PATCH 42/64] force geopandas < 1.* --- environment.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index 77eb2bde..21cc2014 100644 --- a/environment.yml +++ b/environment.yml @@ -2,7 +2,7 @@ name: lidro channels: - conda-forge dependencies: - - python==3.11 + - python==3.11.* - pip - pytest - gdal @@ -11,7 +11,7 @@ dependencies: - scipy - geojson - rasterio - - geopandas + - geopandas<=1.* - pyproj - pdal>=2.6 - python-pdal>=3.2.1 From 3fb6f7d09ac1e4f4361accac27cde4f31e53be3a Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 11 Jul 2024 10:10:14 +0200 Subject: [PATCH 43/64] force geopandas==0.* --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 21cc2014..a44b1992 100644 --- a/environment.yml +++ b/environment.yml @@ -11,7 +11,7 @@ dependencies: - scipy - geojson - rasterio - - geopandas<=1.* + - geopandas==0.* - pyproj - pdal>=2.6 - python-pdal>=3.2.1 From 5fcc4435fd8e3cc8539e85408d7637ae5888f6ed Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 11 Jul 2024 10:19:38 +0200 Subject: [PATCH 44/64] rename function crop_las.py, then refacto other related functions --- lidro/create_virtual_point/pointcloud/crop_las.py | 2 +- lidro/main_create_virtual_point.py | 6 ++++-- test/pointcloud/test_crop_las.py | 6 ++++-- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/lidro/create_virtual_point/pointcloud/crop_las.py b/lidro/create_virtual_point/pointcloud/crop_las.py index e5b0dddc..f75a5241 100644 --- a/lidro/create_virtual_point/pointcloud/crop_las.py +++ b/lidro/create_virtual_point/pointcloud/crop_las.py @@ -7,7 +7,7 @@ from shapely.geometry import MultiPolygon -def crop_pointcloud_by_points(input_points: str, geom: MultiPolygon, classes: List[int:int]) -> np.array: +def read_filter_and_crop_pointcloud(input_points: str, geom: MultiPolygon, classes: List[int:int]) -> np.array: """Filter pointcloud to keep only the ground class Args: diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 4f9abc48..91af322e 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -11,7 +11,9 @@ from pyproj import CRS from shapely.geometry import Point -from lidro.create_virtual_point.pointcloud.crop_las import crop_pointcloud_by_points +from lidro.create_virtual_point.pointcloud.crop_las import ( + read_filter_and_crop_pointcloud, +) from lidro.create_virtual_point.pointcloud.utils_pdal import get_bounds_from_las from lidro.create_virtual_point.vectors.clip_points_with_bounding_box import ( clip_points_with_box, @@ -89,7 +91,7 @@ def extract_points_around_skeleton_points_one_tile(filename): input_pointcloud = os.path.join(input_dir_points, filename) # path to the LAS file logging.info(f"\nCroping filtered pointcloud by Mask Hydro with buffer for tile : {tilename}") # Croping filtered pointcloud by Mask Hydro with buffer for tile - points_clip = crop_pointcloud_by_points(input_pointcloud, str(input_mask_hydro_buffer), classes) + points_clip = read_filter_and_crop_pointcloud(input_pointcloud, str(input_mask_hydro_buffer), classes) logging.info(f"\nCropping skeleton points for tile: {tilename}") # Extract bounding box for clipping points by tile bbox = get_bounds_from_las(input_pointcloud) diff --git a/test/pointcloud/test_crop_las.py b/test/pointcloud/test_crop_las.py index b3269f81..35de067c 100644 --- a/test/pointcloud/test_crop_las.py +++ b/test/pointcloud/test_crop_las.py @@ -4,7 +4,9 @@ import numpy as np -from lidro.create_virtual_point.pointcloud.crop_las import crop_pointcloud_by_points +from lidro.create_virtual_point.pointcloud.crop_las import ( + read_filter_and_crop_pointcloud, +) TMP_PATH = Path("./tmp/create_virtual_point/pointcloud/crop_las") LAS_FILE = "./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz" @@ -31,7 +33,7 @@ def test_crop_pointcloud_default(): 830103.2211538461 6290861.586538463, \ 830210.5288461539 6290861.298076924)))" ) - output = crop_pointcloud_by_points(LAS_FILE, geom, classes) + output = read_filter_and_crop_pointcloud(LAS_FILE, geom, classes) assert isinstance(output, np.ndarray) is True assert isinstance((output == 0).all(), bool) is False From 55e6524eae25e5bbaf0bae3822c58363e8118a59 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 11 Jul 2024 10:53:41 +0200 Subject: [PATCH 45/64] delete function utils_pdal.py for using ign-pdal-tools --- .../pointcloud/utils_pdal.py | 35 ------------------- lidro/main_create_virtual_point.py | 4 +-- test/pointcloud/test_utils_pdal.py | 23 ------------ 3 files changed, 2 insertions(+), 60 deletions(-) delete mode 100644 lidro/create_virtual_point/pointcloud/utils_pdal.py delete mode 100644 test/pointcloud/test_utils_pdal.py diff --git a/lidro/create_virtual_point/pointcloud/utils_pdal.py b/lidro/create_virtual_point/pointcloud/utils_pdal.py deleted file mode 100644 index 670ee636..00000000 --- a/lidro/create_virtual_point/pointcloud/utils_pdal.py +++ /dev/null @@ -1,35 +0,0 @@ -# -*- coding: utf-8 -*- -""" Utils PDAL -""" -import pdal - - -def read_las_file(input_las: str): - """Read a las file and put it in an array""" - pipeline = pdal.Pipeline() | pdal.Reader.las(filename=input_las, nosrs=True) - pipeline.execute() - return pipeline.arrays[0] - - -def get_info_from_las(points): - """Get info from a las to put it in an array""" - pipeline = pdal.Filter.stats().pipeline(points) - pipeline.execute() - return pipeline.metadata - - -def get_bounds_from_las(in_las: str): - """Get bounds=([minx,maxx],[miny,maxy]) from las file - - Args: - in_las (str): Path to the input LAS/LAZ file - - Returns: - tuple : Bounds = ([minx,maxx],[miny,maxy]) from las file - """ - metadata = get_info_from_las(read_las_file(in_las)) - xmin = metadata["metadata"]["filters.stats"]["bbox"]["native"]["bbox"]["minx"] - xmax = metadata["metadata"]["filters.stats"]["bbox"]["native"]["bbox"]["maxx"] - ymin = metadata["metadata"]["filters.stats"]["bbox"]["native"]["bbox"]["miny"] - ymax = metadata["metadata"]["filters.stats"]["bbox"]["native"]["bbox"]["maxy"] - return ([xmin, xmax], [ymin, ymax]) diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 91af322e..0660afcd 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -8,13 +8,13 @@ import hydra import pandas as pd from omegaconf import DictConfig +from pdaltools.las_info import las_get_xy_bounds from pyproj import CRS from shapely.geometry import Point from lidro.create_virtual_point.pointcloud.crop_las import ( read_filter_and_crop_pointcloud, ) -from lidro.create_virtual_point.pointcloud.utils_pdal import get_bounds_from_las from lidro.create_virtual_point.vectors.clip_points_with_bounding_box import ( clip_points_with_box, ) @@ -94,7 +94,7 @@ def extract_points_around_skeleton_points_one_tile(filename): points_clip = read_filter_and_crop_pointcloud(input_pointcloud, str(input_mask_hydro_buffer), classes) logging.info(f"\nCropping skeleton points for tile: {tilename}") # Extract bounding box for clipping points by tile - bbox = get_bounds_from_las(input_pointcloud) + bbox = las_get_xy_bounds(input_pointcloud, 0, crs) points_skeleton_clip = clip_points_with_box(points_list_skeleton, bbox) # Create list with Skeleton's points with Z for step 4 points_skeleton_with_z_clip = [ diff --git a/test/pointcloud/test_utils_pdal.py b/test/pointcloud/test_utils_pdal.py deleted file mode 100644 index 7fb8bce8..00000000 --- a/test/pointcloud/test_utils_pdal.py +++ /dev/null @@ -1,23 +0,0 @@ -import numpy as np - -from lidro.create_virtual_point.pointcloud.utils_pdal import ( - get_bounds_from_las, - read_las_file, -) - -LAS_FILE = "./data/pointcloud/LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST.las" - -# Expected : ([minx,maxx],[miny,maxy]) -BoundsExpected_FILE = ([706044.94, 706292.25], [6626159, 6626324]) - - -def test_read_las_file(): - pts = read_las_file(LAS_FILE) - assert isinstance(pts, np.ndarray) # type is array - - -def test_get_bounds_from_las(): - bounds = get_bounds_from_las(LAS_FILE) - - assert isinstance(bounds, tuple) - assert bounds == BoundsExpected_FILE From 2c0e6a29829e7b081d81af510bce219ae3973672 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 11 Jul 2024 11:20:19 +0200 Subject: [PATCH 46/64] rectify function's comments --- lidro/create_virtual_point/stats/calculate_stat.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/lidro/create_virtual_point/stats/calculate_stat.py b/lidro/create_virtual_point/stats/calculate_stat.py index 6207f73e..5e55453c 100644 --- a/lidro/create_virtual_point/stats/calculate_stat.py +++ b/lidro/create_virtual_point/stats/calculate_stat.py @@ -3,9 +3,9 @@ import numpy as np -def calculate_quartile(points: np.array, q: int) -> float: +def calculate_quartile_(points: np.array, q: int) -> float: """ - Calculates the quartile of the value (ex. altitude from z-coordinate) of points. + Calculates the quartile's value of Z coordinates Parameters: - points (numpy.ndarray): An array of 3D points [x, y, z] @@ -13,25 +13,25 @@ def calculate_quartile(points: np.array, q: int) -> float: Values must be between 0 and 100 inclusive. Returns: - - float: The quartile of the value (ex. z-coordinates). + - float: The quartile of Z coordinates """ altitudes = points[:, 2] # Extract the altitude column - n_quartile = round(np.percentile(altitudes, q), 2) + n_quartile = np.round(np.percentile(altitudes, q), 2) return n_quartile def calculate_median(points: np.array) -> float: """ - Calculates the median of the value (ex. altitude from z-coordinate) of points. + Calculates the median's value of Z coordinates Parameters: - points (numpy.ndarray): An array of 3D points [x, y, z] Returns: - - float: The median of the value (ex. z-coordinates). + - float: The median of Z coordinates """ altitudes = points[:, 2] # Extract the altitude column - n_median = round(np.median(altitudes), 2) + n_median = np.round(np.median(altitudes), 2) return n_median From 8b570158172c0c4c56673d1fa5077c4e44294f6a Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 11 Jul 2024 11:53:01 +0200 Subject: [PATCH 47/64] rectify function's comments v2 --- lidro/create_virtual_point/stats/calculate_stat.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lidro/create_virtual_point/stats/calculate_stat.py b/lidro/create_virtual_point/stats/calculate_stat.py index 5e55453c..576ecb6a 100644 --- a/lidro/create_virtual_point/stats/calculate_stat.py +++ b/lidro/create_virtual_point/stats/calculate_stat.py @@ -16,7 +16,7 @@ def calculate_quartile_(points: np.array, q: int) -> float: - float: The quartile of Z coordinates """ altitudes = points[:, 2] # Extract the altitude column - n_quartile = np.round(np.percentile(altitudes, q), 2) + n_quartile = round(np.percentile(altitudes, q), 2) return n_quartile @@ -32,6 +32,6 @@ def calculate_median(points: np.array) -> float: - float: The median of Z coordinates """ altitudes = points[:, 2] # Extract the altitude column - n_median = np.round(np.median(altitudes), 2) + n_median = round(np.median(altitudes), 2) return n_median From fca39eb708c43aa6032687bea144f2cd8a755ae0 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 11 Jul 2024 11:58:02 +0200 Subject: [PATCH 48/64] rectify function's comments v3 --- lidro/create_virtual_point/stats/calculate_stat.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lidro/create_virtual_point/stats/calculate_stat.py b/lidro/create_virtual_point/stats/calculate_stat.py index 576ecb6a..dd2e3dcc 100644 --- a/lidro/create_virtual_point/stats/calculate_stat.py +++ b/lidro/create_virtual_point/stats/calculate_stat.py @@ -3,7 +3,7 @@ import numpy as np -def calculate_quartile_(points: np.array, q: int) -> float: +def calculate_quartile(points: np.array, q: int) -> float: """ Calculates the quartile's value of Z coordinates @@ -16,7 +16,7 @@ def calculate_quartile_(points: np.array, q: int) -> float: - float: The quartile of Z coordinates """ altitudes = points[:, 2] # Extract the altitude column - n_quartile = round(np.percentile(altitudes, q), 2) + n_quartile = np.round(np.percentile(altitudes, q), 2) return n_quartile @@ -32,6 +32,6 @@ def calculate_median(points: np.array) -> float: - float: The median of Z coordinates """ altitudes = points[:, 2] # Extract the altitude column - n_median = round(np.median(altitudes), 2) + n_median = np.round(np.median(altitudes), 2) return n_median From 927152cbaeb85adac8e048e71a218fba7ff86cf2 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 11 Jul 2024 14:01:57 +0200 Subject: [PATCH 49/64] research KNN from 2D points_skeleton NOT 3D --- .../stats/knn_distance.py | 25 +++++++++++-------- lidro/main_create_virtual_point.py | 2 +- test/stats/test_knn_distance.py | 2 +- test/vectors/test_las_around_point.py | 2 +- 4 files changed, 18 insertions(+), 13 deletions(-) diff --git a/lidro/create_virtual_point/stats/knn_distance.py b/lidro/create_virtual_point/stats/knn_distance.py index fe6b3a46..7362bcbe 100644 --- a/lidro/create_virtual_point/stats/knn_distance.py +++ b/lidro/create_virtual_point/stats/knn_distance.py @@ -4,29 +4,34 @@ from sklearn.neighbors import NearestNeighbors -def find_k_nearest_neighbors(point_2d: np.array, points_3d: np.array, k: int) -> np.array: +def find_k_nearest_neighbors(points_skeleton: np.array, points_ground_3d: np.array, k: int) -> np.array: """Finds the K nearest neighbors of a given point from a list of 3D points Args: - point_2d (np.array): An array of 3D point [x, y, 0] - points_3D (np.array): An array of 3D points [x, y, z] - k (int): The number of nearest neighbors to find + points_skeleton (np.array): a single point representing the skeleton Hydro + an array of 2D point [x, y] + points_ground_3d (np.array): ground's points from LIDAR. + an array of 3D points [x, y, z] + k (int): the number of nearest neighbors to find Returns: - numpy.ndarray: An array of the K nearest 3D points [x, y, z] + numpy.ndarray: an array of the K nearest 3D points [x, y, z] """ - # Convert point_2d to nump.arrat - point_2d_array = np.array(point_2d).reshape(1, -1) + # Extract 2D coordinates (x, y) from the 3D points + points_ground_2d = points_ground_3d[:, :2] + + # Convert point_2d to nump.array + points_skeleton_array = np.array(points_skeleton).reshape(1, -1) # Initialize and fit the NearestNeighbors model - nbrs = NearestNeighbors(n_neighbors=k, algorithm="auto").fit(points_3d) + nbrs = NearestNeighbors(n_neighbors=k, algorithm="auto").fit(points_ground_2d) # Find the distances and indices of the nearest neighbors - indices = nbrs.kneighbors(point_2d_array, return_distance=False) + indices = nbrs.kneighbors(points_skeleton_array, return_distance=False) # Retrieve the points corresponding to the indices # Use indices to retrieve the closest 3D points - k_nearest_points = points_3d[indices[0]] + k_nearest_points = points_ground_3d[indices[0]] # = points_3d[indices.flatten()] return k_nearest_points diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 0660afcd..349cab05 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -98,7 +98,7 @@ def extract_points_around_skeleton_points_one_tile(filename): points_skeleton_clip = clip_points_with_box(points_list_skeleton, bbox) # Create list with Skeleton's points with Z for step 4 points_skeleton_with_z_clip = [ - ([geom.x, geom.y, 0]) for geom in points_skeleton_clip.geometry if isinstance(geom, Point) + ([geom.x, geom.y]) for geom in points_skeleton_clip.geometry if isinstance(geom, Point) ] # Step 4 : Extract a Z elevation value along the hydrographic skeleton logging.info(f"\nExtract a Z elevation value along the hydrographic skeleton for tile : {tilename}") diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py index 8f943eb6..32412591 100644 --- a/test/stats/test_knn_distance.py +++ b/test/stats/test_knn_distance.py @@ -12,7 +12,7 @@ def test_find_k_nearest_neighbors_default(): [830867.61, 6290202.62, 2.89], ] ) - point = np.array([[830574.89, 6290648.53, 0]]) + point = np.array([[830574.89, 6290648.53]]) k = 3 result = find_k_nearest_neighbors(point, points_array, k) diff --git a/test/vectors/test_las_around_point.py b/test/vectors/test_las_around_point.py index 083b0ea2..5cdfb756 100644 --- a/test/vectors/test_las_around_point.py +++ b/test/vectors/test_las_around_point.py @@ -22,7 +22,7 @@ def setup_module(module): def test_las_around_point_default(): # Parameters - points_along_skeleton = [[830864.5181373736, 6290217.943739296, 0], [830866.5957826116, 6290208.162525126, 0]] + points_along_skeleton = [[830864.5181373736, 6290217.943739296], [830866.5957826116, 6290208.162525126]] points_clip = np.array( [ From d4cc007d258f97ea950d75c23ea47edbbddb0be1 Mon Sep 17 00:00:00 2001 From: mdupays Date: Thu, 11 Jul 2024 15:26:55 +0200 Subject: [PATCH 50/64] refacto severals scripts for lauching create virtual point: step3 and 4 --- .../vectors/clip_points_with_bounding_box.py | 10 +-- .../vectors/extract_points_around_skeleton.py | 63 +++++++++++++++++++ .../vectors/identify_upstream_downstream.py | 53 ---------------- lidro/main_create_virtual_point.py | 53 ++++------------ .../test_clip_points_with_bounding_box.py | 5 +- .../test_identify_upstream_downstream.py | 43 ------------- 6 files changed, 80 insertions(+), 147 deletions(-) create mode 100644 lidro/create_virtual_point/vectors/extract_points_around_skeleton.py delete mode 100644 lidro/create_virtual_point/vectors/identify_upstream_downstream.py delete mode 100644 test/vectors/test_identify_upstream_downstream.py diff --git a/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py b/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py index 37e4cd10..77293469 100644 --- a/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py +++ b/lidro/create_virtual_point/vectors/clip_points_with_bounding_box.py @@ -1,17 +1,15 @@ # -*- coding: utf-8 -*- """ Clip Skeleton points by tile (pointcloud) """ -from typing import List - import geopandas as gpd -from shapely.geometry import Point, box +from shapely.geometry import box -def clip_points_with_box(points: List, bbox: tuple) -> gpd.GeoDataFrame: +def clip_points_with_box(gdf_points: gpd.GeoDataFrame, bbox: tuple) -> gpd.GeoDataFrame: """Clip skeleton points by tile (bounding box) Args: - points (List): Points every 2 meters (by default) along skeleton hydro + gdf_points (gpd.GeoDataFrame): Points every 2 meters (by default) along skeleton hydro bbox (tuple): bounding box from tile (pointcloud) Returns: @@ -23,8 +21,6 @@ def clip_points_with_box(points: List, bbox: tuple) -> gpd.GeoDataFrame: ymin = bbox[1][0] ymax = bbox[1][1] - # Create a GeoDataFrame from the points - gdf_points = gpd.GeoDataFrame(geometry=[Point(point) for point in points]) # Create a polygon representing the bounding box bounding_box = box(xmin, ymin, xmax, ymax) diff --git a/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py b/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py new file mode 100644 index 00000000..a7a5cbda --- /dev/null +++ b/lidro/create_virtual_point/vectors/extract_points_around_skeleton.py @@ -0,0 +1,63 @@ +# -*- coding: utf-8 -*- +""" Extract points around skeleton by tile +""" +import logging +import os +from typing import List + +import geopandas as gpd +from pdaltools.las_info import las_get_xy_bounds +from shapely.geometry import Point + +from lidro.create_virtual_point.pointcloud.crop_las import ( + read_filter_and_crop_pointcloud, +) +from lidro.create_virtual_point.vectors.clip_points_with_bounding_box import ( + clip_points_with_box, +) +from lidro.create_virtual_point.vectors.las_around_point import filter_las_around_point + + +def extract_points_around_skeleton_points_one_tile( + filename: str, + input_dir: str, + input_mask_hydro_buffer: gpd.GeoDataFrame, + points_skeleton_gdf: gpd.GeoDataFrame, + classes: List[int:int], + k: int, +): + """Severals steps : + 1- Crop filtered pointcloud by Mask Hydro with buffer + 2- Extract a Z elevation value along the hydrographic skeleton + + Args: + filename (str): filename to the LAS file + input_dir (str): input 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 + 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") + tilename = os.path.splitext(filename)[0] # filename to the LAS file + input_pointcloud = os.path.join(input_dir_points, filename) # path to the LAS file + logging.info(f"\nCroping filtered pointcloud by Mask Hydro with buffer for tile : {tilename}") + # Croping filtered pointcloud by Mask Hydro with buffer for tile + points_clip = read_filter_and_crop_pointcloud(input_pointcloud, str(input_mask_hydro_buffer), classes) + logging.info(f"\nCropping skeleton points for tile: {tilename}") + # Extract bounding box for clipping points by tile + bbox = las_get_xy_bounds(input_pointcloud) + points_skeleton_clip = clip_points_with_box(points_skeleton_gdf, bbox) + # Create list with Skeleton's points with Z for step 4 + points_skeleton_with_z_clip = [ + ([geom.x, geom.y]) for geom in points_skeleton_clip.geometry if isinstance(geom, Point) + ] + # 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) + + return result diff --git a/lidro/create_virtual_point/vectors/identify_upstream_downstream.py b/lidro/create_virtual_point/vectors/identify_upstream_downstream.py deleted file mode 100644 index d97d33d4..00000000 --- a/lidro/create_virtual_point/vectors/identify_upstream_downstream.py +++ /dev/null @@ -1,53 +0,0 @@ -# -*- coding: utf-8 -*- -""" Identify upstream and downstream by skeleton's section -""" -from typing import Dict, List, Tuple - - -def calculate_z_q1_averages(points: List) -> Tuple: - """Calculate the average 'Z_q1' values for the first 10 percent and the last 10 percent of a list of points - - Parameters: - points (list): A list of dictionaries where each dictionary contains a 'geometry' key - with a Point object and a 'z_q1' key with a numeric value. - - Returns: - tuple: A pair of values representing the average 'Z_q1' values - for the first 10 percent and the last 10 percent of the list, respectively. - """ - total_points = len(points) - percent_10 = max(1, total_points // 10) # Ensure at least 1 point to avoid division by zero - - # Extract the first 10 percent and the last 10 percent of points - first_10_percent_points = points[:percent_10] - last_10_percent_points = points[-percent_10:] - - # Calculate the average 'Z_q1' values for each group - average_z_q1_first_10_percent = sum(point["z_q1"] for point in first_10_percent_points) / percent_10 - average_z_q1_last_10_percent = sum(point["z_q1"] for point in last_10_percent_points) / percent_10 - - return average_z_q1_first_10_percent, average_z_q1_last_10_percent - - -def identify_upstream_downstream(points: List) -> Dict: - """Determinate point for upstream and downstream by skeleton's section - - Args: - points (List): A list of dictionaries where each dictionary contains a 'geometry' key - with a Point object and a 'z_q1' key with a numeric value. - - Returns: - Dict : A dictionnary representing geometry (Point 2D) and value of "z_q1" - for upstream and downstream by skeleton's section - """ - # Calculate the average 'Z_q1' values for the first 10 percent and the last 10 percent of a list of points - average_z_q1_first_10_percent, average_z_q1_last_10_percent = calculate_z_q1_averages(points) - - if average_z_q1_first_10_percent < average_z_q1_last_10_percent: - upstream_point = points[0] - downstream_point = points[-1] - else: - upstream_point = points[-1] - downstream_point = points[0] - - return upstream_point, downstream_point diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 349cab05..132e8970 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -8,17 +8,11 @@ import hydra import pandas as pd from omegaconf import DictConfig -from pdaltools.las_info import las_get_xy_bounds from pyproj import CRS -from shapely.geometry import Point -from lidro.create_virtual_point.pointcloud.crop_las import ( - read_filter_and_crop_pointcloud, +from lidro.create_virtual_point.vectors.extract_points_around_skeleton import ( + extract_points_around_skeleton_points_one_tile, ) -from lidro.create_virtual_point.vectors.clip_points_with_bounding_box import ( - clip_points_with_box, -) -from lidro.create_virtual_point.vectors.las_around_point import filter_las_around_point from lidro.create_virtual_point.vectors.mask_hydro_with_buffer import ( import_mask_hydro_with_buffer, ) @@ -73,48 +67,23 @@ def main(config: DictConfig): # Step 2 : Create several points every 2 meters (by default) along skeleton Hydro # Return GeoDataframe - points_gdf = generate_points_along_skeleton(input_skeleton, distance_meters, crs) - points_list_skeleton = [([geom.x, geom.y, 0]) for geom in points_gdf.geometry if isinstance(geom, Point)] - - # Step 3 : Crope filtered pointcloud by Mask Hydro with buffer - def extract_points_around_skeleton_points_one_tile(filename): - """Lauch main.py on one tile - - Args: - filename (str): filename to the LAS file - - Returns: - points_clip (np.array) : Numpy array containing point coordinates (X, Y, Z) after filtering and croping - """ - input_dir_points = os.path.join(input_dir, "pointcloud") - tilename = os.path.splitext(filename)[0] # filename to the LAS file - input_pointcloud = os.path.join(input_dir_points, filename) # path to the LAS file - logging.info(f"\nCroping filtered pointcloud by Mask Hydro with buffer for tile : {tilename}") - # Croping filtered pointcloud by Mask Hydro with buffer for tile - points_clip = read_filter_and_crop_pointcloud(input_pointcloud, str(input_mask_hydro_buffer), classes) - logging.info(f"\nCropping skeleton points for tile: {tilename}") - # Extract bounding box for clipping points by tile - bbox = las_get_xy_bounds(input_pointcloud, 0, crs) - points_skeleton_clip = clip_points_with_box(points_list_skeleton, bbox) - # Create list with Skeleton's points with Z for step 4 - points_skeleton_with_z_clip = [ - ([geom.x, geom.y]) for geom in points_skeleton_clip.geometry if isinstance(geom, Point) - ] - # Step 4 : 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) - - return result + 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) + 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) for file in os.listdir(input_dir_points) + extract_points_around_skeleton_points_one_tile( + file, input_dir, input_mask_hydro_buffer, points_skeleton_gdf, classes, k + ) + for file in os.listdir(input_dir_points) ] # Flatten the list of lists into a single list of dictionaries points_clip = [item for sublist in points_clip_list for item in sublist] diff --git a/test/vectors/test_clip_points_with_bounding_box.py b/test/vectors/test_clip_points_with_bounding_box.py index fa916787..d60b42e0 100644 --- a/test/vectors/test_clip_points_with_bounding_box.py +++ b/test/vectors/test_clip_points_with_bounding_box.py @@ -18,6 +18,8 @@ def test_create_hydro_vector_mask_default(): (830375.000, 6291001.000), # Outside (830400.000, 6291000.000), # Outside ] + # Convert points to GeoDataFrame + points_gdf = gpd.GeoDataFrame(geometry=[Point(point) for point in points]) # Define a bounding box that should clip out some points bbox = ((830375, 830390), (6290970, 6291000)) @@ -29,10 +31,9 @@ def test_create_hydro_vector_mask_default(): (830386.045, 6290980.81), (830389.709, 6290971.505), ] - # Convert expected points to GeoDataFrame for comparison expected_gdf = gpd.GeoDataFrame(geometry=[Point(point) for point in expected_points]) - result_gdf = clip_points_with_box(points, bbox) + result_gdf = clip_points_with_box(points_gdf, bbox) assert (result_gdf.equals(expected_gdf)) is True diff --git a/test/vectors/test_identify_upstream_downstream.py b/test/vectors/test_identify_upstream_downstream.py deleted file mode 100644 index c6216efd..00000000 --- a/test/vectors/test_identify_upstream_downstream.py +++ /dev/null @@ -1,43 +0,0 @@ -from typing import Dict - -from shapely.geometry import Point - -from lidro.create_virtual_point.vectors.identify_upstream_downstream import ( - identify_upstream_downstream, -) - - -def test_identify_upstream_downstream_default(): - # Example list of points - points = [ - {"geometry": Point(830211.190425319829956, 6292000.43919328507036), "z_q1": 2.95}, - {"geometry": Point(830211.950744876405224, 6291980.453650656156242), "z_q1": 3.05}, - {"geometry": Point(830212.711064432864077, 6291960.468108027242124), "z_q1": 3.20}, - {"geometry": Point(830213.471383989439346, 6291940.482565398328006), "z_q1": 3.15}, - {"geometry": Point(830214.231703546014614, 6291920.497022769413888), "z_q1": 3.22}, - {"geometry": Point(830214.992023102473468, 6291900.511480140499771), "z_q1": 4.12}, - {"geometry": Point(830218.033301328658126, 6291820.569309624843299), "z_q1": 4.12}, - {"geometry": Point(830218.793620885233395, 6291800.583766995929182), "z_q1": 4.55}, - {"geometry": Point(830219.553940441692248, 6291780.598224367015064), "z_q1": 4.54}, - {"geometry": Point(830224.663624010980129, 6291680.817003472708166), "z_q1": 5.00}, - {"geometry": Point(830228.343047266826034, 6291661.158370648510754), "z_q1": 5.05}, - {"geometry": Point(830232.022470522555523, 6291641.499737825244665), "z_q1": 5.05}, - {"geometry": Point(830246.740163545822725, 6291562.865206529386342), "z_q1": 5.10}, - {"geometry": Point(830250.419586801668629, 6291543.20657370518893), "z_q1": 5.15}, - {"geometry": Point(830254.099010057398118, 6291523.547940881922841), "z_q1": 2.95}, - {"geometry": Point(830257.778433313244022, 6291503.88930805772543), "z_q1": 2.96}, - {"geometry": Point(830261.457856569089927, 6291484.230675233528018), "z_q1": 6.00}, - {"geometry": Point(830265.137279824819416, 6291464.572042410261929), "z_q1": 6.10}, - {"geometry": Point(830268.81670308066532, 6291444.913409586064517), "z_q1": 6.15}, - {"geometry": Point(830272.496126336511225, 6291425.254776761867106), "z_q1": 6.10}, - {"geometry": Point(830287.621702362317592, 6291346.699770421721041), "z_q1": 7.00}, - ] - - upstream, downstream = identify_upstream_downstream(points) - - assert isinstance(upstream, Dict) - assert isinstance(downstream, Dict) - assert upstream["geometry"] == Point(830211.1904253198, 6292000.439193285) - assert upstream["z_q1"] == 2.95 - assert downstream["geometry"] == Point(830287.6217023623, 6291346.699770422) - assert downstream["z_q1"] == 7.0 From b2091aae53e6cb3a5cf409988f71622386980380 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 15 Jul 2024 14:57:55 +0200 Subject: [PATCH 51/64] modify filter.range for script crop_las.py and configs (keep_classes_ground) --- configs/configs_lidro.yaml | 4 ++-- lidro/create_virtual_point/pointcloud/crop_las.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index 24f49e98..a0eb877f 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -37,8 +37,8 @@ mask_generation: virtual_point: filter: - # Keep ground pointcloud between Hydro Mask and Hydro Mask buffer - keep_classes_ground: "[2:2]" + # Keep ground and water pointclouds between Hydro Mask and Hydro Mask buffer + keep_classes_ground: [2, 9] vector: # distance in meters between 2 consecutive points distance_meters: 2 diff --git a/lidro/create_virtual_point/pointcloud/crop_las.py b/lidro/create_virtual_point/pointcloud/crop_las.py index f75a5241..7ee12d97 100644 --- a/lidro/create_virtual_point/pointcloud/crop_las.py +++ b/lidro/create_virtual_point/pointcloud/crop_las.py @@ -7,13 +7,13 @@ from shapely.geometry import MultiPolygon -def read_filter_and_crop_pointcloud(input_points: str, geom: MultiPolygon, classes: List[int:int]) -> np.array: +def read_filter_and_crop_pointcloud(input_points: str, geom: MultiPolygon, classes: List[int]) -> np.array: """Filter pointcloud to keep only the ground class Args: input_points (str): Path to the input LAS/LAZ file geom (MultiPolygon): An array of WKT or GeoJSON 2D MultiPolygon (Mask Hydro with buffer) - classes (list): List of classes to use for the filtering + classes (List[int]): Values to keep for input points along filter_dimension Returns: np.array : Numpy array containing point coordinates (X, Y, Z) after filtering and croping @@ -22,7 +22,7 @@ def read_filter_and_crop_pointcloud(input_points: str, geom: MultiPolygon, class pipeline = ( pdal.Reader.las(filename=input_points, nosrs=True) | pdal.Filter.range( - limits=f"Classification{classes}", + limits=",".join(f"Classification[{v}:{v}]" for v in classes), ) | pdal.Filter.crop( polygon=geom, From 2b3f6ef3041b0c9dc1b2fc2ee89252ed9665a232 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 15 Jul 2024 15:00:59 +0200 Subject: [PATCH 52/64] create folder (scripts) for putting severals examples scripts --- .../example_create_mask_by_tile.sh | 0 .../example_create_mask_default.sh | 0 .../example_create_virtual_point_by_tile.sh | 0 .../example_create_virtual_point_default.sh | 0 .../example_merge_mask_default.sh | 0 5 files changed, 0 insertions(+), 0 deletions(-) rename example_create_mask_by_tile.sh => scripts/example_create_mask_by_tile.sh (100%) rename example_create_mask_default.sh => scripts/example_create_mask_default.sh (100%) rename example_create_virtual_point_by_tile.sh => scripts/example_create_virtual_point_by_tile.sh (100%) rename example_create_virtual_point_default.sh => scripts/example_create_virtual_point_default.sh (100%) rename example_merge_mask_default.sh => scripts/example_merge_mask_default.sh (100%) diff --git a/example_create_mask_by_tile.sh b/scripts/example_create_mask_by_tile.sh similarity index 100% rename from example_create_mask_by_tile.sh rename to scripts/example_create_mask_by_tile.sh diff --git a/example_create_mask_default.sh b/scripts/example_create_mask_default.sh similarity index 100% rename from example_create_mask_default.sh rename to scripts/example_create_mask_default.sh diff --git a/example_create_virtual_point_by_tile.sh b/scripts/example_create_virtual_point_by_tile.sh similarity index 100% rename from example_create_virtual_point_by_tile.sh rename to scripts/example_create_virtual_point_by_tile.sh diff --git a/example_create_virtual_point_default.sh b/scripts/example_create_virtual_point_default.sh similarity index 100% rename from example_create_virtual_point_default.sh rename to scripts/example_create_virtual_point_default.sh diff --git a/example_merge_mask_default.sh b/scripts/example_merge_mask_default.sh similarity index 100% rename from example_merge_mask_default.sh rename to scripts/example_merge_mask_default.sh From 6b6dec0257e64dd62ff3e10203107187421d0509 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 15 Jul 2024 15:04:59 +0200 Subject: [PATCH 53/64] modify test (crop_las.py) to List[int] --- test/pointcloud/test_crop_las.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/pointcloud/test_crop_las.py b/test/pointcloud/test_crop_las.py index 35de067c..6ec8680f 100644 --- a/test/pointcloud/test_crop_las.py +++ b/test/pointcloud/test_crop_las.py @@ -19,7 +19,7 @@ def setup_module(module): def test_crop_pointcloud_default(): - classes = "[1:2]" + classes = [2] geom = str( "MULTIPOLYGON (((830873.1249999998 6290475.625000002, \ From ac845339afbcb2042760ccc1d310937132115d4f Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 15 Jul 2024 15:08:39 +0200 Subject: [PATCH 54/64] add test (check dimension) from test_knn_distance.py --- test/stats/test_knn_distance.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py index 32412591..3a174776 100644 --- a/test/stats/test_knn_distance.py +++ b/test/stats/test_knn_distance.py @@ -18,3 +18,4 @@ def test_find_k_nearest_neighbors_default(): result = find_k_nearest_neighbors(point, points_array, k) assert isinstance(result, np.ndarray) is True + assert result.ndim == 3 # check dimension numpy.array -> 3D From 1ccd79dc03dc99ad097aea7d8f3e4d4130f6e8cd Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 15 Jul 2024 15:16:32 +0200 Subject: [PATCH 55/64] rename function calculate_quartile to calculate_percentile --- .../stats/calculate_stat.py | 4 +- .../vectors/las_around_point.py | 4 +- test/stats/test_calculate_stats.py | 43 ++++++++----------- 3 files changed, 22 insertions(+), 29 deletions(-) diff --git a/lidro/create_virtual_point/stats/calculate_stat.py b/lidro/create_virtual_point/stats/calculate_stat.py index dd2e3dcc..8bfd0a1c 100644 --- a/lidro/create_virtual_point/stats/calculate_stat.py +++ b/lidro/create_virtual_point/stats/calculate_stat.py @@ -3,13 +3,13 @@ import numpy as np -def calculate_quartile(points: np.array, q: int) -> float: +def calculate_percentile(points: np.array, q: int) -> float: """ Calculates the quartile's value of Z coordinates Parameters: - points (numpy.ndarray): An array of 3D points [x, y, z] - - q (int): Percentage or sequence of percentages for the percentiles to compute. + - q (int): Percentage for the percentiles to compute. Values must be between 0 and 100 inclusive. Returns: diff --git a/lidro/create_virtual_point/vectors/las_around_point.py b/lidro/create_virtual_point/vectors/las_around_point.py index a0be698a..ece0dec5 100644 --- a/lidro/create_virtual_point/vectors/las_around_point.py +++ b/lidro/create_virtual_point/vectors/las_around_point.py @@ -6,7 +6,7 @@ import numpy as np from shapely.geometry import Point -from lidro.create_virtual_point.stats.calculate_stat import calculate_quartile +from lidro.create_virtual_point.stats.calculate_stat import calculate_percentile from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors @@ -29,7 +29,7 @@ def filter_las_around_point(points_skeleton: List, points_clip: np.array, k: int # Calcule Z "Q1" for each points every N meters along skeleton hydro results = [ - ({"geometry": p["geometry"], "z_q1": calculate_quartile(p["points_knn"], 10)}) + ({"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"])) ] diff --git a/test/stats/test_calculate_stats.py b/test/stats/test_calculate_stats.py index 0e0c4c08..3d43af73 100644 --- a/test/stats/test_calculate_stats.py +++ b/test/stats/test_calculate_stats.py @@ -1,49 +1,42 @@ import numpy as np +import pytest from lidro.create_virtual_point.stats.calculate_stat import ( calculate_median, - calculate_quartile, + calculate_percentile, ) -def test_calculate_quartile_25(): +def test_calculate_percentile_25(): points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) - assert calculate_quartile(points, 25) == 6.0 + assert calculate_percentile(points, 25) == 6.0 -def test_calculate_quartile_50(): +def test_calculate_percentile_50(): points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) - assert calculate_quartile(points, 50) == 9.0 + assert calculate_percentile(points, 50) == 9.0 -def test_calculate_quartile_75(): +def test_calculate_percentile_75(): points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) - assert calculate_quartile(points, 75) == 12.0 + assert calculate_percentile(points, 75) == 12.0 -def test_calculate_quartile_100(): +def test_calculate_percentile_100(): points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) - assert calculate_quartile(points, 100) == 15.0 + assert calculate_percentile(points, 100) == 15.0 -def test_calculate_quartile_invalid_q(): - points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) - try: - calculate_quartile(points, 110) - except ValueError: - pass - else: - assert False, "ValueError non levée pour q=110" +def test_calculate_percentile_invalid_q(): + with pytest.raises(ValueError): + points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) + calculate_percentile(points, 110) -def test_calculate_quartile_negative_q(): - points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) - try: - calculate_quartile(points, -10) - except ValueError: - pass - else: - assert False, "ValueError non levée pour q=-10" +def test_calculate_percentile_negative_q(): + with pytest.raises(ValueError): + points = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]]) + calculate_percentile(points, -10) def test_calculate_median(): From d34e15e0ab3f855458159610707f2a506913a32d Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 15 Jul 2024 15:23:06 +0200 Subject: [PATCH 56/64] add test (check dimension) from test_knn_distance.py : V2 --- test/stats/test_knn_distance.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py index 3a174776..65d0c6d9 100644 --- a/test/stats/test_knn_distance.py +++ b/test/stats/test_knn_distance.py @@ -18,4 +18,8 @@ def test_find_k_nearest_neighbors_default(): result = find_k_nearest_neighbors(point, points_array, k) assert isinstance(result, np.ndarray) is True - assert result.ndim == 3 # check dimension numpy.array -> 3D + assert result.ndim == [ + [830438.91, 6290854.32, 2.56], + [830721.84, 6290447.79, 2.23], + [830861.04, 6290242.06, 2.78], + ] From fa4befa5cec4229e0dbec4f889a5de87433aaad4 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 15 Jul 2024 15:52:31 +0200 Subject: [PATCH 57/64] add test (check dimension) from test_knn_distance.py : V2bis --- test/stats/test_knn_distance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py index 65d0c6d9..851c701f 100644 --- a/test/stats/test_knn_distance.py +++ b/test/stats/test_knn_distance.py @@ -18,7 +18,7 @@ def test_find_k_nearest_neighbors_default(): result = find_k_nearest_neighbors(point, points_array, k) assert isinstance(result, np.ndarray) is True - assert result.ndim == [ + assert result == [ [830438.91, 6290854.32, 2.56], [830721.84, 6290447.79, 2.23], [830861.04, 6290242.06, 2.78], From acf7ae752984ccf26176da7cb5d17f9f19792303 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 15 Jul 2024 16:26:01 +0200 Subject: [PATCH 58/64] add test (check dimension) from test_knn_distance.py : V3 --- test/stats/test_knn_distance.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py index 851c701f..66258c69 100644 --- a/test/stats/test_knn_distance.py +++ b/test/stats/test_knn_distance.py @@ -18,8 +18,4 @@ def test_find_k_nearest_neighbors_default(): result = find_k_nearest_neighbors(point, points_array, k) assert isinstance(result, np.ndarray) is True - assert result == [ - [830438.91, 6290854.32, 2.56], - [830721.84, 6290447.79, 2.23], - [830861.04, 6290242.06, 2.78], - ] + assert all(len(sublist) == 3 for sublist in result) # check 3d dimension From 54c35f4f0cbb8f4d55a2422e7a7a3e6055925931 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 17 Jul 2024 10:46:52 +0200 Subject: [PATCH 59/64] refacto create mask buffer hydro (limit bank + n meters) --- lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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 531baa8c..c8ab157b 100644 --- a/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py +++ b/lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py @@ -22,4 +22,7 @@ def import_mask_hydro_with_buffer(file_mask: str, buffer: float, crs: str | int) # Apply buffer (2 meters by default) from Mask Hydro gdf_buffer = gdf.buffer(buffer, cap_style=CAP_STYLE.square) - return gdf_buffer + # Return a polygon representing the limit of the bank with a buffer of N meters + limit_bank = gdf_buffer.difference(gdf) + + return limit_bank From d1b74a047443a2c1065cca13a3eec110dce45912 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 17 Jul 2024 10:49:51 +0200 Subject: [PATCH 60/64] modify names keep_classes_ground to keep_neighbors_classes --- configs/configs_lidro.yaml | 2 +- lidro/main_create_virtual_point.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index a0eb877f..2ad66f58 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -38,7 +38,7 @@ mask_generation: virtual_point: filter: # Keep ground and water pointclouds between Hydro Mask and Hydro Mask buffer - keep_classes_ground: [2, 9] + keep_neighbors_classes: [2, 9] vector: # distance in meters between 2 consecutive points distance_meters: 2 diff --git a/lidro/main_create_virtual_point.py b/lidro/main_create_virtual_point.py index 132e8970..fc610166 100644 --- a/lidro/main_create_virtual_point.py +++ b/lidro/main_create_virtual_point.py @@ -58,7 +58,7 @@ def main(config: DictConfig): 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_classes_ground + classes = config.virtual_point.filter.keep_neighbors_classes k = config.virtual_point.vector.k # Step 1 : Import Mask Hydro, then apply a buffer From 788b6cb577b6fe8267abb7fbd9c064eab04e86a6 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 17 Jul 2024 10:55:27 +0200 Subject: [PATCH 61/64] modify definition : classes --- lidro/create_virtual_point/pointcloud/crop_las.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lidro/create_virtual_point/pointcloud/crop_las.py b/lidro/create_virtual_point/pointcloud/crop_las.py index 7ee12d97..424e0081 100644 --- a/lidro/create_virtual_point/pointcloud/crop_las.py +++ b/lidro/create_virtual_point/pointcloud/crop_las.py @@ -13,7 +13,7 @@ def read_filter_and_crop_pointcloud(input_points: str, geom: MultiPolygon, class Args: input_points (str): Path to the input LAS/LAZ file geom (MultiPolygon): An array of WKT or GeoJSON 2D MultiPolygon (Mask Hydro with buffer) - classes (List[int]): Values to keep for input points along filter_dimension + classes (List[int]): Classification values to keep for input points Returns: np.array : Numpy array containing point coordinates (X, Y, Z) after filtering and croping From fe97c9afef5945990c7220cf4760e8716a18f0a7 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 17 Jul 2024 11:06:19 +0200 Subject: [PATCH 62/64] modify definition q --- lidro/create_virtual_point/stats/calculate_stat.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lidro/create_virtual_point/stats/calculate_stat.py b/lidro/create_virtual_point/stats/calculate_stat.py index 8bfd0a1c..d781f2ef 100644 --- a/lidro/create_virtual_point/stats/calculate_stat.py +++ b/lidro/create_virtual_point/stats/calculate_stat.py @@ -9,7 +9,7 @@ def calculate_percentile(points: np.array, q: int) -> float: Parameters: - points (numpy.ndarray): An array of 3D points [x, y, z] - - q (int): Percentage for the percentiles to compute. + - q (int): Percentage for the percentile to compute. Values must be between 0 and 100 inclusive. Returns: From c920aa5f401428493f906eef060bc03d03e91555 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 17 Jul 2024 11:15:37 +0200 Subject: [PATCH 63/64] refacto test with classes=[1, 2] --- test/pointcloud/test_crop_las.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/pointcloud/test_crop_las.py b/test/pointcloud/test_crop_las.py index 6ec8680f..9f1fac5b 100644 --- a/test/pointcloud/test_crop_las.py +++ b/test/pointcloud/test_crop_las.py @@ -19,7 +19,7 @@ def setup_module(module): def test_crop_pointcloud_default(): - classes = [2] + classes = [1, 2] geom = str( "MULTIPOLYGON (((830873.1249999998 6290475.625000002, \ From a0192f6d08984e1eda15e7f332126854dd320b05 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 17 Jul 2024 11:27:52 +0200 Subject: [PATCH 64/64] add test (check dimension) from test_knn_distance.py : term --- test/stats/test_knn_distance.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/stats/test_knn_distance.py b/test/stats/test_knn_distance.py index 66258c69..122368c8 100644 --- a/test/stats/test_knn_distance.py +++ b/test/stats/test_knn_distance.py @@ -13,9 +13,9 @@ def test_find_k_nearest_neighbors_default(): ] ) point = np.array([[830574.89, 6290648.53]]) - k = 3 + k = 1 result = find_k_nearest_neighbors(point, points_array, k) assert isinstance(result, np.ndarray) is True - assert all(len(sublist) == 3 for sublist in result) # check 3d dimension + assert np.array_equal(result, np.array([[830438.91, 6290854.32, 2.56]]))