如何基于 DEM 和 LCZ 反向求解 DSM

🎓 作者：Clarmy

📢 版权声明：此文的所有文字、图片、代码以及相应数据文件的由本人（和鲸社区：Clarmy）整理、撰写，文责自负。严禁任何商业性质的转载。公益性质转载需联系作者本人获取授权。转载本文时，请务必文字注明“来自：和鲸社区：Clarmy”，并附带本项目超链接:https://www.heywhale.com/mw/project/64f1be79abb16766700babe6。

以下全文代码均已发布至和鲸社区，复制下面链接，可一键fork跑通：https://www.heywhale.com/mw/project/64f1be79abb16766700babe6?shareby=620ce369c1ae5e001747085e

这一篇是一个粉丝的投稿万事屋的需求解答，也是被和鲸社区运营小姐姐选中的一篇约稿，所以优先处理了。

粉丝的原始诉求是用 DEM（Digital Elevation Model） 结合 LCZ（Local Climate Zones）数据，获得一个城市地区更高准确率的 DEM 数据。

根据我后续的一些调研，理解为用上述两个数据反算 DSM（Digital Surface Model） 数据，由于 LCZ 数据所分类的建筑物类型，是有高低的区分的，所以这个需求从理论上是完全可行的。

首先可能需要先解释一下 LCZ 是个什么东西，在我看到这个需求之前我也不知道 LCZ 是什么，根据 claude.ai 的解释：

LCZ 是指地表覆盖区(Land cover zones)或地表覆盖类型(Land cover types)的缩写,主要用于遥感图像分类和地表覆盖研究。

它表示对地表覆盖进行分类和划分不同类型的区域,这在利用遥感影像进行地表覆盖提取和信息提取中是非常重要的一个步骤。常见的LCZ分类有:

• 水体 - 湖泊、河流等
• 植被 - 草地、农田、森林等
• 城市 - 商业区、住宅区、工业区等
• 裸地 - 裸露的土地、岩石等

遥感图像的数字化分类一般会区分出这些典型的LCZ类型,这有助于对 study area 的地表覆盖情况有一个整体概述和理解,并可以用于相关的生态、气候、城市化研究。LCZ的划分也是遥感图像分类的一个重要环节。

后来拿到了粉丝分享的南京市 LCZ 的 GeoTiff 数据以后，直观地认识了这种数据，通俗地说它其实就是根据预设的规则，把城市地区的地貌分为17种类别，每一种类别用一个数字表示，17种类别如下图所示。

图片来源：https://doi.org/10.1175/BAMS-D-11-00019.1

而根据粉丝的需求，我们其实是从这17种分类中，挑选出表示建筑的分类，然后根据它们的高低特征来修正 DEM 数据。事实上这些建筑类型的地貌是没有准确的高度值的，这就需要我们根据一些先验知识对它的高度先做一个预设。

也是根据 claude.ai 阅读了 Local Climatic Zoning and Urban Heat Island in Beirut（DOI:10.1016/j.proeng.2016.10.026） 这篇论文以后，对各种建筑类型地貌的高度进行了估算：

LCZ中表示建筑物的类别主要有:
Compact high-rise (LCZ 1): 紧凑型高层建筑,建筑物高度>25米
Compact midrise (LCZ 2):紧凑型中层建筑,建筑物高度10-25米
Compact low-rise (LCZ 3):紧凑型低层建筑,建筑物高度3-10米
Open high-rise (LCZ 4):开敞型高层建筑,建筑物高度>25米
Open midrise (LCZ 5):开敞型中层建筑,建筑物高度10-25米
Open low-rise (LCZ 6):开敞型低层建筑,建筑物高度3-10米

LCZ中的建筑物高度分类是有明确的量化标准的:
低层建筑:3-10米
中层建筑:10-25米
高层建筑:>25米
这些标准高度阈值有助于根据建筑物高度划分不同的LCZ类型。

那么我们就可以根据这种先验知识，以一种近似的方式解决这个问题，方法就是先按照不同类别把预设的建筑高度赋值给 LCZ 矩阵，形成一个相对地面凸出的高度矩阵，然后把它与 DEM 矩阵进行坐标对齐，最后把二者相加，就得到了被建筑高度近似修正的 DEM 数据（后来查了一下，这种数据其实就是 DSM，即数字地表面模型，它是包含了土地及其上覆的建筑、植被高度的一类地理数据）。

好了，先验知识大概介绍完了，现在介绍如何使用 Python 解决这个问题。

我们需要用的包如下：

numpy
scipy
matplotlib
rasterio==1.3.8

这里需要强调一下 rasterio 包的版本，因为求解过程中涉及到该包的一些函数在不同版本下的行为是不一样的。

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

我们先来设置一个常量，用于设置不同建筑分类的高度范围：

# 不同类型建筑的最高最低高度范围设置
RAISES = {
    1: (25, 50),   # 25m ~ 50m 之间
    2: (10, 25),   # 10m ~ 25m 之间
    3: (3, 10),    # 以此类推
    4: (25, 50),
    5: (10, 25),
    6: (3, 10),
    7: (3, 10),
    8: (3, 10),
    9: (3, 10)
}

此外，由于不同分类还涉及到建筑物分布的特征，例如“密集”或“稀疏”的描述，因此我们还应该增加一个密度参数用于对建筑物分布密集程度做调节。

# 不同类型建筑的分布密度值
DENSITY = {
    1: 0.8,  # 紧密分布
    2: 0.8,
    3: 0.8,
    4: 0.5,  # 相对松散
    5: 0.5,
    6: 0.5,
    7: 0.8,
    8: 0.3,  # 稀疏分布
    9: 0.3
}

然后，我们可以先开始加载 DEM 数据及其对应的经纬度网格：

