使用Python读取las点云，写入las点云，无损坐标精度

1 为什么要写这个博文

搜索使用python读写las点云数据，可以找到很多结果。但是！ 有些只是简单的demo，且没有发现/说明可能遇到的问题；有些晦涩难懂，不够实用主义；有些须付费观看，没有开源精神。

本人在使用laspy v2.3读写点云文件时，着实被坐标精度问题难到了，顺便就仔细学习了下las格式到底是什么来头。

您猜怎么着，如果能打开这个网址，都在这里面详细说明了：laspy document

我还是种个树吧，直接拔了就能用。文末贴源码。

2 提出一些关键问题

注意：
* 实测版本 laspy v2.3，2.x版本应该都可用，但不适合1.x版本。
* las 存储数据时，需要设置 scales 和 offsets，否则会出现精度问题。
* las 存储颜色时，数值类型为 16 位无符号整型。rgb = (normal_rgb * 65535).astype(np.uint16)
las 格式原生支持的属性字段，与las版本密切相关。官方说明：https://laspy.readthedocs.io/en/latest/intro.html#point-records
对 scales 和 offsets 的理解：
比例scales 表明数据的准确性。 0.001 是毫米精度。这意味着如果您的坐标是 0.123456，它将被限制为 0.123。
偏移offset 的目的是避免整数溢出。假设您要存储 123456789.123。在 LAS 文件中，您实际上将存储一个整数：123456789123，该整数将在读取时使用比例因子转换为 123456789.123。但 123456789123 比 32 位整数所能存储的要大得多。因此存储时，将该值偏移 123450000，实际存的是6789123。
(6789123 * 0.001 + 123450000 = 123456789.123)
段落出处

3 给出全部代码

安装依赖

使用 pip （# 选择一个安装就好）：

# Install _without_ LAZ support
pip install laspy# Install with LAZ support via lazrs
pip install laspy[lazrs]# Install with LAZ support via laszip
pip install laspy[laszip]# Install with LAZ support via both lazrs & laszip
pip install laspy[lazrs,laszip]

使用 conda （# 选择一个安装就好）：

conda install -c conda-forge laspy
conda install -c conda-forge lazrs-python

源码（laspy v2.x）

三个工具函数：
read_las_fit() : 读取 las 文件
write_las_fit(): 保存 las 文件
get_las_header_attrs() : 获取不同 las 版本支持的固有属性
备注：_fit 的意思是，可以支持各种属性信息的读取和写入

