I’m working on a number of large .las files (>300M points) from LiDAR scans where I have to perform some calculations on a subset of the points in those files. Reading the files all at once is problematic due to the large memory use when reading all the data into memory, making the processing extremely slow. I’m not looking for a solution that writes the files to disk (e.g. chunk writing), but rather something that returns a LasData object with the same dimensions/point format as the original .las file, but with a subset of the points.

Often I only need a smaller amount of the total points, and the selection can be based on filtering based on dimension values (such as a range of intensity), clipping to a polygon area or simply thinning the point cloud. The size of the point cloud after filtering is not known beforehand, so I cannot preallocate an array of a known, final, size.

I’ve come up with two solutions that appear to produce the correct result (see check_expected_points and the expected results from read_las_baseline, set the variable TEST_LASDATA = True to run test), and are both more memory efficient and faster compared to reading all at once as in read_las_baseline. It’s exemplified with decimating the point cloud, but the goal is to use this is cases with additional filtering steps and without being able to know the final point count. The function read_las_chunks_filter_preallocate performs best, but being new to .las data and a python novice in general I wonder if there are still better/faster ways to do this and if this is a proper way to handle .las data.

The example can be ran using simple.las from the laspy repository to check for correctness, but for performance testing one needs larger files. Run generate_las to create a larger file on disk.

import timeit
import gc

import laspy # version 2.5.4
import numpy as np # version 2.0.2

LAS_FILE = r"large_dummy_file.las" # generated using generate_las_100M(), or use "simple.las" from https://github.com/laspy/laspy/tree/master/tests/data
DECIMATE_FACTOR=5
TEST_LASDATA = False # whether to run the tests comparing the data with check_expected_points(), requires extra runs and might lead to less accurate timing results

def generate_las(output_path='large_dummy_file.las', n_points=100_000_000):
    """ Creates a test .las file with 100M points (about 3.5GB disk space required) """
    # Taken from https://laspy.readthedocs.io/en/latest/examples.html#creating-a-new-lasdata
    SHAPE = int(n_points**0.5)
    # 0. Creating some dummy data
    my_data_xx, my_data_yy = np.meshgrid(np.linspace(-20, 20, SHAPE), np.linspace(-20, 20, SHAPE))
    my_data_zz = my_data_xx ** 2 + 0.25 * my_data_yy ** 2
    my_data = np.hstack((my_data_xx.reshape((-1, 1)), my_data_yy.reshape((-1, 1)), my_data_zz.reshape((-1, 1))))
    
    # 1. Create a new header
    header = laspy.LasHeader(point_format=3, version="1.2")
    header.add_extra_dim(laspy.ExtraBytesParams(name="random", type=np.int32))
    header.offsets = np.min(my_data, axis=0)
    header.scales = np.array([0.1, 0.1, 0.1])
    
    # 2. Create a Las
    las = laspy.LasData(header)
    
    las.x = my_data[:, 0]
    las.y = my_data[:, 1]
    las.z = my_data[:, 2]
    las.random = np.random.randint(-1503, 6546, len(las.points), np.int32)
    
    las.write(output_path)

def check_expected_points(true_las: laspy.LasData, las_to_test: laspy.LasData):
    """ Compares two laspy.LasData objects. Based on https://github.com/laspy/laspy/blob/master/tests/test_chunk_read_write.py """
    assert true_las.header.point_count == las_to_test.header.point_count
    assert true_las.header.point_format == las_to_test.header.point_format
    
    np.testing.assert_array_equal(true_las.header.offsets, las_to_test.header.offsets)
    np.testing.assert_array_equal(true_las.header.scales, las_to_test.header.scales)
    
    expected_points = true_las.points
    to_test_points = las_to_test.points
    for dim_name in to_test_points.array.dtype.names:
        assert np.allclose(
            expected_points[dim_name], to_test_points[dim_name]
        ), f"{dim_name} not equal"
            
        
def read_las_baseline(las_file, decimate):
    """ Read and decimate without reading in chunks"""
    las = laspy.read(las_file)
    las.points = las.points[::decimate]

    return las

