河南省30米DEM数据处理实战：GDAL拼接、裁剪与坡度计算

简介：河南省30米分辨率DEM数字高程模型，适合GIS、遥感或地理测绘领域的从业者与学生使用。数据基于ASTER GDEM V3拼接而成，采用GeoTiff格式与WGS84坐标系，覆盖河南全省范围，可直接用于地形渲染、坡度坡向提取、流域分析与地貌特征研究。压缩包共10个文件，核心为HeNan_DEM_30m_ASTGTMV003.tif栅格数据，同时包含河南省边界Shapefile全套文件（shp/shx/dbf等），以及tfw、xml等空间参考辅助文件，整体约176.18MB。该资源提供河南省行政边界和DEM栅格一体化配套，用户无需再自行拼接裁剪，可节省预处理时间。目前已有590人学习下载，适合需要快速获取区域地形数据进行授课、课题研究或项目演示的GIS用户。

1. 这份河南省DEM不是一张普通的tif，先分清栅格与矢量再动手

拿到河南DEM.rar解压后，你看到的不是一张图，而是一组配套文件：主文件HeNan_DEM_30m_ASTGTMV003.tif是30米分辨率的栅格高程，旁边的河南省.shp、.dbf、.prj才是行政区边界。如果只把 tif 拖进查看器，很容易误以为“这就是一张地形灰度图”，但真正做坡度、裁剪、面积统计时，栅格和矢量必须分开处理。这份数据来自 ASTER GDEM V3，按河南省范围拼接好，并附了边界 shp，适合省级尺度的地形分析、流域提取、灾害普查和制图底图。下面按文件拆解、拼接裁剪、地形分析和数据核查四个阶段展开，每一步给出可复现的命令和参数，避免你在 WGS84 经纬度坐标系上踩“算出来的坡度明显偏小”这类坑。

2. ASTER GDEM V3的文件清单：GeoTIFF、TFW、AUX.XML与shp边界

2.1 30米分辨率与ASTER GDEM V3的精度边界

30米分辨率表示每个像素对应地面约30m×30m，像素值是该单元内的地表高程。河南全省面积约16.7万平方公里，按30m栅格计算理论上有上亿个像素，所以这里的HeNan_DEM_30m_ASTGTMV003.tif体积不会是几MB的小图，而是一个需要分层金字塔才能流畅浏览的大栅格。ASTER GDEM V3 是 NASA 与 METI 联合发布的第三代全球高程模型，采样间隔1弧秒，在赤道约等于30米，纬度越高经度方向的实际距离越短。V3 相对 V2 主要在空洞填补、水体掩膜和异常像元上做了修正，但不同区域的验证精度差异很大，平原地区垂直中误差可能在10米内，山区会更大。所以把它用在省级规划、地形分类、流域提取上是合适的，用于工程勘察或建筑单体设计就不够。

河南地形本身也值得注意：西部伏牛山、北部太行山、南部桐柏山与大别山，东部是黄淮平原。这套30米数据能分辨出山脊、沟谷和河流阶地的宏观走向，但到了单条冲沟或小规模切坡的尺度，像素平均效应会让微地形被抹平。理解这一点，后面做坡度分级和地貌分割时就不会对细节抱有过高期待。

2.2 文件包里的每个文件是干什么的

解压后的文件列表看起来杂乱，但可以归成三组：栅格主文件、矢量边界文件、辅助元数据文件。栅格主文件就是HeNan_DEM_30m_ASTGTMV003.tif，它内部已经嵌入了坐标信息和仿射变换参数，这是绝大多数 GIS 软件读取位置的第一来源。河南省.shp、.shx、.dbf、.sbn、.sbx是一套完整的 Shapefile 边界数据，.sbn/.sbx是 ESRI 的空间索引，在只读边界时不参与，但被移动后会触发警告。.prj是文本坐标定义，用 WKT 格式记录边界使用的坐标系。

辅助文件里有三个容易被忽略：HeNan_DEM_30m_ASTGTMV003.tfw是外部世界文件，.tif.aux.xml是 GDAL 生成的元数据缓存，.tif.vat.dbf是 ArcGIS 的栅格属性表。它们的职责如下表：