# 读取 DEM 及相应的经纬网格数据
dem = np.load('dem.npy')
dem_lons = np.load('lons.npy')
dem_lats = np.load('lats.npy')
dem[dem<-5] = -5  # 将小于-5 的值替换为 -5

上面这段所使用的 dem 数据是使用 pyterrain 包提前下载得到的，经纬度变量也是二维网格形式。

先浅画一下：

plt.contourf(dem_lons, dem_lats, dem, levels=np.linspace(-5, dem.max(), 50), cmap=plt.cm.terrain)

熟悉南京地形的朋友应该可以认出来。

下面，我们来读取和解析 LCZ 数据

dataset = rasterio.open('./nj.tif')

lcz = dataset.read(1)

transform = dataset.transform

nx, ny = dataset.width, dataset.height

x_pixel = np.arange(nx)
y_pixel = np.arange(ny)

x_coords, y_coords = rasterio.transform.xy(transform, y_pixel[:, np.newaxis], x_pixel)  # 通过仿射参数计算经纬网格坐标

lcz_lons = np.stack(x_coords)
lcz_lats = np.stack(y_coords)

这里得到的 lcz 就是原始数据的 LCZ 矩阵，它是由 0-17 的数字组成的二维矩阵，lcz_lons 和 lcz_lats 是其对应的二维经纬网格。

直接画一张图来大概长这样（使用离散型 colormap）：

plt.contourf(lcz_lons, lcz_lats, lcz, cmap=plt.cm.tab20)

得到带有经纬坐标的 LCZ 矩阵以后，还需要将其向 DEM 数据的坐标对齐才能进行后续的融合操作，所以我们现在对 LCZ 矩阵进行一次临近点插值。

# 将 lcz 数据按照 dem 数据的坐标插值
lcz1d = lcz.flatten()
lczlons1d = lcz_lons.flatten()
lczlats1d = lcz_lats.flatten()

coords = np.stack([lczlons1d, lczlats1d]).T

new_lcz = griddata(coords, lcz1d, (dem_lons, dem_lats), method='nearest')  # 使用最临近插值
assert new_lcz.shape == dem.shape

插值完画图检查一下：

plt.contourf(dem_lons, dem_lats, new_lcz, cmap=plt.cm.tab20)

好的，看起来与之前的图并没有太大区别，说明可以排除插值过程出现错误的情况。现在的 new_lcz 虽然跟 dem 矩阵的形状相同了，但是 new_lcz 仍然是分类值而非高度值，所以我们需要对其进行一次转换，根据之前预设RAISES和 DENSITY 常量生成高度偏移矩阵，在这里我们需要先定义一个生成随机零值索引的函数，用于调整不同类型的地貌分布密度：

from typing import Tuple, Any

def make_random_zindices(array: np.ndarray, density: float) -> Tuple[Any, ...]:
    """
    根据给定的密度生成随机索引，用于将输入数组中对应位置的元素替换为0。

    参数:
    array (np.ndarray): 输入的numpy数组。
    density (float): 保持非零元素的密度，范围应在0到1之间，
                     其中1表示全部保留，0表示全部替换为0。

    返回:
    Tuple[Any, ...]: 一个元组，包含了用于替换原数组中元素为0的索引。

    示例:
    >>> arr = np.array([[1, 2], [3, 4]])
    >>> indices = make_random_zindices(arr, 0.5)
    >>> arr[indices] = 0
    """
    # 计算数组中元素的总数量
    num_elements = array.size
    
    # 根据密度计算要替换为0的元素数量
    num_zero_elements = int(num_elements * (1 - density))

    # 生成要替换为0的元素的随机索引
    # np.unravel_index 函数将一维的索引转换为多维索引，以匹配原数组的形状
    zindices = np.unravel_index(np.random.choice(num_elements, num_zero_elements, replace=False), array.shape)
    
    return zindices

然后，我们再开始生成高度偏移矩阵：

lcz_diff = np.full(new_lcz.shape, 0)  # 创建建筑相对地面的偏移程度

# 计算建筑高度偏移矩阵
for i, (imin, imax) in RAISES.items():
    single_kind = new_lcz[new_lcz==i]
    zindices = make_random_zindices(single_kind, DENSITY[i])
    single_kind = np.random.randint(imin,imax,single_kind.shape)
    single_kind[zindices] = 0
    lcz_diff[new_lcz==i] = single_kind

把得到的结果浅画一下：

plt.contourf(dem_lons, dem_lats, lcz_diff)

完成这一步，基本上就算完成了，最后我们只需要将 dem 与 lcz_diff 相加就可以得到近似的 dsm 了。

# 计算 DSM
dsm = dem + lcz_diff

浅画一下：

plt.contourf(dem_lons, dem_lats, dsm, levels=np.linspace(-5, dsm.max(), 50), cmap=plt.cm.terrain)

这样看与前面直接画的 dem 已经有了明显的区别，主要是城市居民区的 dem 分布颜色有明显的变化。

现在我们可以筛选城中比较小的区域“放大”来看，这里选的大致是栖霞区范围，囊括了钟山风景区和栖霞山以及长江部分区域，具体的经纬范围是：

118.80156,32.184035   # 左上
119.010254,32.004412  # 右下

计算选取范围坐标索引的代码在这里就暂时略过，我们直接看画出的结果：

plt.contourf(dem_lons[idx], dem_lats[idx], dsm[idx], levels=np.linspace(-5, dsm.max(), 50), cmap=plt.cm.terrain)

可以看到居民区部分有很多蓝色的“斑点”，这些就是基于 LCZ 估算的建筑高度叠加的效果。

我们再来对比一下原始 DEM 和估算 DSM 的 3D 效果图：

最后，得到的 dsm 可以保存为 npy 文件:

np.save('dsm.npy', dsm)

以上。