第5章 GIS操作与DEM处理

gis.py 是OGGM中代码量较大的单文件模块之一，约2180行代码，承担着所有冰川地理空间处理的入口职责。从原始DEM（Digital Elevation Model）的读取与平滑，到冰川掩码的栅格化（rasterization）生成，再到网格参数的选择和DEM重投影，gis.py 将遥感与GIS领域的基础操作封装为OGGM工作流中的标准entity_task。本章按数据处理链路逐步拆解这些关键函数，揭示其内部算法和参数含义。

5.1 模块概览与数据流

在OGGM的完整工作流中，gis.py 处于最前端的位置。一条冰川的处理从GIS任务开始，其整体数据流为：

dem.tif + outlines.shp
    │
    ▼
[process_dem] ──── gridded_data.nc (topo, topo_smoothed, topo_valid_mask)
    │
    ▼
[glacier_masks] ──── gridded_data.nc (glacier_mask, glacier_ext)
    │
    ▼
[glacier_grid_params] ──── glacier_grid.json
    │
    ▼
[gridded_attributes] ──── gridded_data.nc (slope, aspect, dis_from_border, ...)
    │
    ▼
[compute_hypsometry_attributes] ──── hypsometry.csv
    │
    ▼
几何（outline）导出 ──── geometries.pkl

整个流程由 gis_prepro_tasks 函数（约第1850行）编排，它将各步骤注册为GlobalGCModel的任务图节点，确保在正确的工作流阶段执行。

进阶

理解GIS模块的关键

gis.py 中的几乎所有函数都操作两种核心数据结构：glacier_grid.json（标量参数，如地图投影、网格分辨率）和 gridded_data.nc（二维栅格数组，如高程场、掩码）。理解这两类数据的读写模式是阅读本章的先决条件。

5.2 GriddedNcdfFile 上下文管理器

在深入各处理函数之前，必须先了解 GriddedNcdfFile（约第682-763行），它是所有栅格数据的统一容器。这个类基于 NetCDF4 实现了一个上下文管理器，负责单个NetCDF文件中多个二维变量的读写。

class GriddedNcdfFile():
    """单个NetCDF文件中管理所有gridded_data变量的上下文管理器"""

    def __init__(self, gdir):
        self.f = None
        self.fid = None
        self.gdir = gdir

    def __enter__(self):
        # 以 'a' 模式打开NetCDF（不存在则创建）
        self.fid = netCDF4.Dataset(self.gdir.get_filepath('gridded_data'), 'a')
        return self

    def create_dim(self, nx, ny):  ...
    def create_var(self, name, data): ...
    def read_var(self, name, **kw): ...
    def __exit__(self, ...):  ...

关键设计决策：使用NetCDF而非多个独立GeoTIFF，是因为单个冰川的栅格数据维度通常不超过200x200像素，将高程、掩码、坡度、坡向等变量储存在一个文件中减少了I/O开销，并且NetCDF自带的元数据机制便于后续分析时的变量发现。

5.3 DEM处理工作流（process_dem）

process_dem（约第767-899行）是整个GIS预处理的核心函数。它接收原始DEM和冰川轮廓，生成平滑后的高程场，是整个模型的地形基础。按处理顺序可分为四个阶段。

5.3.1 第一阶段：读取GeoTIFF DEM

DEM的读取通过 read_geotiff_dem 函数委托给 rasterio 库完成。这一函数不仅读取高程数组，还同步提取地理参考元数据：

def read_geotiff_dem(gdir):
    # 使用rasterio读取DEM GeoTIFF
    with rasterio.open(dem_path) as ds:
        topo = ds.read(1)  # 读取第一波段
        # 提取仿射变换参数 (左上角坐标 + 像素分辨率)
        transform = ds.transform
    return topo, transform

topo 是一个二维float64数组，单位为米。当DEM中存在NaN或负值（表示海洋/空洞）时，后续的空洞填充算法将被触发。

5.3.2 第二阶段：空洞填充（Void Filling）

DEM空洞填充采用两级策略（约第820-845行）。首先对小于cfg.PARAMS['max_thin_void']（默认50像素）的小空洞使用双线性插值（scipy.interpolate.griddata 的 linear 模式）。对于更大的空洞，使用最近邻外推（nearest 模式）。如果空洞覆盖面积超过max_void_percentage（默认1.0），则整个DEM被拒绝并返回错误。

进阶

空洞填充的边界效应

当大空洞出现在冰川区域边缘时，最近邻外推可能引入虚假地形。OGGM在填充前先计算 topo_valid_mask（原始有效像素掩码），存于gridded_data.nc 中，供后续坡度计算时排除填充区域。

5.3.3 第三阶段：高斯平滑

平滑是DEM处理中最关键的步骤之一，通过 gaussian_blur 函数（约第91-115行）实现：

