Merge pull request #103 from IGNF/add_points_with_pdal and version 1.9.1

Add points with pdal

Author: gliegard

CHANGELOG.md

@@ -1,3 +1,7 @@
+
+# 1.9.1
+- las_add_points_to_pointcloud: Fix add points to LAS (use PDAL instead of Laspy)
+
 # 1.9.0
 - custom PDAL: in the docker image, compile custom PDAL (waiting for PDAL 2.9)
 

pdaltools/_version.py

@@ -1,4 +1,4 @@
-__version__ = "1.9.0"
+__version__ = "1.9.1"
 
 
 if __name__ == "__main__":

pdaltools/add_points_in_pointcloud.py

@@ -1,5 +1,6 @@
 import argparse
 from shutil import copy2
+import tempfile
 
 import geopandas as gpd
 import laspy
@@ -10,6 +11,8 @@
 
 from pdaltools.las_info import get_epsg_from_las, get_tile_bbox
 
+import pdal
+
 
 def parse_args(argv=None):
     parser = argparse.ArgumentParser("Add points from GeoJSON in LIDAR tile")
@@ -127,13 +130,12 @@ def add_points_to_las(
         crs (str): CRS of the data.
         virtual_points_classes (int): The classification value to assign to those virtual points (default: 66).
     """
-    # Copy data pointcloud
-    copy2(input_las, output_las)
 
     if input_points_with_z.empty:
         print(
             "No points to add. All points of the geojson file are outside the tile. Copying the input file to output"
         )
+        copy2(input_las, output_las)
         return
 
     # Extract XYZ coordinates and additional attribute (classification)
@@ -148,24 +150,30 @@ def add_points_to_las(
         header = las.header
         if not header:
             header = laspy.LasHeader(point_format=8, version="1.4")
+
+    new_points = laspy.ScaleAwarePointRecord.zeros(nb_points, header=header)  # use header for input_las
+    # then fill in the gaps (X, Y, Z an classification)
+    new_points.x = x_coords.astype(new_points.x.dtype)
+    new_points.y = y_coords.astype(new_points.y.dtype)
+    new_points.z = z_coords.astype(new_points.z.dtype)
+    new_points.classification = classes.astype(new_points.classification.dtype)
+
+    with tempfile.NamedTemporaryFile(suffix="_new_points.las") as tmp:
+        with laspy.open(tmp.name, mode="w", header=header) as las_file:
+            las_file.write_points(new_points)
+
         if crs:
-            try:
-                crs_obj = CRS.from_user_input(crs)  # Convert to a pyproj.CRS object
-            except CRSError:
-                raise ValueError(f"Invalid CRS: {crs}")
-            header.add_crs(crs_obj)
-
-    # Add the new points with 3D points
-    with laspy.open(output_las, mode="a", header=header) as output_las:  # mode `a` for adding points
-        # create nb_points points with "0" everywhere
-        new_points = laspy.ScaleAwarePointRecord.zeros(nb_points, header=header)  # use header for input_las
-        # then fill in the gaps (X, Y, Z an classification)
-        new_points.x = x_coords.astype(new_points.x.dtype)
-        new_points.y = y_coords.astype(new_points.y.dtype)
-        new_points.z = z_coords.astype(new_points.z.dtype)
-        new_points.classification = classes.astype(new_points.classification.dtype)
-
-        output_las.append_points(new_points)
+            a_srs = crs
+        else:
+            a_srs = get_epsg_from_las(input_las)
+
+        # Use pdal to merge the new points with the existing points
+        pipeline = pdal.Pipeline()
+        pipeline |= pdal.Reader.las(filename=input_las)
+        pipeline |= pdal.Reader.las(filename=tmp.name)
+        pipeline |= pdal.Filter.merge()
+        pipeline |= pdal.Writer.las(filename=output_las, forward="all", a_srs=a_srs)
+        pipeline.execute()
 
 
 def line_to_multipoint(line, spacing: float, z_value: float = None):

test/test_add_points_in_pointcloud.py

@@ -272,6 +272,22 @@ def test_add_points_from_geometry_to_las(input_file, input_points, epsg, expecte
         input_points, input_file, OUTPUT_FILE, 68, epsg, 1000, spacing, altitude_column
     )
     assert Path(OUTPUT_FILE).exists()  # check output exists
+    
+    # Read input and output files to compare headers
+    input_las = laspy.read(input_file)
+    output_las = laspy.read(OUTPUT_FILE)
+    
+    # Compare headers
+    assert input_las.header.version == output_las.header.version
+    assert input_las.header.system_identifier == output_las.header.system_identifier
+    assert input_las.header.extra_header_bytes == output_las.header.extra_header_bytes
+    assert input_las.header.extra_vlr_bytes == output_las.header.extra_vlr_bytes
+    assert input_las.header.number_of_evlrs == output_las.header.number_of_evlrs
+    assert input_las.header.point_format == output_las.header.point_format
+    assert np.array_equal(input_las.header.scales, output_las.header.scales)
+    assert np.array_equal(input_las.header.offsets, output_las.header.offsets)
+    assert input_las.header.vlrs[0].string == output_las.header.vlrs[0].string
+    
     point_count = compute_count_one_file(OUTPUT_FILE)["68"]
     assert point_count == expected_nb_points  # Add all points from geojson