这里需要说明，.tfw和 TIF 内嵌坐标信息一旦不一致，结果以 TIF 内部为准。很多人在 ArcGIS 里看到“栅格没有金字塔”的提示，实际是因为.aux.xml丢失，重新计算即可，不会损坏数据。.aux.xml里还会缓存STATISTICS_MINIMUM和STATISTICS_MAXIMUM，如果你用gdalinfo -stats发现统计值明显不对，先删掉.aux.xml再算一次。

2.3 用gdalinfo检查GeoTIFF与WGS84的实际情况

无论你用什么软件，第一步都应该用 GDAL 命令行读一下真实元数据，避免被文件命名误导。打开终端执行：

gdalinfo HeNan_DEM_30m_ASTGTMV003.tif

重点关注几行输出：Size is后面的宽高像素数；Origin是左上角经纬度；Pixel Size = (0.000277777777778, -0.000277777777778)代表大约1弧秒，即赤道附近30.9米；Coordinate System is: GEOGCS["WGS 84"...]表示目前是经纬度坐标，单位是度而不是米。如果看到NoData Value=*这一行，记住这个值，后续拼接和裁剪时统一使用。

注意一个容易被忽略的坑：在河南所在的北纬32—36度区域，经度方向1弧秒对应的地面距离约25—26米，纬度方向仍是约30.9米。也就是说，这份数据在WGS84经纬度下严格说并不是正方形的地面分辨率。相应地，直接用度为单位计算坡度、坡向时，必须引入尺度因子。这个细节我放到第4章展开。

3. 用GDAL把ASTER GDEM V3分幅数据拼成河南省30米DEM并裁剪

3.1 拼接前先统一坐标系：gdalbuildvrt与gdal_translate

ASTER GDEM V3 原始发布是按 1°×1° 分幅的，每幅一个GeoTIFF，文件名通常包含经纬度信息，比如ASTGTMV003_N33E112_dem.tif。你拿到的这份河南DEM如果来自“自己拼接”，大概率就是先把覆盖河南的相邻分幅找到，然后合并。合并的常见做法是先建 VRT 虚拟栅格，不复制数据，检查范围无误后再输出为正式的GeoTIFF：

# 构建虚拟拼接，只记录文件位置和范围，不复制像素 gdalbuildvrt -srcnodata -9999 -vrtnodata -9999 henan_tiles.vrt ASTGTMV003_N3*E1*_dem.tif # 转成正式GeoTIFF，统一数据类型和无值标记 gdal_translate -ot Float32 -a_nodata -9999 -co COMPRESS=DEFLATE -co TILED=YES henan_tiles.vrt henan_dem_merged.tif

这里-srcnodata -9999 -vrtnodata -9999是很多教程不会写的关键参数。原始分幅边缘的NoData值可能各不相同，有的写0，有的写-9999，在VRT阶段统一设置为-9999后，后续裁剪和统计才不会被边缘黑边干扰。-ot Float32保证高程精度，COMPRESS=DEFLATE是无损压缩，TILED=YES让大文件按块读取，分析速度比按行读取快很多。通配符ASTGTMV003_N3*E1*_dem.tif要按你实际下载的分幅文件路径修改，我这里是按河南覆盖范围示例。

如果你手头已经有拼好的HeNan_DEM_30m_ASTGTMV003.tif，可以跳过拼接，直接进入裁剪。但如果发现最终结果有接缝线或重叠重影，就要回头检查参与拼接的分幅是否都来自V3版本，V2和V3混用会出现明显的高程阶跃。

提示：拼接前先用gdalinfo抽查两个分幅的NoData值，确保一致，否则边缘会出现异常条纹。

3.2 用河南省shp进行裁剪：gdalwarp的cutline参数

拿到河南省边界 shp 后，最常见的错误是用gdal_translate -projwin按矩形裁剪，这样只切出外接矩形，河南省以外的区域仍然留在图里，后续面积统计会虚高。正确做法是用gdalwarp的-cutline参数做矢量边界裁剪：

gdalwarp -cutline 河南省.shp -crop_to_cutline -dstnodata -9999 \ -tr 0.000277777777778 0.000277777777778 -r bilinear \ henan_dem_merged.tif henan_dem_cut.tif

