Add virtual points to tiles

configs/configs_lidro.yaml

@@ -11,6 +11,7 @@ io:
   input_mask_hydro: null
   input_skeleton: null
   input_dir: null
+  input_dir_point_virtual: null
   output_dir: null
   dir_points_skeleton: null # folder to contains files (.GeoJSON) by LIDAR tiles : neighboring points for each point of the skeleton
   srid: 2154

lidro/create_virtual_point/pointcloud/add_virtual_points_to_pointcloud.py

@@ -0,0 +1,34 @@
+# -*- coding: utf-8 -*-
+""" Add virtual points by tiles"""
+import json
+
+import pdal
+from shapely.geometry import Polygon
+
+
+def add_virtual_points_by_tiles(input_file: str, input_las: str, output_laz_file: str, geom: Polygon):
+    """Add virtual points (3D point grid in LAZ format) by LIDAR tiles (tiling file)
+
+    Args:
+        input_file (str): Path to the input virtual points file (LAZ)
+        input_las (str): Path to the LIDAR tiles (LAZ)
+        output_dir (str): Path to the output LIDAR tiles (LAZ) with virtual points
+        geom (Polygon):  Shapely Polygon (LIDAR tiles)
+    """
+    # Add virtual points by LIDAR tiles
+    pipeline = [
+        {"type": "readers.las", "filename": input_file, "nosrs": True},
+        {"type": "filters.crop", "polygon": geom.wkt},  # Convert the Shapely Polygon to WKT
+        {"type": "readers.las", "filename": input_las, "nosrs": True},
+        {"type": "filters.merge"},
+        {
+            "type": "writers.las",
+            "filename": output_laz_file,
+            "compression": True,
+            "minor_version": 4,
+            "dataformat_id": 6,
+        },
+    ]
+    # Create and execute the pipeline
+    p = pdal.Pipeline(json.dumps(pipeline))
+    p.execute()

lidro/create_virtual_point/pointcloud/auto_tiling_from_las.py

@@ -2,6 +2,7 @@
 """ This function automatically creates the tiling of lidar tiles (.GeoJSON)
 """
 import logging
+import math
 import os
 from typing import List
 
