使用laspy高效分块读取并过滤大型LAS文件的优化方案咨询及实现合理性验证

我目前在处理超大型LiDAR扫描生成的LAS文件（点数超过3亿），需要对文件中的部分点进行计算。一次性读取整个文件会占用大量内存，导致处理速度极慢。我不需要将处理后的文件写入磁盘的方案（比如分块写入），而是希望得到一个和原LAS文件维度/点格式一致，但仅包含子集点的LasData对象。

通常我只需要总点数中的一小部分，筛选条件可能基于维度值范围（比如强度区间）、多边形区域裁剪，或者直接对点云进行抽稀。过滤后的点云大小无法预先得知，所以我不能预分配已知最终大小的数组。

我已经想出了两种似乎能得到正确结果的方案（可通过设置TEST_LASDATA = True运行测试，对比check_expected_points和read_las_baseline的预期结果），这两种方案相比read_las_baseline中的一次性读取方式，内存效率更高、速度更快。示例中用点云抽稀来演示，但目标是将其用于包含更多过滤步骤、且无法提前知道最终点数的场景。read_las_chunks_filter_preallocate的性能表现最佳，但由于我是LAS数据新手，同时也是Python初学者，我想知道是否还有更好/更快的实现方式，以及这种处理LAS数据的方式是否规范。

可以使用laspy仓库中的simple.las来验证结果的正确性，但性能测试需要更大的文件。可以运行generate_las函数在磁盘上生成大文件。

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 laspy repository 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 laspy's official examples
    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 laspy's 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()

性能测试结果 #

当内存被占满时，处理时间差异非常大：

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

当内存充足时，差异相对较小：

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

我没有做正式的内存使用基准测试，只是在代码运行时观察内存占用情况。read_las_chunks_filter_preallocate函数的内存占用优势非常明显。

核心疑问 #

1. 关于read_las_chunks_filter_preallocate：
   这个函数预分配一个和输入文件大小相同的数组，然后切片到过滤后的实际点数。这是速度最快的方法，可能主要是因为内存占用极低，且仅与保留的点数成正比。但我想知道，这是否是使用laspy处理LAS数据的规范方式？有没有更高效的预分配数据容器的方法？这种方式感觉有点“取巧”，我担心处理LAS数据时会遇到没注意到的陷阱。

2. 关于read_las_chunks_filter_list_concat：
   这个函数把过滤后的points.array存储在列表中，最后将列表中的数组拼接起来，再创建新的LasData对象。相比前一种方法，它的内存占用高很多，处理大型点云时速度也稍慢。有没有办法让这种方式的内存效率更高？或者有没有更好的数组拼接方式？

备注：内容来源于stack exchange，提问作者rhkarls