利用las格式点云数据生成DSM影像

本方法使用las格式的点云数据生成DSM。

首先来了解DSM是什么。DSM是地学中常用到的数字表面模型。要从点云生成数字表面模型 (DSM)，首先要将感兴趣的点云区域划分为规则网格。随后，为每个网格单元分配一个高程值，该高程值源自单元内或附近点的3D坐标（X、Y、Z）。详细介绍查看博文从点云创建 DSM：网格化和可视化实用指南_点云网格化-优快云博客。

需要导入的库为

from scipy.interpolate import griddata
import laspy

其中，SciPy的interpolate模块提供了许多对数据进行插值运算的函数，其中griddata()函数可以在二维平面进行插值操作。其插值方法有三中，分别为“linear“（表示线性插值）、”nearest“（表示最邻近插值）、”cubic“表示三次样条插值。

laspy是一个Python库，主要用于读取、修改和创建LAS点云文件。本代码仅使用该库进行las格式点云的读取，其他操作参照官方用户手册Installation — laspy 2.5.0 documentation。

本文所使用到的数据为ISPRS 提供的Vauhingen数据集中的area1区域点云数据。 

生成DSM影像代码如下：

import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import laspy

path = r"las点云数据路径"
las = laspy.read(path)

points = las.xyz

x = points[:, 0]
y = points[:, 1]
z = points[:, 2]

# Define grid parameters
xmin, xmax = np.nanmin(x), np.nanmax(x)
ymin, ymax = np.nanmin(y), np.nanmax(y)

grid_resolution =1 #格网分辨率
 
# Create a grid
xi, yi = np.meshgrid(np.arange(xmin, xmax, grid_resolution), np.arange(ymin, ymax, grid_resolution))

# Interpolate elevation values
zi = griddata((x, y), z , (xi, yi), method='linear')#线性内插

plt.figure(figsize=(10, 8))
new_zi = np.nan_to_num(zi, nan=np.nanmax(zi))
new_zi1 = np.clip(new_zi, a_min=None, a_max=z.max()+1)
new_zi1 = np.flipud(new_zi1)#flipud使数组上下翻转

#显示
plt.imshow(new_zi1, extent=[xi.min(), xi.max(), yi.min(), yi.max()], cmap="terrain")
plt.colorbar(label='Elevation')
plt.title('Digital Elevation Model (DEM) with imshow')
plt.xlabel('X-coordinate')
plt.ylabel('Y-coordinate')
plt.show()

 生成结果：