使用Python读取las点云，写入las点云，无损坐标精度
"""
update: 2025/3/7
已经发布到 **** ：/qq_27816785/article/details/138634496
使用 laspy v2.3
1. 注意：
* 实测版本 laspy v2.3，版本应该都可用，但不适合版本。
* las 存储数据时，需要设置 scales 和 offsets，否则会出现精度问题。
* las 存储颜色时，数值类型为 16 位无符号整型。rgb = (normal_rgb * 65535).astype(np.uint16)
2. las 格式原生支持的属性字段，与las版本密切相关。官方说明：/en/latest/#point-records
3. 对 scales 和 offsets 的理解：
/questions/346679/if-converting-xyz-to-las-using-e-g-laspy-then-how-to-specify-offset-and-scal
比例scales 表明数据的准确性。 0.001 是毫米精度。这意味着如果您的坐标是 0.123456，它将被限制为 0.123。
偏移offset 的目的是避免整数溢出。假设您要存储 123456789.123。在 LAS 文件中，您实际上将存储一个整数：123456789123，
该整数将在读取时使用比例因子转换为 123456789.123。但 123456789123 比 32 位整数所能存储的要大得多。
因此存储时，将该值偏移 123450000，实际存的是6789123。(6789123 * 0.001 + 123450000 = 123456789.123)
所以 xmin、ymax、zmin 是有效的偏移值选择，但我们通常更喜欢将这些值四舍五入到最接近的整数米/公里
4. 请注意，读取las后，再重新保存las时，最好使用原文件的 scales 和 offsets，以保持精度。因为这2个参数都会影响新处理的点是否与原始点一致。
scales 保证了在 0.001 毫米级别下的精度一致。
但 offsets 可能会导致毫米级别以下的精度损失，因此造成小于毫米级的精度损失。
"""
import warnings
import laspy
import numpy as np
def read_las_fit(filename, user_attrs=None, ignore_default=True):
"""
读取 las 文件，获取三维坐标 xyz, 颜色 rgb, 属性 attr_dict。
当文件没有 RGB 信息时，返回全0的 RGB 信息
当 user_attrs 缺省时，返回所有的属性（忽略无效值）
Args:
filename: <str> las 文件路径
user_attrs: <list> 需要额外获取的属性信息 如 ['label']
Returns:
xyz, rgb, attr_dict
"""
# 读取点云数据
inFile = laspy.read(filename)
# 读取header信息
version, point_format = get_las_info(filename)
# 获取支持的属性字段
all_fields = get_las_header_attrs(point_format, version)
print(f" --读取 {filename.split('/')[-1]} 文件，版本={version}, 点格式={point_format}")
print(f" --支持的属性字段：{all_fields}")
if user_attrs is None:
# 获取所有可用的维度信息，但不包括 X Y Z red green blue
not_need = ["X", "Y", "Z", "red", "green", "blue"]
user_attrs = [dim.name for dim in inFile.point_format.dimensions if dim.name not in not_need]
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))
''' 注意。如果是大写的 X Y Z，需要转换后才是真实坐标: real_x = scale[0] * + offset[0] '''
# 初始化 rgb 全是 0
rgb = 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 = (, (N_points, 1))
rgb = np.hstack((r, g, b))
else:
print(f"注意：{filename.split('/')[-1]} 没有RGB信息，返回全0的RGB信息！")
# 无论 user_attrs 是否需要，都会返回 point_format、version、scales、offsets
attr_dict = {
"point_format": point_format,
"version": version,
"scales": inFile.header.scale,
"offsets": inFile.header.offset
}
for attr in user_attrs:
value = None
# 先判断 是否为额外属性
if hasattr(inFile, attr):
value = getattr(inFile, attr)
print(f" -获取有效属性：{attr}") # 包括固有属性，如 gps_time, intensity, 以及额外属性，如 label, pred
# 再判断 header 中是否有该属性
elif hasattr(inFile.header, attr):
value = getattr(inFile.header, attr)
print(f" --获取header属性：{attr}")
else:
warnings.warn(f"{filename.split('/')[-1]} 没有属性= {attr} 的信息！", category=Warning) # 使用 warnning 警告
# 获取到值后，如果是 n * 1 的 ndarray 时，转换为 1 维数组
if hasattr(value, "array"):
value = np.array(value)
# 判断值是否有效，如果是无效值，则不添加到 attr_dict 中
if ignore_default:
if value is not None and not np.all(value == 0) and user_attrs is not None:
attr_dict[attr] = value
else:
attr_dict[attr] = value
print()
return xyz, rgb, attr_dict
def write_las_fit(out_file, xyz, rgb=None, attrs=None):
"""
将点云数据写入 las 文件，支持写入 坐标xyz, 颜色rgb, 属性attrs
Args:
out_file: 输出文件路径
xyz: 点云坐标 ndarray (N, 3)
rgb: 点云颜色 ndarray (N, 3)
attrs: 字典{key:value, ...}
固有属性：file_source_id, gps_time, Intensity, Number of Returns, ....
额外属性：label, pred, ...
注意：如果不传入 scales 和 offsets，则会自动计算
Returns:
"""
if attrs is None:
attrs = {}
''' 创建 las 文件头。point_format和version决定了las支持哪些固有属性 '''
# 详情见 /en/latest/?highlight=red#point-records
point_format = attrs.pop("point_format") if "point_format" in attrs else 7 # 认为 7 支持rgb
version = attrs.pop("version") if "version" in attrs else "1.4" # 默认 1.4
header = laspy.LasHeader(point_format=point_format, version=version)
''' 自动计算 scales 和 offsets，确保坐标精度无损 '''
# /questions/77308057/conversion-accuracy-issues-of-e57-to-las-in-python-using-pye57-and-laspy
header.scale = attrs.pop("scales") if "scales" in attrs else [0.001, 0.001, 0.001] # 0.001 是毫米精度
# 如果保存时与原las文件的offsets不同，可能会出现偏移（偏移精度 < scales精度限制）
header.offset = attrs.pop("offsets") if "offsets" in attrs else np.floor(np.min(xyz, axis=0))
''' 初始化一些需要保存的属性值。 '''
extra_attr = [] # 如果是额外属性，添加到 header 中, 后续赋值
# 获取当前版本的保留字段（与 point_format 和 version 有关）
all_reserved_fields = get_las_header_attrs(point_format, version)
for attr, value in attrs.items():
if attr not in all_reserved_fields: # 添加额外属性，在 las 初始化后赋值，如 label, pred
header.add_extra_dim(laspy.ExtraBytesParams(name=attr, type=np.float32))
extra_attr.append(attr)
''' 创建 las 文件对象 '''
las = laspy.LasData(header)
# 添加xyz坐标
las.x = xyz[:, 0]
las.y = xyz[:, 1]
las.z = xyz[:, 2]
''' 添加颜色信息 '''
# 如果RGB全是0, 则不存储颜色。如果是归一化的颜色，则需要乘以 65535，转为 uint16
if (rgb is not None) and (np.max(rgb) > 0):
if np.max(rgb) <= 1:
rgb = (rgb * 65535).astype(np.uint16) # 65535 = 2^16 - 1, las存储颜色是16位无符号整型
elif np.max(rgb) <= 255:
rgb = (rgb / 255 * 65535).astype(np.uint16)
las.red = rgb[:, 0]
las.green = rgb[:, 1]
las.blue = rgb[:, 2]
''' 设置固有属性 '''
for attr, value in attrs.items():
if attr in all_reserved_fields: # 如果是保留字段，则不添加额外属性。如 X, Y, Z, Red, Green, Blue
if value.ndim == 2 and value.shape[1] == 1:
value = value.flatten()
las[attr] = value
''' 设置额外属性 '''
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.write(out_file)
def get_las_info(las_file):
"""
读取 LAS 文件的格式和版本信息（适用于 laspy ）
:param las_file: LAS 文件路径
:return: 版本号 () 和 点格式 ID
"""
with laspy.open(las_file) as las: # 以只读方式打开 LAS 文件
header = las.header # 读取头部信息
# 获取 LAS 版本号
version_major = header.version.major
version_minor = header.version.minor
version = f"{version_major}.{version_minor}"
# 获取点格式（point format ID）
point_format = header.point_format.id
return version, point_format
def get_las_header_attrs(point_format=7, version="1.4"):
"""
根据 point_format 和 version 获取 las 文件的 header 属性
说明文档：/en/latest/#point-records
Args:
point_format: 点格式
version: 版本
Returns:
"""
dimensions = []
header = laspy.LasHeader(point_format=point_format, version=version) # 7 支持rgb
for dim in header.point_format.dimensions:
dimensions.append(dim.name)
return dimensions
if __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 = "/home/gisleung/dataset2/Epoch_March2018/LiDAR/Mar18_test_small.las"
xyz_, rgb_, attrs_ = read_las_fit(read_las_path, ["scales", "offsets", "leung"])
print(attrs_)
# 测试3： 写入LAS文件
# save_las_path = "/home/gisleung/dataset2/Epoch_March2018/LiDAR/Mar18_test_small2_fit.las"
# write_las_fit(save_las_path, xyz_, rgb_, {
# # "scales": attrs_["scales"],
# # "offsets": attrs_["offsets"]
# })