在使用激光雷達(dá)進(jìn)行樣地調(diào)查時，我們常常會面臨一個問題：
明明只掃描了一個樣地，導(dǎo)出的點云數(shù)據(jù)卻覆蓋了周圍一大片區(qū)域。
如下圖所示：
為什么會這樣？因為激光雷達(dá)的掃描距離通常很遠(yuǎn)，即使是手持設(shè)備，掃描范圍也輕松超過幾十米，有些甚至可達(dá)200米。這意味著你掃描一個小樣地，結(jié)果卻得到了一大片“用不到”的數(shù)據(jù)。
這些多余的數(shù)據(jù)不僅會占用大量存儲空間，還會讓后續(xù)的樹木提取、單木結(jié)構(gòu)參數(shù)提取等處理過程變得又慢又麻煩。
所以，我們第一步要做的，就是： 把只屬于樣地范圍的點云“裁剪”出來。
無論你是機載數(shù)據(jù)還是手持?jǐn)?shù)據(jù)，處理的范圍是樣地還是林分，首先第一步就是根據(jù)你的研究區(qū)對采集數(shù)據(jù)進(jìn)行裁剪。
在這篇推文中，我將介紹兩種常用的裁剪方法：
- 根據(jù)矩形坐標(biāo)范圍裁剪
- 根據(jù)樣地矢量范圍裁剪
此外，我還會提醒大家一個實用經(jīng)驗：
?? 裁剪范圍最好比實際樣地大5米左右，這樣可以保留邊界附近樹木的完整冠幅，若樣地具有較大的坡度嚴(yán)格按樣地邊界裁剪，甚至?xí)堰吔缒镜臉淠局鞲汕械簟?/p>
一、根據(jù)矩形坐標(biāo)范圍裁剪
根據(jù)矩形四至進(jìn)行裁剪，即根據(jù)最大和最小的 x、y 坐標(biāo)進(jìn)行限制，對點云數(shù)據(jù)進(jìn)行裁剪。
1.先導(dǎo)入包和點云數(shù)據(jù)
'''
@微信公眾號：生態(tài)遙感學(xué)習(xí)記錄
@作者：大馬哈余
按照四至范圍采集點云
'''
import laspy
import numpy as np
import pyvista as pv
file=r"C:\Users\YSYmc\Desktop\近期內(nèi)容\lidar_data\sample_l.las"
output_las = r'C:\Users\YSYmc\Desktop\近期內(nèi)容\lidar_data\clipped_small.las'
las=laspy.read(file)
2.根據(jù)四至對點數(shù)據(jù)進(jìn)行篩選（切片）
# 裁剪范圍
x_min, x_max = 534450, 534465
y_min, y_max = 4689430, 4689445
# 提取點的 X、Y 坐標(biāo)
x = las.x
y = las.y
# 先進(jìn)行bool篩選，篩選位于范圍內(nèi)的點
#這一步類似于dataframe的bool索引
mask = (x >= x_min) & (x <= x_max) & (y >= y_min) & (y <= y_max)
clipped_points = las.points[mask]
3.保存裁剪后的las點云數(shù)據(jù)
# # 創(chuàng)建新的 las 文件并寫入裁剪后的點
header = las.header
clipped = laspy.LasData(header)#寫入頭文件信息
clipped.points = clipped_points#寫入點數(shù)據(jù)到clipped中
# 保存之前可先進(jìn)行可視化
clipped.write(output_las)
4.可視化點云數(shù)據(jù)
points = np.vstack((clipped.points.x, clipped.points.y, clipped.points.z)).T
cloud = pv.PolyData(points)
cloud["Elevation"]=clipped_points.z
color_field = "Elevation" # 你可以改成其他字段
cloud.plot(
scalars=color_field #基于該屬性字段進(jìn)行著色
,cmap="viridis" # 其他的色帶：'viridis', 'plasma', 'terrain', 'coolwarm'
,point_size=2 #控制點的顯示大小
,show_scalar_bar=True #顯示colorbar
,render_points_as_spheres=True)#該參數(shù)用于控制點云在渲染時是否以球體的形式顯示，是的話顯示效果更好
二、根據(jù)樣地矢量范圍裁剪
理解上面的四至裁剪點云數(shù)據(jù)，下面的自然也很好理解。重要的就是判斷點是不是在矢量（多邊形）范圍內(nèi)，這里我使用的是shapely.contains_xy，判斷的比較快。
import laspy
import geopandas as gpd
import shapely
file=r"C:\Users\YSYmc\Desktop\近期內(nèi)容\lidar_data\sample_l.las"
shp=r"C:\Users\YSYmc\Desktop\近期內(nèi)容\lidar_data\clip_polygon01.shp"
output_las = r'C:\Users\YSYmc\Desktop\近期內(nèi)容\lidar_data\clipped_shp.las'
las=laspy.read(file)
x = las.x
y = las.y
# 讀取裁剪矢量（Shapefile）
clip_gdf = gpd.read_file(shp)
clip_polygon = clip_gdf.geometry.iloc[0] # 假設(shè)只有一個面
# 判斷每個點是否在多邊形內(nèi)
# 使用 shapely.contains_xy 判斷每個點是否在多邊形內(nèi)，這個函數(shù)判斷速度比較快
mask = shapely.contains_xy(clip_polygon, x, y)
# 裁剪點
clipped_points = las.points[mask]
# 寫出裁剪后的 las 文件
clipped_las = laspy.LasData(las.header)
clipped_las.points = clipped_points
clipped_las.write(output_las)