参数含义分别是：-cutline指定裁切矢量面；-crop_to_cutline把输出范围收缩到裁剪面的外接矩形，并让边界外像素变成无效值；-dstnodata -9999统一输出无效值；-tr手动指定输出分辨率，这里沿用原始约1弧秒值，避免重采样改变地面单元大小；-r bilinear是重采样方式，对高程数据来说，双线性插值比最邻近法平滑，但在地形断裂处可能轻微拉低山峰值。如果数据已经投影到平面坐标，-tr就换成30 30（米）。

另一个常见问题是 shp 坐标是 WGS84 而 DEM 已经是投影坐标。gdalwarp 会自动根据 shp 内带的.prj做重投影，但依赖网络查找栅格；为了可控，建议先在原数据上检查是否一致。若不一致，先用ogr2ogr把 shp 转成对应投影再裁剪。

3.3 裁剪后的金字塔与NoData设置

裁剪完成不代表可以直接拿去分析，还要做两件事：清理0值伪有效像素和建立金字塔。某些重采样算法会把边界外的NoData写成0，而0对高程数据可能是真实海拔，但河南省最低点不可能为0，所以出现0值大概率是NoData被覆盖了。用gdal_calc.py清洗：

# 把0值重写为-9999，避免误当真实高程 gdal_calc.py -A henan_dem_cut.tif --outfile=henan_dem_clean.tif \ --calc="where(A==0, -9999, A)" --NoDataValue -9999 # 建立金字塔，提升ArcGIS/QGIS缩放浏览速度 gdaladdo -r average henan_dem_clean.tif 2 4 8 16

gdal_calc.py的表达式where(A==0, -9999, A)把等于0的像素改成NoData；如果你确认数据中存在真实0海拔点，就把条件改成结合其他无效值，或者干脆不处理。gdaladdo -r average在栅格内建多层降采样金字塔，缩放渲染时性能提升明显。注意：gdaladdo是对原文件追加内部金字塔，不是生成新文件，所以运行后 tif 体积会略微增大。

4. 基于30米DEM的坡度、山体阴影与高程重分类实战

4.1 对WGS84经纬度DEM计算坡度时，必须显式指定水平比例因子

如果你把这份河南DEM直接丢给 ArcGIS 的 Slope 工具，大概率会得到一张坡度整体偏小的图，因为默认梯度计算假设水平和垂直单位一致。原始数据是 WGS84 经纬度，水平单位是度，垂直单位是米，所以必须提供比例因子。GDAL 的gdaldem slope提供了-s参数，表示水平单位与垂直单位的比例。在纬度约 34° 的河南，1 度约等于 111 公里，即 111000 米左右，所以：

# 计算坡度，-p输出百分比，-s指定水平/垂直比例 gdaldem slope henan_dem_clean.tif henan_slope.tif -p -s 111320 # 生成山体阴影，用于可视化底图 gdaldem hillshade henan_dem_clean.tif henan_hillshade.tif -z 2.0 -az 315 -alt 45

解释：-p让输出值变成百分比坡度，不熟悉百分比的读者注意，45° 对应的百分比是100%；-s 111320是把经纬度水平坐标换算成米的系数，如果不加，坡度会被低估约5—10倍。hillshade里的-z 2.0是垂直拉伸系数，河南省东部平原海拔变化小，拉伸到2.0能看出微起伏；-az 315表示光源来自西北方向，-alt 45是太阳高度角，这两个参数只影响渲染效果，不做分析时不必严谨。

如果你预先用gdalwarp -t_srs EPSG:32650（UTM 50N）把DEM转成米制投影，gdaldem slope就不需要-s了，这是更保险的做法，尤其在做坡度分级统计和流域提取时，建议先投影再计算。

4.2 DEM、DSM、DOM和DTM：这份数据到底适合做什么

热词里“dsm生成dem”和“dem裁剪”经常放在一起讨论，但首先得分清DEM和DSM。DEM（Digital Elevation Model）通常指地表裸地高程，也叫DTM；DSM（Digital Surface Model）则记录地表物体顶面，比如树冠、屋顶；DOM是正射影像，不记录高程。它们的差别决定了你能从这份30米数据里提取什么，不能提取什么。