def gaussian_blur(grid, sigma):
    """使用FFT卷积进行高斯平滑"""
    # 1. 生成高斯核
    l = int(sigma * 4)
    x = np.arange(-l, l+1)
    kern = np.exp(-(x**2)/(2*sigma**2))
    kern /= kern.sum()

    # 2. 在频域执行二维分离卷积（先x后y，或先y后x）
    from scipy.signal import fftconvolve
    smoothed_x = fftconvolve(grid, kern[np.newaxis, :], mode='same')
    smoothed = fftconvolve(smoothed_x, kern[:, np.newaxis], mode='same')
    return smoothed

平滑窗口大小由 cfg.PARAMS['smooth_window'] 控制。这里的 sigma 被设置为 smooth_window / dx（其中dx是像素分辨率），默认值为3.5格网单位。FFT卷积相比空域卷积在处理大核时性能更优。平滑后的高程场写入topo_smoothed变量。

5.3.4 第四阶段：裁剪与输出

可选的裁剪参数 clip_dem0 将低于海平面（0 m a.s.l.）的像素设为0。最终输出三个变量到 gridded_data.nc：topo（原始高程）、topo_smoothed（平滑高程）、topo_valid_mask（有效像素掩码）。

参数名	默认值	作用
`smooth_window`	3.5	高斯平滑核宽度（像素单位）
`topo_interp`	`'linear'`	小空洞填充方法（linear/nearest/spline）
`max_thin_void`	50	"小型"空洞的最大像素数（超过则用nearest）
`max_void_percentage`	1.0	空洞面积占DEM的最大允许比例（%）
`clip_dem0`	`False`	是否将低于0m的像素强制设为0

5.4 冰川掩码计算（glacier_masks）

glacier_masks（约第902-1032行）将矢量形状文件（outlines.shp）中的冰川多边形转换为栅格掩码。这是连接矢量GIS数据和栅格模型的桥梁。

5.4.1 多边形到像素坐标的转换

转换分两步进行。首先 _interp_polygon 将冰川多边形的顶点等距插值到 grid_dx 的分辨率级别——这对于保持小冰川的几何精度至关重要，因为原始RGI多边形可能只有十几个顶点。然后 _polygon_to_pix 将地理坐标（经度/纬度或投影米）映射到格网像素索引。

def _interp_polygon(poly_pix, dx):
    """在规则步长下插值多边形顶点"""
    # 计算相邻顶点间距，在超过dx距离的顶点间线性插值新点
    dist = np.sqrt(np.sum(np.diff(poly_pix, axis=0)**2, axis=1))
    for i in range(len(dist)):
        if dist[i] > dx:
            n_interp = int(dist[i] / dx)
            # 在poly_pix[i]和poly_pix[i+1]之间线性插入n_interp个点
    ...

5.4.2 栅格化（rasterization）与冰原岛处理

多边形栅格化（rasterization）使用 skimage.draw.polygon 将多边形边界内的所有像素标记为冰川。对于含有内环（holes/nunataks）的复杂多边形，OGGM通过Shapely分别提取外环和内环——外环填充冰川掩码，内环在掩码中除外（设回0）。此外，函数使用 scipy.ndimage.label 识别连通组件，移除面积小于最大连通区5%的孤立小斑块（ice islands，即与主体冰川掩码分离的像素团）。

最终输出写入 gridded_data.nc：glacier_mask（二值掩码，glacier=1 else=0）和 glacier_ext（冰川外轮廓内的所有像素，包括nunatak）。

5.4.3 简化掩码方法（simple_glacier_masks）

simple_glacier_masks（约第1034-1163行）是替代方案，当不需要亚像素精度时使用 rasterio.features.rasterize（riomask模式）直接栅格化（rasterization）多边形，省去了坐标变换中的插值步骤，速度更快但几何精度略低。

5.5 网格参数选择（glacier_grid_params）

glacier_grid_params（约第226-302行）决定每个冰川在其地图投影中的网格分辨率（dx）。这是OGGM中最精妙的空间尺度自适应算法之一。

5.5.1 四种分辨率方法

方法	公式/逻辑	适用场景
`linear`	dx = d1 * sqrt(area) + d2	默认方法，小冰川细网格、大冰川粗网格
`square`	dx = sqrt(area) / 50	适用于对测高数据敏感的冰厚反演
`fixed`	dx = grid_dx_fixed（用户指定）	对区域同质网格进行批量处理
`by_bin`	按面积分箱，每箱一个固定dx	便于与其他模型进行区间比较

linear 方法的两个系数 d1 和 d2（默认14.0和0.0）由cfg.PARAMS设置，其背后的物理直觉是：冰川的宽度大致与其长度的平方根成正比，因此像素数应保持恒定。参数dx_scale 用于全局缩放所有derived的分辨率。

最终选定的dx、地图投影参数（由 salem 库基于冰川质心位置确定合适的UTM区或Albers等面积投影）、以及网格的nx和ny都写入glacier_grid.json。

