upgrade laspy to 2.2.0 with laszip support, allows to directly read and write laszip file

opendm/dem/commands.py

@@ -50,36 +50,13 @@ def rectify(lasFile, debug=False, reclassify_threshold=5, min_area=750, min_poin
     start = datetime.now()
 
     try:
-        # Currently, no Python 2 lib that supports reading and writing LAZ, so we will do it manually until ODM is migrated to Python 3
-        # When migration is done, we can move to pylas and avoid using PDAL for conversion
-        tempLasFile = os.path.join(os.path.dirname(lasFile), 'tmp.las')
-
-        # Convert LAZ to LAS
-        cmd = [
-            'pdal',
-            'translate',
-            '-i %s' % lasFile,
-            '-o %s' % tempLasFile
-        ]
-        system.run(' '.join(cmd))
-
         log.ODM_INFO("Rectifying {} using with [reclassify threshold: {}, min area: {}, min points: {}]".format(lasFile, reclassify_threshold, min_area, min_points))
         run_rectification(
-            input=tempLasFile, output=tempLasFile, debug=debug, \
+            input=lasFile, output=lasFile, debug=debug, \
             reclassify_plan='median', reclassify_threshold=reclassify_threshold, \
             extend_plan='surrounding', extend_grid_distance=5, \
             min_area=min_area, min_points=min_points)
 
-        # Convert LAS to LAZ
-        cmd = [
-            'pdal',
-            'translate',
-            '-i %s' % tempLasFile,
-            '-o %s' % lasFile
-        ]
-        system.run(' '.join(cmd))
-        os.remove(tempLasFile)
-
     except Exception as e:
         raise Exception("Error rectifying ground in file %s: %s" % (lasFile, str(e)))
 

opendm/dem/ground_rectification/extra_dimensions/distance_dimension.py

@@ -35,7 +35,7 @@ class DistanceDimension(Dimension):
         return 'distance_to_ground'
 
     def get_las_type(self):
-        return 10
+        return 'float64'
 
     def __calculate_angle(self, model):
         "Calculate the angle between the estimated plane and the XY plane"

opendm/dem/ground_rectification/extra_dimensions/extended_dimension.py

@@ -20,4 +20,4 @@ class ExtendedDimension(Dimension):
         return 'extended'
 
     def get_las_type(self):
-        return 3
+        return 'uint16'

opendm/dem/ground_rectification/extra_dimensions/partition_dimension.py

@@ -22,4 +22,4 @@ class PartitionDimension(Dimension):
         return self.name
 
     def get_las_type(self):
-        return 5
+        return 'uint32'

opendm/dem/ground_rectification/io/las_io.py

@@ -1,52 +1,39 @@
 # TODO: Move to pylas when project migrates to python3
 
-from laspy.file import File
-from laspy.header import Header
+import laspy
 import numpy as np
 from ..point_cloud import PointCloud
 
 def read_cloud(point_cloud_path):
     # Open point cloud and read its properties
-    las_file = File(point_cloud_path, mode='r')
-    header = (las_file.header.copy(), las_file.header.scale, las_file.header.offset,las_file.header.evlrs, las_file.header.vlrs)
-    [x_scale, y_scale, z_scale] = las_file.header.scale
-    [x_offset, y_offset, z_offset] = las_file.header.offset
+    las_file = laspy.read(point_cloud_path)
+    header = las_file.header
 
-    # Calculate the real coordinates
-    x = las_file.X * x_scale + x_offset
-    y = las_file.Y * y_scale + y_offset
-    z = las_file.Z * z_scale + z_offset
+    x = las_file.x.scaled_array()
+    y = las_file.y.scaled_array()
+    z = las_file.z.scaled_array()
 
-    cloud = PointCloud.with_dimensions(x, y, z, las_file.Classification, las_file.red, las_file.green, las_file.blue)
-
-    # Close the file
-    las_file.close()
+    cloud = PointCloud.with_dimensions(x, y, z, las_file.classification.array, las_file.red, las_file.green, las_file.blue)
 
     # Return the result
     return header, cloud
 
 def write_cloud(header, point_cloud, output_point_cloud_path, write_extra_dimensions=False):
-    (h, scale, offset, evlrs, vlrs) = header
-
     # Open output file
-    output_las_file = File(output_point_cloud_path, mode='w', header=h, evlrs=evlrs, vlrs=vlrs)
+    output_las_file = laspy.LasData(header)
 
     if write_extra_dimensions:
-        # Create new dimensions
-        for name, dimension in point_cloud.extra_dimensions_metadata.items():
-            output_las_file.define_new_dimension(name=name, data_type=dimension.get_las_type(), description="Dimension added by Ground Extend")
-
+        extra_dims = [laspy.ExtraBytesParams(name=name, type=dimension.get_las_type(), description="Dimension added by Ground Extend") for name, dimension in point_cloud.extra_dimensions_metadata.items()]
+        output_las_file.add_extra_dims(extra_dims)
         # Assign dimension values
         for dimension_name, values in point_cloud.extra_dimensions.items():
             setattr(output_las_file, dimension_name, values)
 
     # Adapt points to scale and offset
-    [x_scale, y_scale, z_scale] = scale
-    [x_offset, y_offset, z_offset] = offset
     [x, y] = np.hsplit(point_cloud.xy, 2)
-    output_las_file.X = (x.ravel() - x_offset) / x_scale
-    output_las_file.Y = (y.ravel() - y_offset) / y_scale
-    output_las_file.Z = (point_cloud.z - z_offset) / z_scale
+    output_las_file.x = x.ravel()
+    output_las_file.y = y.ravel()
+    output_las_file.z = point_cloud.z
 
     # Set color
     [red, green, blue] = np.hsplit(point_cloud.rgb, 3)
@@ -55,11 +42,6 @@ def write_cloud(header, point_cloud, output_point_cloud_path, write_extra_dimens
     output_las_file.blue = blue.ravel()
 
     # Set classification
-    output_las_file.Classification = point_cloud.classification.astype(np.uint8)
+    output_las_file.classification = point_cloud.classification.astype(np.uint8)
 
-    # Set header
-    output_las_file.header.scale = scale
-    output_las_file.header.offset = offset
-
-    # Close files
-    output_las_file.close()
+    output_las_file.write(output_point_cloud_path)

requirements.txt

@@ -7,7 +7,7 @@ ODMExifRead==3.0.4
 Fiona==1.8.17 ; sys_platform == 'linux' or sys_platform == 'darwin'
 https://github.com/OpenDroneMap/windows-deps/raw/main/Fiona-1.8.19-cp38-cp38-win_amd64.whl ; sys_platform == 'win32'
 joblib==1.1.0
-laspy==1.7.0
+laspy[laszip]==2.3.0
 lxml==4.6.1
 matplotlib==3.3.3
 networkx==2.5