refacto severals scripts for lauching create virtual point: step3 and 4

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)
 

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

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]

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