@@ -78,7 +79,8 @@ def create_geojson_from_laz_files(laz_files: List, output_geojson_path: str, crs
     # Create a list of dictionaries representing the tiles
     tiles = [
         {
-            "tile_id": os.path.basename(laz_file),
+            "tile_id": str(int(min_x)) + "_" + str(math.ceil(max_y / 1000) * 1000),
+            "tilename_las": os.path.basename(laz_file),
             "geometry": Polygon([(min_x, min_y), (max_x, min_y), (max_x, max_y), (min_x, max_y)]),
         }
         for laz_file in laz_files

lidro/create_virtual_point/vectors/run_add_virtual_points_by_tile.py

@@ -0,0 +1,37 @@
+# -*- coding: utf-8 -*-
+""" Run function "create virtual points by tiles"
+"""
+import json
+import logging
+import os
+
+from shapely.geometry import shape
+
+from lidro.create_virtual_point.pointcloud.add_virtual_points_to_pointcloud import (
+    add_virtual_points_by_tiles,
+)
+
+
+def compute_virtual_points_by_tiles(input_virtual_points: str, tile_geojson: str, input_dir: str, output_dir: str):
+    """
+    Loops through the tiles in the GeoJSON file and clips a LAZ file into multiple files based on polygons.
+
+    Args:
+        input_virtual_points (str): Path to the input virtual points file
+        tile_geojson (str): Path to the GeoJSON file containing tile polygons
+        input_dir (str): folder to input LIDAR tiles
+        output_dir (str): folder to output virtual points by LIDAR tiles
+    """
+    # Load the GeoJSON tile geometries
+    with open(tile_geojson, "r") as f:
+        geojson_data = json.load(f)
+
+    # Clip the virtual points by each tile
+    for feature in geojson_data["features"]:
+        if feature.get("geometry", {}).get("type") == "Polygon":
+            name = feature["properties"]["tilename_las"]
+            input_las = os.path.join(input_dir, name)
+            output_file = os.path.join(output_dir, name)
+            shapely_polygon = shape(feature["geometry"])
+            add_virtual_points_by_tiles(input_virtual_points, input_las, output_file, shapely_polygon)
+            logging.info(f"Clip virtual points by tiles : {input_las}")

lidro/create_virtual_point/vectors/run_create_virtual_points.py

@@ -21,7 +21,7 @@
 )
 
 
-def launch_virtual_points_by_section(
+def compute_virtual_points_by_section(
     points: gpd.GeoDataFrame,
     line: gpd.GeoDataFrame,
     mask_hydro: gpd.GeoDataFrame,

lidro/main_clip_virtual_point_by_tile.py

@@ -0,0 +1,59 @@
+""" Main script for create virtuals points
+"""
+
+import logging
+import os
+
+import hydra
+from omegaconf import DictConfig
+from pyproj import CRS
+
+from lidro.create_virtual_point.pointcloud.auto_tiling_from_las import (
+    create_geojson_from_laz_files,
+)
+from lidro.create_virtual_point.vectors.run_add_virtual_points_by_tile import (
+    compute_virtual_points_by_tiles,
+)
+
+
+@hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2")
+def main(config: DictConfig):
+    """Create a virtual point inside hydro surfaces (3D grid) from the points classification of
+    the input LAS/LAZ file and the Hyro Skeleton (GeoJSON) and save it as LAS file.
+
+    It can run either on a single file, or on each file of a folder
+
+    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)
+
+    # Parameters for clip virtual point by tiles
+    crs = CRS.from_user_input(config.io.srid)
+    input_dir_point_virtual = config.io.input_dir_point_virtual
+
+    # Clip virtual points file by LIDAR tiles
+    # Create the tiling of lidar tiles
+    json_tiles = os.path.join(output_dir, "tiles_from_las.GeoJSON")
+    create_geojson_from_laz_files([os.path.join(input_dir, file) for file in os.listdir(input_dir)], json_tiles, crs)
+    # Clip virtual points (3D point grid in LAZ format) by LIDAR tiles (tiling file)
+    virtul_points_file = os.path.join(input_dir_point_virtual, "virtual_points.laz")
+    compute_virtual_points_by_tiles(virtul_points_file, json_tiles, input_dir, output_dir)
+
+
+if __name__ == "__main__":
+    main()

lidro/main_create_virtual_point.py

@@ -4,23 +4,21 @@
 import ast
 import logging
 import os
-import sys
+
 import geopandas as gpd
 import hydra
 import pandas as pd
 from omegaconf import DictConfig
 from pyproj import CRS
 
-sys.path.append('../lidro')
-
 from lidro.create_virtual_point.pointcloud.convert_list_points_to_las import (
     list_points_to_las,
 )
 from lidro.create_virtual_point.vectors.merge_skeleton_by_mask import (
     merge_skeleton_by_mask,
 )
 from lidro.create_virtual_point.vectors.run_create_virtual_points import (
-    launch_virtual_points_by_section,
+    compute_virtual_points_by_section,
 )
 
 
@@ -84,7 +82,7 @@ def process_points_knn(points_knn):
 
         # Step 3 : Generate a regular grid of 3D points spaced every N meters inside each hydro entity
         list_virtual_points = [
-            launch_virtual_points_by_section(
+            compute_virtual_points_by_section(
                 points_gdf,
                 gpd.GeoDataFrame([{"geometry": row["geometry_skeleton"]}], crs=crs),
                 gpd.GeoDataFrame([{"geometry": row["geometry_mask"]}], crs=crs),
@@ -96,6 +94,7 @@ def process_points_knn(points_knn):
             for idx, row in gdf_merged.iterrows()
         ]
         logging.info("Calculate virtuals points by mask hydro and skeleton")
+
         # Step 4 : Save the virtual points in a file (.LAZ)
         list_points_to_las(list_virtual_points, output_dir, crs, classes)
     else:

scripts/example_clip_virtual_point_default .sh

@@ -0,0 +1,8 @@
+# Launch hydro mask merging
+python -m lidro.main_clip_virtual_point_by_tile \
+io.input_dir=./data/pointcloud/ \
+io.input_dir_point_virtual=./data/point_virtual/ \
+io.output_dir=./tmp/ \
+
+
+

scripts/example_create_mask_default.sh

@@ -2,7 +2,7 @@
 python -m lidro.main_create_mask \
 io.input_dir=./data/pointcloud/ \
 io.output_dir=./tmp/ \
-io.mask_generation.raster.dilatation_size=3 \
-io.mask_generation.filter.keep_classes="[0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67]" \ 
+mask_generation.raster.dilatation_size=3 \
+mask_generation.filter.keep_classes="[0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67]" \ 
 
 

scripts/example_create_virtual_point_by_tile.sh

@@ -5,6 +5,6 @@ io.input_filename=Semis_2021_0830_6291_LA93_IGN69.laz \
 io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geojson \
 io.input_skeleton=./data/skeleton_hydro/Skeleton_Hydro.geojson \
 io.output_dir=./tmp/ \
-io.virtual_point.pointcloud.points_grid_spacing=0.5 \
+virtual_point.pointcloud.points_grid_spacing=0.5 \
 
 

scripts/example_create_virtual_point_default.sh

@@ -5,7 +5,7 @@ io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geojson \
 io.input_skeleton=./data/skeleton_hydro/Skeleton_Hydro.geojson \
 io.dir_points_skeleton=./tmp/point_skeleton/ \
 io.output_dir=./tmp/ \
-io.virtual_point.pointcloud.points_grid_spacing=0.5 \
+virtual_point.pointcloud.points_grid_spacing=0.5 \
 
 
 

scripts/example_extract_points_around_skeleton_default.sh

@@ -4,10 +4,10 @@ io.input_dir=./data/ \
 io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geosjon \
 io.input_skeleton=./data/skeleton_hydro/Skeleton_Hydro.geojson \
 io.output_dir=./tmp/points_skeleton/ \
-io.virtual_point.vector.distance_meters=5 \
-io.virtual_point.vector.buffer=3 \
-io.virtual_point.vector.k=3 \
-io.virtual_point.filter.keep_neighbors_classes="[2, 9]" \ 
+virtual_point.vector.distance_meters=5 \
+virtual_point.vector.buffer=3 \
+virtual_point.vector.k=3 \
+virtual_point.filter.keep_neighbors_classes="[2, 9]" \ 
 
 
 

scripts/example_merge_mask_default.sh

@@ -2,10 +2,10 @@
 python -m lidro.main_merge_mask \
 io.input_dir=./data/mask_hydro/ \
 io.output_dir=./tmp/merge_mask_hydro/ \
-io.mask_generation.vector.min_water_area=150 \
-io.mask_generation.vector.buffer_positive=1 \
-io.mask_generation.vector.buffer_negative=-1.5 \
-io.mask_generation.vector.tolerance=1 \
+mask_generation.vector.min_water_area=150 \
+mask_generation.vector.buffer_positive=1 \
+mask_generation.vector.buffer_negative=-1.5 \
+mask_generation.vector.tolerance=1 \
 
 
 

test/pointcloud/test_add_virtual_points_to_pointcloud.py

@@ -0,0 +1,59 @@
+import os
+import shutil
+from pathlib import Path
+
+import pdal
+from shapely.geometry import Polygon
+
+from lidro.create_virtual_point.pointcloud.add_virtual_points_to_pointcloud import (
+    add_virtual_points_by_tiles,
+)
+
+# Path to temporary directory for output files
+TMP_PATH = Path("./tmp/virtual_points/pointcloud/add_virtual_points_by_tiles")
+input_file = str("./data/point_virtual/virtual_points.laz")
+input_las = str("./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz")
+output_laz_file = str(TMP_PATH / "output_0830_6291_with_virtual_points.laz")
+geom = Polygon(
+    [
+        (830000.0, 6290000.0),
+        (831000.0, 6290000.0),
+        (831000.0, 6291000.0),
+        (830000.0, 6291000.0),
+        (830000.0, 6290000.0),
+    ]
+)
+
+
+def setup_module():
+    """Setup function to ensure the temporary directory is clean before tests."""
+    if TMP_PATH.is_dir():
+        shutil.rmtree(TMP_PATH)
+    os.makedirs(TMP_PATH)
+
+
+def test_add_virtual_points_by_tiles():
+    """Test that the output file is in the correct LAZ format."""
+    add_virtual_points_by_tiles(input_file, input_las, output_laz_file, geom)
+
+    # Check that the output file exists
+    assert os.path.exists(output_laz_file), "Output LAZ file should be created."
+
+    # Verify the format using PDAL
+    pipeline_input = pdal.Pipeline() | pdal.Reader.las(input_file, nosrs=True)
+    pipeline_input.execute()
+
+    pipeline_virtual_point = pdal.Pipeline() | pdal.Reader.las(input_las, nosrs=True)
+    pipeline_virtual_point.execute()
+
+    pipeline_output = pdal.Pipeline() | pdal.Reader.las(output_laz_file, nosrs=True)
+    pipeline_output.execute()
+
+    # Check if the point count for input, virtual points and output
+    point_count_input = pipeline_input.metadata["metadata"]["readers.las"]["count"]
+    point_count_virtual_point = pipeline_virtual_point.metadata["metadata"]["readers.las"]["count"]
+    point_count_output = pipeline_output.metadata["metadata"]["readers.las"]["count"]
+
+    assert point_count_output == int(
+        point_count_input + point_count_virtual_point
+    ), "The output LAZ file should contain points."

test/pointcloud/test_auto_tiling_from_las.py

@@ -6,7 +6,6 @@
 import pytest
 from shapely.geometry import Polygon
 
-# Import the function from your module
 from lidro.create_virtual_point.pointcloud.auto_tiling_from_las import (
     create_geojson_from_laz_files,
     extract_bounds_from_laz,

test/test_main_clip_virtual_point_by_tile.py

@@ -0,0 +1,49 @@
+import os
+import subprocess as sp
+from pathlib import Path
+
+from hydra import compose, initialize
+
+from lidro.main_clip_virtual_point_by_tile import main
+
+INPUT_DIR = Path("data/pointcloud")
+INPUT_DIR_POINT_VIRTUAL = Path("data/point_virtual")
+OUTPUT_DIR = Path("tmp") / "clip_virtual_point_by_tile/main"
+
+
+def setup_module(module):
+    os.makedirs("tmp/clip_virtual_point_by_tile/main", exist_ok=True)
+
+
+def test_main_run_okay():
+    repo_dir = Path.cwd().parent
+    cmd = f"""python -m lidro.main_clip_virtual_point_by_tile \
+        io.input_dir="{repo_dir}/lidro/data/pointcloud/"\
+        io.input_dir_point_virtual="{repo_dir}/lidro/data/point_virtual/"\
+        io.output_dir="{repo_dir}/lidro/tmp/clip_virtual_point_by_tile/main/"
+        """
+    sp.run(cmd, shell=True, check=True)
+
+
+def test_main_clip_virtual_point_by_tile():
+    input_dir = INPUT_DIR
+    input_dir_point_virtual = INPUT_DIR_POINT_VIRTUAL
+    output_dir = OUTPUT_DIR / "main_clip_virtual_point_by_tile"
+    input_filename = "Semis_2021_0830_6291_LA93_IGN69.laz"
+    srid = 2154
+
+    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.input_dir_point_virtual={input_dir_point_virtual}",
+                f"io.output_dir={output_dir}",
+                f"io.srid={srid}",
+            ],
+        )
+    main(cfg)
+    assert (Path(output_dir) / "tiles_from_las.GeoJSON").is_file()
+    assert (Path(output_dir) / input_filename).is_file()