4. 基于地形湿度的空间分布制图

任务拆解与数据准备

地形湿度指数（TWI）刻画的是“某像元汇集上游来水的能力与其排水能力之比”，经典定义为：

TWI = ln( a / tanβ )

其中 a 为单位等高线宽度的汇流面积（比汇流面积，m²/m），β 为坡度角（弧度）。a 越大、坡度越缓，TWI 越高，越容易形成湿润带。本例在 GIS 中完成三件事：由 DEM 派生坡度、由填洼 DEM 派生汇流累积量并换算为比汇流面积、两者相除取对数得到 TWI 栅格，最后与土地利用矢量做分区统计。

关键约束：a 与 β 必须来自同一投影坐标系（本例用 UTM，像元 30 m），否则单位不一致会让 TWI 失去物理意义。汇流累积量需先乘以像元边长再除以像元边长，等价于“累积像元数 × 30 m / 30 m”，即比汇流面积在栅格上可直接用“累积像元数”近似（单位宽度的假设下），但严谨做法是 a = flowacc × cellsize / cellsize，保留该步骤便于核查。

核心计算流程

下面给出可直接在 QGIS 图形化处理或 Python（rasterio + numpy）中复现的步骤。以 Python 为例，假设已有填洼后的 DEM 与 D8 汇流累积栅格：

import numpy as np
import rasterio

CELL = 30.0          # 像元边长，m
Z_MIN = 0.001        # 坡度下限，防止除零

with rasterio.open("dem_fill.tif") as src:
    dem = src.read(1).astype("float64")
    profile = src.profile.copy()
    nodata = src.nodata

with rasterio.open("flowacc_d8.tif") as src:
    facc = src.read(1).astype("float64")

# 1) 坡度（度）→ 弧度，用中心差分近似
gy, gx = np.gradient(dem, CELL)
slope_rad = np.arctan(np.hypot(gx, gy))
slope_rad = np.maximum(slope_rad, np.arctan(Z_MIN))

# 2) 比汇流面积 a = 累积像元数 × 像元边长 / 像元边长
a = facc * CELL / CELL
a = np.maximum(a, 1.0)

# 3) TWI
twi = np.log(a / np.tan(slope_rad))

# 4) 掩膜无效值
mask = (dem == nodata) | (facc <= 0)
twi[mask] = np.nan

profile.update(dtype="float32", nodata=np.nan, count=1)
with rasterio.open("twi.tif", "w", **profile) as dst:
    dst.write(twi.astype("float32"), 1)

计算完成后先做一次合理性检查：TWI 理论范围通常在 3～20 之间，若出现大面积负值或超过 30，多半是坡度下限过小或汇流累积单位错位。可用 np.nanpercentile(twi, [2, 50, 98]) 快速查看分布，再决定分级断点。

湿度分级与土地利用叠加

按分位数分级比等间距更稳健，因为 TWI 分布常右偏。示例采用 5 级：低湿（<P20）、较干（P20–P40）、中等（P40–P60）、较湿（P60–P80）、高湿（>P80）。分级后与土地利用矢量做分区统计，得到每类用地的平均 TWI 与高湿面积占比，用于判断“哪些用地类型偏向分布在湿润地形部位”。

表中数值为教学示例，实际结果取决于研究区 DEM 与用地分类。判读时注意：高湿占比高不等于“该用地造成湿润”，而是“该用地倾向占据湿润地形位”，二者因果方向不能混淆。

结果验证与常见问题

• 量纲核对：若把汇流累积像元数直接当作 a 而不乘像元边长，TWI 会整体偏移一个常数 ln(30)≈3.4，分级断点随之失效。本例保留 × CELL / CELL 写法，便于一眼看出单位处理。
• 坡度下限：平坦区 tanβ→0 会让 TWI 发散。设 0.001 相当于约 0.057°，可抑制极端值，但会人为抬高平地 TWI，需在报告中说明。
• 叠加一致性：土地利用矢量必须先投影到与 TWI 相同的 UTM 坐标系，再做分区统计，否则面积占比会因投影差异产生偏差。
• 验证方法：抽取若干高湿像元，回到 DEM 检查是否位于河谷或坡脚；再抽取低湿像元检查是否位于山脊，二者应呈现明显的地形位置差异。