ASTER GDEM V3 属于 DEM，由于光学生成方式，它实际上在一些地形突变区域会混入树冠或建筑顶部的信号，但整体还是以地表为主。所以在河南省范围内做山地平原分类、河流提取、洪水淹没模拟都是可用的；如果要计算建筑高度、电线塔净空这类精细场景，必须换 LiDAR 生成的 DSM 或更高精度 DEM。如果手头只有DSM想得到DEM，通常要做地面滤波，常见方法包括渐进形态学滤波和布料模拟滤波，但这不在ASTER GDEM的工作范围内。

4.3 把连续高程分割成平原、丘陵和山地：栅格重分类

“gis中dem如何分割”本质是把连续高程值离散成类别。常见做法是gdal_calc.py做多条件组合，这也比在 ArcGIS 里点选 Reclassify 更可复现：

# 分三类：<200m平原，200-500m丘陵，>=500m山地 gdal_calc.py -A henan_dem_clean.tif --outfile=henan_relief_class.tif \ --calc="(A<200)*1 + (A>=200)*(A<500)*2 + (A>=500)*3" \ --NoDataValue -9999

这个表达式的逻辑是：每个布尔条件在成立时返回1，不成立返回0，三类条件互斥，相加后得到1、2、3。阈值可以根据项目调整，比如做大别山灾害分区时把山地下限提到800米。分割结果出来后，再用gdalinfo -hist查看直方图：

gdalinfo -hist henan_relief_class.tif

输出里的Count字段是每个类别的像素数。用像素数乘以单个像素面积即可估算各类面积。需要注意，WGS84下单个像素在不同纬度面积不同，做省级面积统计前最好先投影到 Albers 等积投影，否则平原丘陵比例会有偏差。要给分类结果出一张彩色图，可以用颜色表配合gdaldem color-relief：

# 1=平原(绿)，2=丘陵(黄)，3=山地(棕) printf "1 180 220 120\n2 255 230 120\n3 180 120 80\n" > relief_color.txt gdaldem color-relief henan_relief_class.tif relief_color.txt henan_relief_class_color.tif

5. 数据核查与再加工：TFW、VAT和坐标系一致性检查技巧

5.1 别让TFW的旧值带偏你

.tfw是外部世界文件，但 GeoTIFF 内部同样存有仿射变换参数。多数软件（QGIS、ArcGIS）优先读内部坐标，TFW只在某些老工具和脚本里被引用。所以当你用gdalwarp重采样后，根目录里残留的旧.tfw会与新TIF不匹配。检查方法是打开TFW文件，第1行是X方向像素宽度，第4行是Y方向像素高度（通常为负），第5行和第6行是左上角坐标。然后执行：

# 查看TIF内嵌仿射变换和尺寸，与tfw对比 gdalinfo -json HeNan_DEM_30m_ASTGTMV003.tif | grep -E "geoTransform|size"

对比geoTransform[1]和geoTransform[5]与TFW前两行的数值是否一致。不一致时，用文件管理器删除旧TFW，再运行gdal_translate -co TFW=YES重新生成。同时注意，.tif.aux.xml在重采样后同样需要删除，否则会沿用旧的统计最小值和最大值，导致渲染拉伸异常。

5.2 删掉过期的VAT，再校验边界与栅格范围

.vat.dbf是 ArcGIS 为分类栅格生成的属性表。用 GDAL 做裁剪、重投影后，像素值和像素数量改变，旧的VAT文件已经失效，却仍然躺在目录里。ArcGIS 打开时检测到同名.vat.dbf会使用它，造成“属性表和栅格不匹配”的假报错。我的做法是：凡是经过 gdal_translate 或 gdalwarp 处理的输出目录，如果原目录存在.vat.dbf，直接删除：

# 删除过期栅格属性表 rm -f HeNan_DEM_30m_ASTGTMV003.tif.vat.dbf

最后校验边界是否匹配。查看 shp 范围：

ogrinfo -so 河南省.shp

查看 DEM 范围：

gdalinfo henan_dem_clean.tif | grep "Corner Coordinates"

比较两组四角坐标：如果 shp 是投影坐标而 DEM 是经纬度，则先执行ogr2ogr -t_srs EPSG:4326 河南_wgs84.shp 河南省.shp再对比。只有边界范围和坐标系都对齐，后续裁剪、坡度计算、面积统计才可信。这是交付给下游同事前最值得花的两分钟。

本文还有配套的精品资源，点击获取