import laspy
import numpy as npdef read_las_fit(filename, attrs=None):"""读取 las 文件，获取三维坐标 xyz, 颜色 rgb, 属性 attr_dict。当文件没有 RGB 信息时，返回全0的 RGB 信息Args:filename: <str> las 文件路径attrs: <list> 需要额外获取的属性信息 如 ['label']Returns:xyz, rgb, attr_dict"""if attrs is None:attrs = []# 默认返回 scales, offsets ，合并 ["scales", "offsets"]attrs = list(set(attrs + ["scales", "offsets"]))# 读取点云inFile = laspy.read(filename)# inFile.point_format.dimensions可以获取所有的维度信息N_points = len(inFile)x = np.reshape(inFile.x, (N_points, 1))y = np.reshape(inFile.y, (N_points, 1))z = np.reshape(inFile.z, (N_points, 1))xyz = np.hstack((x, y, z))# TODO 注意。如果是大写的 X Y Z，需要转换后才是真实坐标: real_x = scale[0] * inFile.X + offset[0]# 初始化 rgb 全是 0rgb = np.zeros((N_points, 3), dtype=np.uint16)if hasattr(inFile, "red") and hasattr(inFile, "green") and hasattr(inFile, "blue"):r = np.reshape(inFile.red, (N_points, 1))g = np.reshape(inFile.green, (N_points, 1))b = np.reshape(inFile.blue, (N_points, 1))# i = np.reshape(inFile.Reflectance, (N_points, 1))rgb = np.hstack((r, g, b))else:print(f"注意：{filename.split('/')[-1]} 没有RGB信息，返回全0的RGB信息！")# 组织其他属性信息attr_dict = {}for attr in attrs:# 先判断 header 中是否有该属性if hasattr(inFile.header, attr):value = getattr(inFile.header, attr)if hasattr(value, "array"):attr_dict[attr] = np.array(value)else:attr_dict[attr] = value# 再判断 是否为额外属性elif hasattr(inFile, attr):value = getattr(inFile, attr)if hasattr(value, "array"):attr_dict[attr] = np.array(value)else:attr_dict[attr] = valueelse:attr_dict[attr] = Noneprint(f"注意：{filename.split('/')[-1]} 没有属性 {attr} 信息！")return xyz, rgb, attr_dictdef write_las_fit(out_file, xyz, rgb=None, attrs=None):"""将点云数据写入 las 文件，支持写入 坐标xyz, 颜色rgb, 属性attrsArgs:out_file: 输出文件路径xyz: 点云坐标 ndarray (N, 3)rgb: 点云颜色 ndarray (N, 3)attrs:固有属性：file_source_id, gps_time, Intensity, Number of Returns, ....额外属性：label, pred, ...注意：如果不传入 scales 和 offsets，则会自动计算Returns:"""if attrs is None:attrs = {}# 1. 创建 las 文件头。point_format和version决定了las支持哪些固有属性# 详情见 https://pylas.readthedocs.io/en/latest/intro.html?highlight=red#point-recordsheader = laspy.LasHeader(point_format=7, version="1.4")  # 7 支持rgb# 自动计算 scales 和 offsets，确保坐标精度无损# https://stackoverflow.com/questions/77308057/conversion-accuracy-issues-of-e57-to-las-in-python-using-pye57-and-laspyif "offset" not in attrs:min_offset = np.floor(np.min(xyz, axis=0))attrs["offset"] = min_offsetif "scales" not in attrs:attrs["scales"] = [0.001, 0.001, 0.001]  # 0.001 是毫米精度# 初始化一些需要保存的属性值。如果是固有属性，直接赋值; 如果是额外属性，添加到 header 中, 后续赋值extra_attr = []for attr, value in attrs.items():if hasattr(header, attr):  # 设置固有的属性的值, 如 scales, offsetsheader.__setattr__(attr, value)else:  # 添加额外属性，在 las 初始化后赋值header.add_extra_dim(laspy.ExtraBytesParams(name=attr, type=np.float32))extra_attr.append(attr)# 2. 创建 las 文件las = laspy.LasData(header)# 添加xyz坐标las.x = xyz[:, 0]las.y = xyz[:, 1]las.z = xyz[:, 2]# 添加RGB颜色，如果是归一化的颜色，则需要乘以 65535，转为 uint16if rgb is not None:if np.max(rgb) <= 1:rgb = (rgb * 65535).astype(np.uint16)  # 65535 = 2^16 - 1, las存储颜色是16位无符号整型las.red = rgb[:, 0]las.green = rgb[:, 1]las.blue = rgb[:, 2]# 添加额外属性for attr in extra_attr:# 当 value 是 n * 1 的 ndarray 时，转换为 1 维数组value = attrs[attr]if value.ndim == 2 and value.shape[1] == 1:value = value.flatten()las[attr] = value# 保存LAS文件las.write(out_file)def get_las_header_attrs(point_format=7, version="1.4"):"""根据 point_format 和 version 获取 las 文件的 header 属性说明文档：https://laspy.readthedocs.io/en/latest/intro.html#point-recordsArgs:point_format: 点格式version: 版本Returns:"""dimensions = []header = laspy.LasHeader(point_format=point_format, version=version)  # 7 支持rgbfor dim in header.point_format.dimensions:dimensions.append(dim.name)return dimensionsif __name__ == '__main__':# 测试1： 获取 las 文件头属性fields = get_las_header_attrs(7, "1.4")print(f"point_format=7, version=1.4, 头文件包含字段： {fields}")# 测试2： 读取LAS文件read_las_path = "/path_2_data/one_point_cloud.las"xyz_, rgb_, attrs_ = read_las_fit(read_las_path, ["scales", "offsets"])print(attrs_)# 测试3： 写入LAS文件save_las_path = "/path_2_data/one_point_cloud_fit.las"write_las_fit(save_las_path, xyz_, rgb_, {# "scales": attrs_["scales"],# "offsets": attrs_["offsets"]})