QGIS中如何从DEM提取点位高程值？

Question

在使用QGIS进行地形分析时，常需从DEM（数字高程模型）中提取特定点位的高程值。常见问题如下：当加载了DEM数据和包含兴趣点的矢量图层后，如何准确地将DEM中的高程信息提取并赋值到点要素的属性表中？用户在操作过程中常遇到插值不准确、坐标系不匹配或结果字段为空等问题。特别是在点位于DEM边缘或存在投影差异时，提取结果易出错。此外，对“采样”工具（如“栅格采样”或“Extract Raster Values”）的参数设置不当，也会导致高程值无法正确写入点图层。如何确保空间参考一致，并正确使用QGIS内置工具高效、精确地完成高程提取，是实际应用中的关键技术难点。

Answer 1

一、基础概念与数据准备

在使用QGIS进行地形分析时，从DEM（数字高程模型）中提取点位高程值是常见操作。其核心原理是基于点的地理坐标，在栅格数据中查找对应像元值或通过插值计算获取高程。

• DEM：通常是浮点型栅格图层，代表地表高程信息。
• 点矢量图层：包含兴趣点（如观测站、采样点等）的几何位置。
• 目标：将每个点所对应的DEM高程值写入其属性表中的新字段。

为确保后续处理顺利，应首先检查以下三项：

1. DEM和点图层是否均已正确加载；
2. 两者空间参考系统（SRS）是否一致；
3. 点要素是否位于DEM覆盖范围内。

二、空间参考一致性验证与处理

坐标系不匹配是导致高程提取失败的主要原因之一。即使视觉上点与DEM重叠，若投影不同，实际坐标存在偏移。

检查项	操作方法	建议工具
查看图层CRS	右键图层 → 属性 → 信息 → 坐标参考系统	QGIS 图层属性面板
统一投影	使用“重新投影图层”工具转换至相同坐标系	Processing Toolbox → Reproject Layer
动态重投影启用状态	项目设置 → CRS → 启用“在项目范围内启用'on-the-fly’变换”	Project Settings

三、高程提取的核心工具对比

QGIS提供多种方式实现栅格值采样，以下是常用工具及其适用场景：

• SAGA GIS Raster Values at Points：支持双线性插值，适合高精度需求。
• GDAL Sample Raster Values：基于最近邻法，速度快但精度较低。
• GRASS v.what.rast：适用于复杂拓扑环境。
• 内置“采样栅格”工具（Sample Raster Values）：最常用，集成于主菜单。

推荐优先使用“采样栅格值”工具（Vector → Analysis Tools → Sample Raster Values），因其用户友好且可自动创建输出字段。

四、操作流程详解

以“采样栅格值”为例，执行步骤如下：

1. 打开菜单栏：Vector → Analysis Tools → Sample Raster Values；
2. 选择输入点图层（Input Point Vector Layer）；
3. 选择待采样的DEM图层（Raster Layer to Sample）；
4. 设置输出字段名称（如 "elevation"）；
5. 确认所有参数后点击“运行”；
6. 完成处理后，打开点图层属性表，验证是否新增字段并填充数值。

注意：若结果字段为空，需排查点是否超出DEM范围或存在NoData区域。

五、边缘点与NoData处理策略

当点位于DEM边缘或处于NoData像元上方时，提取结果可能为空或异常。

# 示例：使用Python控制台检查NoData值
layer = QgsProject.instance().mapLayersByName('dem_tif')[0]
provider = layer.dataProvider()
print("NoData Value:", provider.sourceNoDataValue(1))
  

解决方案包括：

• 扩展DEM边界（使用Raster → Miscellaneous → Warp补边）；
• 对点图层进行缓冲区筛选，排除边缘可疑点；
• 在提取前使用“Fill NoData Cells”工具填补空洞。

六、高级技巧：批处理与自动化脚本

对于大规模点集或多DEM场景，手动操作效率低下。可通过Processing Modeler构建自动化流程。

graph TD A[加载DEM与点图层] --> B{坐标系一致?} B -- 是 --> C[运行采样工具] B -- 否 --> D[重新投影点或DEM] D --> C C --> E[检查输出字段是否完整] E --> F{是否存在空值?} F -- 是 --> G[执行插值填补或标记异常] F -- 否 --> H[导出结果Shapefile]

此外，可结合PyQGIS编写脚本批量处理多个区域：

import processing
for point_layer in point_layers:
    result = processing.run("qgis:samplerastervalues", {
        'INPUT': point_layer,
        'RASTERCOPY': dem_layer,
        'COLUMN_PREFIX': 'elev_',
        'OUTPUT': f'memory:{point_layer.name()}_with_elev'
    })
    QgsProject.instance().addMapLayer(result['OUTPUT'])