def read_las_chunks_filter_preallocate(las_file, decimate=None):
    """ This function uses pre-allocated PackedPointRecord that of the full size, then slices is to the reduced
    size afterwards """
    
    CHUNK_SIZE = 1_000_000
    
    with laspy.open(las_file) as f:
        point_record = laspy.PackedPointRecord.zeros(f.header.point_count, f.header.point_format)
        new_header = f.header

        current_insert_index = 0
        for points in f.chunk_iterator(CHUNK_SIZE):
            # can manipulate points here..
            # e.g. filter on angle, intensity, in polygon etc
            # the size of points after filtering is not known beforehand in final application
            if decimate:
                points = points[::decimate]
            
            chunk_arr_len = points.array.shape[0]
            point_record.array[current_insert_index:current_insert_index+chunk_arr_len] = points.array
            current_insert_index += chunk_arr_len
         
        # slice to the actual size of inserted data, and update the header
        point_record = point_record[:current_insert_index]
        new_header.point_count=len(point_record)
 
        output_las = laspy.LasData(header=new_header, points=point_record)
    
    return output_las


def read_las_chunks_filter_list_concat(las_file, decimate=None):
    """ This function stores the filtered points.array in a list, then in the end concatenates the points
    and uses these to create a new LasData object.
    """

    CHUNK_SIZE = 1_000_000
    with laspy.open(las_file) as f:
        filtered_points = []
        final_point_record = laspy.PackedPointRecord.empty(f.header.point_format)
        
        for points in f.chunk_iterator(CHUNK_SIZE):
            # can manipulate points here..
            # e.g. filter on angle, intensity, in polygon etc
            # the size of points after filtering is not known beforehand in final application
            
            if decimate:
                points = points[::decimate]
            
            filtered_points.append(points.array)
            
        concatenated_points = np.concatenate(filtered_points)
        final_point_record.array = concatenated_points
        output_las = laspy.LasData(header=f.header)
        output_las.points = final_point_record # setting points here instead of LasData call will set correct point_count
    
    return output_las

def main():
    
    methods = [
        ('read_las_baseline', read_las_baseline),
        ('read_las_chunks_filter_preallocate', read_las_chunks_filter_preallocate),
        ('read_las_chunks_filter_list_concat', read_las_chunks_filter_list_concat),
    ]
    
    if TEST_LASDATA:
        expected = read_las_baseline(LAS_FILE, DECIMATE_FACTOR)
   
    for name, method in methods:
        
        if TEST_LASDATA:
            result = method(LAS_FILE, decimate=DECIMATE_FACTOR)
            check_expected_points(expected, result)
            del result
            gc.collect()
        # timing a single run
        t = timeit.Timer(lambda: method(LAS_FILE, DECIMATE_FACTOR))
        print(f"{name}: {t.timeit(number=1):.1f} seconds")

if __name__ == '__main__':
    main()

The timing difference is large when hitting memory limits, if there is enough memory there are no large differences in timing. However, the memory use is quite different (not shown here).

Case when maxing out available RAM, large difference in processing time:

read_las_baseline: 38.3 seconds
read_las_chunks_filter_preallocate: 4.1 seconds
read_las_chunks_filter_list_concat: 5.9 seconds

Sufficient RAM available, less of a difference:

read_las_baseline: 6.0 seconds
read_las_chunks_filter_preallocate: 2.9 seconds
read_las_chunks_filter_list_concat: 3.4 seconds

I’ve not formally benchmarked the memory use, but simply observed it as the code runs. The read_las_chunks_filter_preallocate function is far superior when it comes to memory use.

read_las_chunks_filter_preallocate: The function which preallocates an array of the input size, and then slices to the size of the filtered points. This is the fastest method, probably mainly due to the very low memory use which seems to be proportional to the amount of points that are kept. I wonder though, if this is a proper way to do this with laspy, and if there is a more efficient way to preallocate the data container? It feels a bit hacky, and I wonder if there are pitfalls I’m not seeing when dealing with .las data.

read_las_chunks_filter_list_concat: Stores the filtered points.array in a list when doing chunk reading, and in the end of concatenates the arrays in the list. This uses much more memory compared to read_las_chunks_filter_preallocate, and is a little slower for large point clouds. Is there a way to make this more memory efficient, or a better way to concatenate the arrays?