5.6 DEM重投影（reproject_dem）

reproject_dem（约第351-453行）使用 rasterio.warp.reproject 将地理坐标系的DEM重投影到每条冰川的特定投影坐标系。重采样算法默认为双线性（bilinear），对于坡度等派生变量可用三次（cubic）。

def reproject_dem(gdir):
    # 确定目标投影（UTM / Albers）和网格参数
    proj_out = gdir.grid.center_grid  # 从 glacier_grid.json 读取

    with rasterio.open(src_path) as src:
        reproject(
            source=src.read(1),
            destination=dest_array,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=dst_transform,
            dst_crs=proj_out,
            resampling=Resampling.bilinear)

重投影后的高程存储在GlacierDirectory的本地GeoTIFF中，并作为后续所有二维计算的底图。这一步的精度直接决定了网格属性（坡度、坡向）的可靠性。

5.7 网格属性计算（gridded_attributes）

gridded_attributes（约第1407-1523行）计算多个二维场，这些场是后续中心线搜索和冰厚反演算法的必要输入：

属性	计算方法	用途
`slope`	`np.gradient(topo_smoothed, dx)` 的模	冰厚反演（驱动力）、中心线搜索
`aspect`	`np.arctan2(dy, dx)`	海拔带分析、太阳辐射
`slope_factor`	max(tan(slope), min_slope)	反演中防止近零坡度导致的数值爆炸
`dis_from_border`	`scipy.ndimage.distance_transform_edt` 取欧氏距离	中心线成本网格、按海拔分配厚度
`ice_divides`	分水岭算法在水文流向图上执行	支流划分、集水区面积计算

其中 slope_factor 的 min_slope 默认值为 tan(1.5°)，低于此坡度的像素被提升到此阈值，避免冰厚反演中除零或产生无穷大厚度。

5.8 测高曲线计算（compute_hypsometry_attributes）

compute_hypsometry_attributes（约第1165-1274行）以50m为固定高程箱计算每个海拔段内的面积及其占冰川总面积的千分比（permil）。输出为 hypsometry.csv，表头包括topo_bin、area、permil三列。

# 简化逻辑
bins = np.arange(min_elev // 50 * 50, max_elev + 50, 50)
for i in range(len(bins)-1):
    area_in_bin = np.sum((topo >= bins[i]) & (topo < bins[i+1])
                         & glacier_mask) * dx**2
    hyps[i, :] = [bins[i]+25, area_in_bin, area_in_bin / total_area * 1000]

测高曲线是物质平衡校准和面积-高度分布匹配的关键输入，在第8章的物质平衡模型中将反复引用。

5.9 GIS任务的完整工作流

gis_prepro_tasks 函数（约第1850行）将所有上述步骤编排为任务图。这允许OGGM以任意顺序构建依赖关系并并行化处理：

def gis_prepro_tasks(gdirs):
    task_list = [
        taskprocess_dem,                     # 1. DEM处理
        taskglacier_masks,                   # 2. 冰川掩码
        taskglacier_grid_params,             # 3. 网格参数
        taskglacier_mu_candidates,            # 4. 物质平衡参考期
        taskcompute_centerlines,             # 5. → 第6章
        ...
    ]

5.10 关键参数速查

参数名所属函数默认值描述

border process_dem 40 冰川周围缓冲区的像素数

grid_dx_method glacier_grid_params 'linear' 网格分辨率选择方法

d1 glacier_grid_params 14.0 linear方法的缩放系数

d2 glacier_grid_params 0.0 linear方法的偏移量

dx_scale glacier_grid_params 1.0 全局网格分辨率缩放

topo_interp process_dem 'linear' 空洞填充插值方法

smooth_window process_dem 3.5 高斯平滑窗口（像素）

clip_dem0 process_dem False 低于海平面裁剪

min_slope gridded_attributes tan(1.5°) 最小坡度因子阈值

5.11 总结

gis.py 是OGGM的基础模块——每一寸冰川模型的构建都始于这里的DEM处理。从空洞填充到高斯平滑，从多边形栅格化（rasterization）到网格分辨率自适应选择，这些看似基础的操作中蕴含着丰富的冰川学工程权衡：如何在精度和计算效率之间取得平衡，如何处理真实DEM中的各种瑕疵，如何让同一套代码适用于从0.1km²到10000km²的冰川。掌握GIS模块的细节，对理解后续所有模型输出的物理含义至关重要。

下一章将进入OGGM最具原创性的算法领域——冰川中心线提取和几何参数化。在那里，将看到本章生成的地形属性（坡度、坡向、边界距离）如何被转化为中心线优化算法的成本函数。

第4章：GlacierDirectory 数据中枢第6章：冰川中心线与几何参数化

OGGM 源码深度解析 © 2025 — 一本面向冰川研究者的开源模型技术指南