Python处理LiDAR点云：用Laspy库从读取到可视化的保姆级实战（附完整代码）
当你第一次拿到一个LAS格式的LiDAR点云文件时，可能会感到无从下手。这些密密麻麻的三维坐标点背后藏着丰富的地理信息，但如何用Python打开它们、理解它们，最终让这些数据"说话"？本文将带你从零开始，使用laspy这个强大的Python库，一步步完成从数据读取到3D可视化的完整流程。
1. 环境准备与数据获取
在开始之前，我们需要确保Python环境已经配置好所有必要的工具。推荐使用Python 3.8或更高版本，并创建一个干净的虚拟环境来管理依赖。
安装核心库非常简单，只需运行以下命令：
pip install laspy numpy matplotlib
注意：如果你计划处理大型点云文件（超过1GB），建议额外安装
laspy[laz]
以支持LAZ压缩格式：
pip install laspy[laz]
关于测试数据，可以从以下几个公开资源获取：
- OpenTopography ：提供全球多个地区的免费LiDAR数据
- USGS 3DEP ：美国地质调查局的高程数据计划
- 当地政府开放数据平台 ：许多城市会公开本地的LiDAR扫描数据
为了本教程的连贯性，我们假设你已经获得了一个名为
sample.las
的点云文件。这个文件可能包含数十万到数百万个空间点，每个点都有x、y、z坐标以及可能的附加属性（如强度、分类、颜色等）。
2. 理解LAS文件结构
在直接操作数据之前，了解LAS文件的内部结构至关重要。LAS格式是LiDAR数据的行业标准，它不仅仅包含原始点坐标，还有丰富的元数据信息。
一个典型的LAS文件由两部分组成：
-
文件头(Header) ：包含全局元数据
- 文件版本信息
- 坐标参考系统
- 点的数量
- 空间范围（最小/最大坐标）
- 缩放因子和偏移量（用于优化存储）
-
点记录(Point Records) ：实际的三维点数据
- 必选：X,Y,Z坐标
- 可选：强度(Intensity)、回波次数(Return Number)、分类(Classification)、颜色(RGB)等
让我们用laspy打开文件并查看这些信息：
import laspy
# 读取LAS文件
las = laspy.read("sample.las")
# 访问文件头信息
header = las.header
print(f"文件版本: {header.major_version}.{header.minor_version}")
print(f"点数量: {header.point_count}")
print(f"空间范围: X({header.min[0]}, {header.max[0]}), Y({header.min[1]}, {header.max[1]}), Z({header.min[2]}, {header.max[2]})")
print(f"缩放因子: {header.scale}")
print(f"偏移量: {header.offset}")
关键概念解释 ：
-
缩放因子(Scale)
：为了节省存储空间，LAS文件中的坐标值实际上是整数，需要通过
实际坐标 = 存储值 × 缩放因子 + 偏移量的公式转换 -
点格式(Point Format)
：决定了每个点包含哪些属性，常见的有：
- 格式0：仅包含位置信息
- 格式1：包含位置和GPS时间
- 格式2：包含位置和RGB颜色
- 格式3：包含位置、GPS时间和RGB颜色
3. 深入探索点云数据
现在我们已经了解了文件的整体结构，接下来让我们深入查看具体的点数据。laspy提供了多种访问点属性的方式，最直接的是通过点对象的属性访问器。
# 访问前5个点的坐标
print("X坐标示例:", las.x[:5])
print("Y坐标示例:", las.y[:5])
print("Z坐标示例:", las.z[:5])
# 检查点云是否包含额外属性
if hasattr(las, 'intensity'):
print("强度值范围:", las.intensity.min(), "到", las.intensity.max())
if hasattr(las, 'classification'):
print("分类代码:", set(las.classification))
if hasattr(las, 'rgb'):
print("RGB颜色通道可用")
常见问题排查 ：
-
如果遇到
AttributeError提示某个属性不存在，说明你的点格式不支持该属性 - 坐标值看起来非常大或非常小？记得它们可能还应用了缩放和偏移
-
内存不足？考虑使用laspy的
chunked_read功能分批处理大型文件
为了更好地理解数据分布，我们可以使用numpy进行一些统计分析：
import numpy as np
# 将坐标转换为numpy数组
coords = np.vstack((las.x, las.y, las.z)).transpose()
# 计算基本统计量
print("高程(Z)统计:")
print(f" - 平均值: {np.mean(las.z):.2f}")
print(f" - 标准差: {np.std(las.z):.2f}")
print(f" - 中位数: {np.median(las.z):.2f}")
print(f" - 95%分位数: {np.percentile(las.z, 95):.2f}")
4. 数据预处理技巧
原始LiDAR数据通常需要一些预处理才能用于分析或可视化。下面介绍几种常见的操作：
4.1 坐标转换
由于LAS文件使用相对坐标和缩放存储，有时我们需要将其转换为绝对坐标：
def scale_and_offset(coords, scales, offsets):
"""将LAS内部坐标转换为真实世界坐标"""
return coords * scales + offsets
# 转换所有点
real_coords = scale_and_offset(
coords,
np.array([header.scale[0], header.scale[1], header.scale[2]]),
np.array([header.offset[0], header.offset[1], header.offset[2]])
)
4.2 点云过滤
根据特定条件筛选点可以大幅减少数据量，提高处理效率：
# 示例1：基于高程过滤
ground_points = las.z < (np.median(las.z) + 1)
# 示例2：基于强度过滤
high_intensity = las.intensity > np.percentile(las.intensity, 75)
# 示例3：基于分类代码过滤（2通常代表地面）
ground_class = las.classification == 2
# 组合条件
filtered_points = ground_class & (las.z < np.percentile(las.z, 50))
4.3 降采样处理
对于超高密度的点云，适当的降采样可以提高处理速度：
# 随机选择10%的点
random_indices = np.random.choice(len(las.points), size=len(las.points)//10, replace=False)
downsampled = las.points[random_indices]
5. 3D可视化实战
理解了数据结构并完成必要的预处理后，现在是时候将点云可视化了。我们将使用matplotlib的3D功能来创建交互式视图。
5.1 基础可视化
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 绘制点云 - 使用高程值作为颜色映射
sc = ax.scatter(
las.x[::100], # 每100个点取一个，避免过度绘制
las.y[::100],
las.z[::100],
c=las.z[::100],
cmap='viridis',
s=0.1,
alpha=0.5
)
# 添加颜色条
plt.colorbar(sc, ax=ax, label='高程值')
# 设置标签和标题
ax.set_xlabel('X坐标')
ax.set_ylabel('Y坐标')
ax.set_zlabel('高程')
ax.set_title('LiDAR点云3D视图')
plt.tight_layout()
plt.show()
5.2 高级可视化技巧
基础散点图虽然直观，但有时我们需要更专业的可视化方式：
按分类着色 ：
# 创建分类颜色映射
class_colors = {
0: [0.5, 0.5, 0.5], # 未分类 - 灰色
1: [0.9, 0.9, 0.9], # 未分配 - 浅灰
2: [0.2, 0.6, 0.2], # 地面 - 绿色
3: [0.8, 0.3, 0.3], # 低植被 - 红色
4: [0.4, 0.8, 0.4], # 中植被 - 浅绿
5: [0.1, 0.4, 0.1], # 高植被 - 深绿
6: [0.7, 0.7, 0.1], # 建筑 - 黄色
}
# 为每个点分配颜色
colors = np.array([class_colors.get(cls, [1,0,0]) for cls in las.classification])
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(las.x[::100], las.y[::100], las.z[::100], c=colors[::100], s=0.2)
ax.set_title('按分类着色的LiDAR点云')
plt.show()
交互式可视化 ：
虽然matplotlib的3D视图可以旋转缩放，但对于大型点云，更推荐使用专门的工具如
pyvista
或
open3d
：
# 需要先安装：pip install pyvista
import pyvista as pv
# 创建PyVista点云对象
cloud = pv.PolyData(np.vstack((las.x, las.y, las.z)).T)
cloud["elevation"] = las.z # 添加高程作为标量字段
# 创建绘图窗口
plotter = pv.Plotter()
plotter.add_mesh(cloud, point_size=2, scalars="elevation")
plotter.show()
6. 实用技巧与性能优化
处理大型LiDAR数据集时，性能和内存管理变得至关重要。以下是几个实用建议：
6.1 内存高效处理
# 使用chunked_read处理大文件
with laspy.open('large.las') as reader:
for chunk in reader.chunk_iterator(2_000_000): # 每次处理200万个点
process_chunk(chunk) # 你的处理函数
6.2 并行处理
from multiprocessing import Pool
def process_points(points):
# 处理点数据的函数
return results
# 将点云分成4部分并行处理
with Pool(4) as p:
results = p.map(process_points, np.array_split(las.points, 4))
6.3 数据导出
处理后的数据可能需要导出为其他格式：
# 导出为CSV
np.savetxt('points.csv', np.vstack((las.x, las.y, las.z)).T, delimiter=',')
# 导出为新的LAS文件
with laspy.open('output.las', mode='w', header=las.header) as writer:
writer.write_points(las.points)
7. 实际应用案例
为了帮助你将所学应用到实际问题中，这里展示一个完整的LiDAR数据处理流程示例：
import laspy
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def process_lidar_file(input_path, output_path=None):
"""完整的LiDAR处理流程"""
# 1. 读取数据
las = laspy.read(input_path)
# 2. 提取地面点（分类代码2）
ground = las.classification == 2
ground_points = las.points[ground]
# 3. 计算地面高程模型
x = las.x[ground]
y = las.y[ground]
z = las.z[ground]
# 4. 可视化
fig = plt.figure(figsize=(12, 6))
# 原始点云
ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(las.x[::100], las.y[::100], las.z[::100],
c=las.z[::100], s=0.1, cmap='terrain')
ax1.set_title('原始点云')
# 地面点
ax2 = fig.add_subplot(122, projection='3d')
ax2.scatter(x[::10], y[::10], z[::10],
c=z[::10], s=1, cmap='terrain')
ax2.set_title('提取的地面点')
plt.tight_layout()
plt.show()
# 5. 可选：保存地面点
if output_path:
with laspy.open(output_path, mode='w', header=las.header) as writer:
writer.write_points(ground_points)
return ground_points
# 使用示例
ground = process_lidar_file('sample.las', 'ground_points.las')
扫一扫
8364
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
