使用Python处理Las文件，提取激光雷达离散点云数据信息（自动选取植被点云样本）

最近获取了一批激光雷达离散点云数据，但是不知如何下手处理，于是查阅资料，发现python的laspy库可以直接处理，于是对此展开了一定的研究。

由于点云数量众多，如果要利用点云进行一些计算或者标注，数据量很庞大，费时费力。

本文主要目的是提取出离散点云数据中的植被区点云，即自动选取植被点云样本，避免手动标注的过程。这样的处理过程直接利用点云信息，比直接利用软件分类来的快得多。

话不多说，直接上码：

1. Las文件处理，点云信息提取

首先导入laspy和numpy库

import laspy
import numpy as np

接着读取存放激光雷达点云数据的las文件

las_file="las0.las"
las=laspy.read(las_file)

查看头文件，这里显示点云格式是7，不懂的同学可以查看laspy的文档，7号格式有介绍

#获取文件头
header=las.header

<LasHeader(1.4, <PointFormat(7, 0 bytes of extra dims)>)>

接着我们查看一下点云的属性字段名

#属性字段名
dimension_names=point_format.dimension_names
all_elements = list(dimension_names)
print(all_elements)

可以看到，有经纬度高程（XYZ），回波强度，回波数量，RGB等等；

接下来查看点云个数：

#点个数
point_num=header.point_count
print(f'总的点云个数为:{point_num}个')

可以看到，2.5e+8个点云，数量也是非常庞大，因此下面我们就要从这里面提取一部分植被信息了；

2. 自动化提取植被点云样本

这里利用VDVI来提取植被点云，VDVI是NDVI的变体，由于这里激光雷达数据仅含有RGB三个波段信息，没有近红外信息，因此用VDVI来代替NDVI，计算方法如下：

las_x=np.array(las.x)
las_y=np.array(las.y)
las_z=np.array(las.z)
las_r=np.array(las.red)
las_g=np.array(las.green)
las_b=np.array(las.blue)
np.seterr(divide='ignore',invalid='ignore')
# 计算VDVI
def VDVI_com(R,G,B):
    R = np.array(R, dtype = float)
    G = np.array(G, dtype = float)
    B = np.array(B, dtype = float)
    VDVI = (2*G-R-B) / (2*G+R+B)
    return VDVI

VDVI = VDVI_com(las_r,las_g,las_b)

然后自行判断阈值，这里选取VDVI为0.1，即VDVI大于等于0.1的为植被点云：

veg_ind = (VDVI >= 0.1) # 根据阈值找到VDVI大于等于0.1的植被区的索引

然后我们看下效果，由于数量巨大，因此只查看前1000个点云（veg_indices）：

veg_indices = np.where(VDVI >= 0.1)[0]
# 组合,提取植被区的1000个color,veg,用作绘图
selected_indices = veg_indices[:1000]
veg = np.stack([las_x[selected_indices], las_y[selected_indices], las_z[selected_indices]], axis=1)
colors = np.stack([las_r[selected_indices], las_g[selected_indices], las_b[selected_indices]], axis=1)
veg.shape

绘图，查看效果：

# 绘图
import matplotlib.pyplot as plt

# 创建一个三维图形对象
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# 绘制散点图
ax.scatter(veg[:,0], veg[:,1], veg[:,2],c=colors/255,marker='.')

# 设置坐标轴标签
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

# 显示图形
plt.show()

可以看到，图中点云均为植被，说明效果还是不错的。

3. 保存提取后的植被点云

这一步是保存原有的属性名，如RGB,XYZ这种：

# 这一步用于取出所有元素，便于下面保存
# 创建一个字典来存储变量
las_data = {}

# 循环遍历元素名，生成相应的变量
for element_name in all_elements:
    # 使用getattr函数根据元素名从las对象中获取对应的属性
    element_data = getattr(las, element_name.lower(), None)
    
    # 如果属性存在，将其转换为NumPy数组并存储在字典中
    if element_data is not None:
        las_data[element_name.lower()] = np.array(element_data)

这一步是保存las文件，首先创建一个point_record来存储点云，然后设置一个循环，关键为data_to_assign = las_data[i_lower][veg_ind]这一步，即将植被区的点云数据提取到data_to_assign中，然后利用setattr函数，将提取的数据信息赋值到point_record，最后利用writer导出，las1.las即为保存好的植被区点云。

# 保存las文件,只保存植被区域
with laspy.open("las1.las", mode="w", header=header) as writer:
    point_record = laspy.ScaleAwarePointRecord.zeros(num, header=header)
    for i in all_elements:
        # 将属性名称转换为小写，以匹配字典中的键
        i_lower = i.lower()
        data_to_assign = las_data[i_lower][veg_ind]
        # 使用setattr函数将las_data中的值赋给point_record对象的属性
        setattr(point_record, i_lower, data_to_assign)
    
    writer.write_points(point_record)