赞
踩
最近获取了一批激光雷达离散点云数据,但是不知如何下手处理,于是查阅资料,发现python的laspy库可以直接处理,于是对此展开了一定的研究。
由于点云数量众多,如果要利用点云进行一些计算或者标注,数据量很庞大,费时费力。
本文主要目的是提取出离散点云数据中的植被区点云,即自动选取植被点云样本,避免手动标注的过程。这样的处理过程直接利用点云信息,比直接利用软件分类来的快得多。
话不多说,直接上码:
首先导入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个点云,数量也是非常庞大,因此下面我们就要从这里面提取一部分植被信息了;
这里利用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)
![](https://csdnimg.cn/release/blogv2/dist/pc/img/newCodeMoreWhite.png)
然后自行判断阈值,这里选取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()
![](https://csdnimg.cn/release/blogv2/dist/pc/img/newCodeMoreWhite.png)
可以看到,图中点云均为植被,说明效果还是不错的。
这一步是保存原有的属性名,如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)
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。