从点云创建 DSM：网格化和可视化实用指南

今天我将向您展示如何从点云创建数字表面模型（DSM）。首先，我们将尝试了解 DSM 是什么，然后我们将进入讨论的更实际部分。

什么是 DSM？

DSM 是一个描述表面及其表面所有内容的模型。现在，为了更清楚地了解 DSM，您不仅对某些区域的地貌特征感兴趣，而且还对位于该表面的所有树木、灌木丛和大岩石感兴趣。然而，在 DTM（数字地形模型）中，只有地球表面，没有树木、岩石或灌木丛。这里我将创建一个DSM。要从点云生成数字表面模型 (DSM)，第一步是将感兴趣的区域划分为网格。随后，为每个网格单元分配一个高程值，该高程值源自单元内或附近点的 3D 坐标（X、Y、Z）。此过程可全面呈现地球表面，捕获地形、建筑物和植被等特征。为了更好地理解网格化的工作原理，我将创建一个合成数据并自行实现网格化。

让我们首先导入必要的库：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import ticker, cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.spatial import cKDTree

创建网格区域和一些样本点：

# grid range in x and y direction
x = np.arange(8)
y = np.arange(8)
# creating high number of points (70) for this grid area 
# with x, y, z coordinates
xp = np.random.random(size=70)*8 - 0.5
yp = np.random.random(size=70)*8 - 0.5
zp = np.random.randint(10, size = 70)
# grid coodinates
xc = np.arange(8) # x coordinates 
yc = np.arange(8) # y coordinates
grc = np.meshgrid(xc, yc) # grid coordinates

现在，下一步我们将使用最近邻插值来找到距离网格中心最近的点：

tree = cKDTree(np.c_[xp, yp])
dd, ii = tree.query(np.c_[grc[0].ravel(), grc[1].ravel()], k=1)
x_close = xp[ii]
y_close = yp[ii]

可视化以便更好地理解：

plt.figure(figsize = (12, 12))
plt.xticks(x-0.5)
plt.yticks(y-0.5)
plt.scatter(xp, yp, label='Observed Points')
plt.scatter(grc[0], grc[1], label='Grid Centers')
plt.grid(color='k', linestyle='-', linewidth=2)
plt.xlim((-.5,7.5))
plt.ylim((-.5,7.5))for i in range(len(ii)):plt.plot([grc[0].ravel()[i], xp[ii[i]]], [grc[1].ravel()[i], yp[ii[i]]], c='g')#plt.text(-0.11, 7.5-0.31, 'A )', bbox=dict(fill=False, edgecolor='r', linewidth=1.25), fontsize=18)plt.legend(bbox_to_anchor=(.95, 1.02), loc='upper left')

观察到蓝色点，橙色点是网格中心，绿线显示网格中心和距离该中心最近的点之间的连接。

然后，一旦为每个网格分配一个值，您就得到了栅格。让我们想象一下：

from matplotlib import ticker, cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
fig, ax = plt.subplots(1, figsize = (12, 12))
ax.set_xticks(x-0.5)
ax.set_yticks(y-0.5)
elv = zp[ii]
elv.shape = (x.shape[0], y.shape[0])
im = ax.imshow(np.flip(elv, 0), cmap=cm.get_cmap('PuBu_r', 9))
ax.grid(color='k', linestyle='-', linewidth=2)
ax.tick_params('both', bottom=False, top=False, left=False, right=False,labelbottom=False, labeltop=False, labelleft=False, labelright=False)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = fig.colorbar(im, ticks=np.arange(10), ax = ax, cax = cax)
cbar.ax.set_yticklabels(['0 m    ', '1 m    ', '2 m    ', '3 m     ', '4 m', '5 m', '6 m', '7 m', '8 m', '9 m'])  
cbar.ax.set_ylabel('Elevation of interpolated Points', rotation=270, fontsize=15)

上述合成点的栅格结果。

对现实世界数据进行网格化。

下一步，我将使用从 OpenTopography 下载的真实世界数据。可以从这里获取。让我们看看数据是什么样的：

来自 Colrado (OpenTopography) 的点云数据。

现在让我们对现实世界的点云数据进行网格化，但这次我将使用 scipy griddata 而不是编写自己的网格化操作：

import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
import laspy
path = "pointsColoradowithoutBuilding.las"
las = laspy.read(path)
points = las.xyz
amplitude = las.intensity
red = las.red
green = las.green
blue = las.blue
from scipy.interpolate import griddata
import matplotlib.pyplot as pltx = 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 = .5# 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_zi = np.clip(new_zi, a_min=None, a_max=z.max()+1)
plt.imshow(new_zi, 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()