北京30m地形地貌栅格处理：GDAL解包、投影与面积统计

简介：北京市最新30m精度地形地貌数据包，依据海拔、起伏程度与成因形态，将北京市划分为低海拔至极高海拔、丘陵至极大起伏、平原山脉沟壑等地貌类型，并区分海积、湖积、冲积、洪积、风积、冰碛等成因。面向GIS专业学生、规划人员和地理爱好者，适合用于制图、空间分析与区域地貌研究。数据采用WGS84与Albers_Conic_Equal_Area坐标，精度达30m，已按省整理为TIF栅格格式。压缩包共30个文件，约1.41MB，核心为TIF栅格数据，配套TFW坐标文件、XML元数据、CPG编码与DBF属性表，另有使用说明RAR、概览PNG，结构完整便于直接加载。目前已有193人学习，可作为北京市地貌专题图、土地规划或教学案例的基础数据，直接导入ArcGIS、QGIS等软件查看与后续分析。

1. 北京 30m 地形地貌栅格到底能用来做什么

做场地选址、防洪评估或者风环境模拟的时候，DEM 只告诉你“海拔多高、坡度多陡”，但它不说这块地是冲积平原还是冰碛垄，是低海拔丘陵还是中起伏山地。拿到这份北京地形地貌 30m 精度栅格数据，等于把宏观地貌类型、按海拔划分的等级、按起伏度划分的等级一次性嵌进栅格属性里。数据以 TIF 组织，附带世界文件和 DBF 属性表，可以直接进 ArcGIS/QGIS，也可以用 GDAL 批处理。适合做国土空间规划分析、工程建设适宜性评价的从业者，以及所有需要在 Python 里处理分类栅格的读者。

2. 解包与栅格元数据核查

拿到压缩包后，第一件事不是拖进 GIS 软件，而是先看栅格的投影、像元尺寸、波段数和属性表是否配对。这套数据目录里有不少同名文件，例如 landform_北京市.tif、landform_北京市.tfw、landform_北京市.tif.aux.xml、landform_北京市.tif.vat.dbf，它们分工不同，少一个都可能影响符号化和属性读取。

2.1 解压后先分清六个关键文件

压缩包解开后，按功能可以分成三层。第一层是栅格本体，如 landform_北京市.tif、海拔分类_北京市.tif、起伏程度分类_北京市.tif、陆地地貌类型_北京市.tif、中国地貌类型_北京市.tif。第二层是空间定位文件，包括 .tfw 和 .aux.xml，前者是 ESRI 风格的世界文件，记录左上角坐标和像元尺寸；后者是 GDAL 自动生成的辅助元数据，记录色彩解释和统计信息。第三层是属性表文件，.vat.dbf 和 .vat.cpg，前者是带分类代码和计数值的 dBASE 表，后者是表编码，通常是 UTF-8。

下面用命令解压并检查文件类型：

mkdir -p beijing_landform 7z x 北京市地形地貌最新30m精度.rar -obeijing_landform cd beijing_landform file landform_北京市.tif landform_北京市.tfw landform_北京市.tif.vat.dbf

Windows 下用 7-Zip 或 WinRAR 直接解压，Linux 下用 7z 即可。中文文件名在部分 Linux 环境里解压后可能显示为乱码，可以用convmv -f GBK -t UTF-8 --notest *批量转码。解开后先用file确认 tif 是 GeoTIFF，tfw 是 ASCII 文本，vat.dbf 是 dBASE 表，三者类型正常再继续，避免后续 GDAL 读取时报格式错误。

2.2 用 gdalinfo 核对投影、像元尺寸与分类栅格特征

gdalinfo landform_北京市.tif gdalinfo 海拔分类_北京市.tif

输出里重点看Size is，确认 x 方向和 y 方向的像元数；接着看Pixel Size，期望接近 30 米；然后看Coordinate System is，这套数据同时提供 WGS84 和 Albers_Conic_Equal_Area 两套坐标描述，gdalinfo 会列出具体的投影参数和基准面信息。

gdalinfo -proj4 landform_北京市.tif

如果是 Albers_Conic_Equal_Area，Proj4 里会出现+proj=aea，同时带+lat_1、+lat_2、+lon_0等参数。这里有一个常见的坑：这类分类栅格经常是调色板文件，ColorInterp显示为Palette，所以 gdalinfo 看到的Band 1 Block=... Type=Byte是正常现象，不要误判成普通 RGB 影像。类型为 Byte 意味着分类码在 0 到 255 之间，大概率是地貌分类代码，而不是高程真值。

2.3 用 Python 读取 vat.dbf 验证分类码完整性

GDAL 把 .vat.dbf 当作栅格属性表管理，严格说是 ESRI 的 Raster Attribute Table。不需要手动打开 DBF，直接读栅格，再用分层统计值交叉验证。

import numpy as np try: from osgeo import gdal gdal.UseExceptions() except ImportError: raise RuntimeError("请安装 GDAL 的 Python 绑定") src = gdal.Open("landform_北京市.tif", gdal.GA_ReadOnly) band = src.GetRasterBand(1) data = band.ReadAsArray() print("行列数:", data.shape) print("最小分类码:", int(data.min()), "最大分类码:", int(data.max()))

读过以后，最小值和最大值应落在 1 到 30 左右的区间，如果出现 255 或者 0，说明存在 NoData 或背景值，后面统计面积前要单独处理。.vat.dbf里的Value和Count字段正好对应分类码和该分类的像元数量，用这两个字段可以在不扫描全图的情况下预先了解各地类的占比。由于这张图是 Byte 分类，数据量不大，直接 NumPy 扫描也很快，但生产环境里第一遍仍然建议先读 dbf，速度快且能发现属性表与栅格不同步的问题。

到这里，栅格本体、定位文件、属性表都验过了。下一步要弄清楚分类码到底代表什么，这是做后续分析的前提。

3. 分类体系与属性表连接逻辑

这套数据里的“地形地貌”不是单张 DEM 派生的坡度切片，而是把地貌学里的海拔分级、起伏分级与成因类型拆成四个主题栅格。它们之间通过分类码关联，理解编码体系后，你才能把 landform_北京市.tif 里的综合代码翻译成“低海拔冲积平原”这种可读语义。

3.1 海拔分类：从低海拔到极高海拔

海拔分类_北京市.tif 把北京市海拔分成低海拔、中海拔、中高海拔、高海拔、极高海拔五档。严格的分级界值在各学科的用法不完全一致，这类产品常见做法是：

北京的地势整体西北高、东南低，西部和北部是西山、军都山，所以市区和平原区基本落在低海拔档，门头沟、延庆部分区域进入中海拔以上。使用这份分类图的关键是把分类码映射成可读名称，而不是直接去读 TIF 的像元值。

3.2 起伏程度分类与陆地地貌类型的组合关系

起伏程度分类_北京市.tif 是另一套独立分级：丘陵、小起伏、中起伏、大起伏、极大起伏。陆地地貌类型_北京市.tif 则更接近地貌成因和形态的细分，包括山地、丘陵、平原、台地等，成因属性里还会出现冲积、洪积、湖积、海积、风积、冰碛、剥蚀侵蚀等字段。实际做工程评估时，往往是“海拔 + 起伏 + 成因”三表叠加得到最终地貌单元，例如“低海拔冲积平原”和“中海拔侵蚀丘陵”的工程意义完全不同，前者适合建设用地平整，后者需要边坡治理。

这三张图的像元尺寸与范围一致，可以直接做逐像元叠加。叠加前用一个 Python 字典维护分类码到中文语义的映射，比反复查 dbf 更直观。

import csv class_map = { 1: ("低海拔", "平原", "冲积"), 2: ("低海拔", "丘陵", "侵蚀"), 3: ("中海拔", "小起伏", "洪积"), # 实际编码以 vat.dbf 为准，这里只是演示结构 } rows = [] for v, names in class_map.items(): rows.append([v] + list(names)) with open("beijing_landform_class.csv", "w", newline="", encoding="utf-8") as f: writer = csv.writer(f) writer.writerow(["code", "elevation", "relief", "genesis"]) writer.writerows(rows) print(rows)

把 class_map 导出成 CSV 后，可以直接在 ArcGIS 里对栅格做 Lookup 或者 Join，也可以用 QGIS 的“栅格唯一值”功能挂接 .vat.dbf。分类码对应的名称在 dbf 的 Value 字段和别名属性里，先打印再写映射，不要凭经验猜。

3.3 用 GDAL 统计各分类的像元数与占比

属性表里的 Count 字段可以直接预统计，但为了和后续的投影转换衔接，这里用 GDAL 自带工具最省事：

gdalinfo -hist landform_北京市.tif

-hist会输出 0-256 共 256 个桶的直方图，这是一个比较粗糙的分布预览。精确统计用 Python 更快：

import numpy as np from osgeo import gdal src = gdal.Open("海拔分类_北京市.tif") data = src.GetRasterBand(1).ReadAsArray().astype(np.uint8) valid = (data != 0) & (data != 255) unique, counts = np.unique(data[valid], return_counts=True) for cls, cnt in zip(unique, counts): print(f"分类码 {cls}: {cnt} 个像元，占比 {cnt / valid.sum():.2%}")

这段代码先把 0 和 255 当无效值剔除，然后用np.unique统计有效分类码的像元数。占比计算时除以valid.sum()，只统计有效区域，避免背景值把比例拉低。如果发现某个分类码数量明显异常，回到 .vat.dbf 查一下它对应的名称，往往能发现是 NoData 值混进了分类。

到这里，编码体系和统计路径都已经跑通。下一步需要解决投影与面积统计的问题，因为 WGS84 经纬度坐标下的像元面积随纬度变化，不能在经纬度投影下直接算平方公里。

4. Albers 投影下的重投影、裁剪与面积统计

这份资料的坐标信息有两套：WGS84 经纬度和 Albers_Conic_Equal_Area。前者用于在互联网地图和 GPS 数据里对齐位置，后者用于面积量算和省级制图。两套坐标系统各有分工，以下处理中会把 Albers 当作分析基准，避免按经纬度像元算面积引起的错误。

4.1 为什么面积统计必须用等积投影

如果栅格停留在 WGS84，经纬度像元在地面的实际宽度是“赤道约 111 公里，北纬 40 度约 85 公里”，也就是说，一个 0.0003 度的像元在东西方向和南北方向的地面距离不一样，而且随着纬度变化。北京在北纬 39°26′ 到 41°03′，纬度跨度约 1.5 度，直接用 WGS84 像元数乘以固定常数会带来明显面积误差。Albers_Conic_Equal_Area 是等积投影，投影后像元面积和地面面积成常数比例，因而统计各分类面积才靠得住。

常见的 Albers 参数是中央经线 105°E、双标准纬线 25°N 和 47°N，这是中国省级和全国制图常用的配置。具体数值先读取 tif 的投影描述再作参考。下面是重投影命令：

gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +datum=WGS84" \ -r near \ -tr 30 30 \ -overwrite \ 海拔分类_北京市.tif 海拔分类_beijing_aea.tif

-t_srs指定目标投影；-r near指定重采样算法为最邻近，分类栅格必须用 near，不能使用 bilinear 或 cubic，否则会在类别边界插值出不存在的新分类码；-tr 30 30表示输出像元分辨率 30 米。如果数据包内有 .prj 文件，投影参数要优先以它为准。

4.2 按区界或图幅裁剪

北京区域在某些分析里只需中心城区六区，或者按街道边界裁剪。使用矢量边界文件裁剪栅格时，注意先把边界转到同样的投影，避免坐标对齐问题：

ogr2ogr -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +datum=WGS84" bj_aea.shp bj_bound.shp gdalwarp -cutline bj_aea.shp \ -crop_to_cutline \ -dstalpha \ -tr 30 30 \ 海拔分类_beijing_aea.tif 海拔分类_downtown.tif

-cutline指定裁剪矢量，-crop_to_cutline表示输出范围严格贴合边界，-dstalpha在输出文件里增加一个 Alpha 波段，把边界以外的像元标成透明。对分类栅格保留 alpha 通道有利于后续渲染，但面积统计时要通过掩膜过滤掉该区域。

4.3 基于像元数和 Albers 像元面积统计地类面积

等积投影下每个 30m×30m 像元面积是 900 平方米，等于 0.0009 平方千米。统计公式可以简化为“像元数 × 0.0009 平方千米”。完整代码如下：

import numpy as np from osgeo import gdal def class_area(tif_path, pixel_width=30.0): src = gdal.Open(tif_path) band = src.GetRasterBand(1) data = band.ReadAsArray() gt = src.GetGeoTransform() x_res, y_res = abs(gt[1]), abs(gt[5]) cell_area_km2 = (x_res * y_res) / 1_000_000.0 if data.dtype == np.uint8: valid = (data != 0) & (data != 255) else: valid = data > 0 vals, counts = np.unique(data[valid], return_counts=True) return vals, counts * cell_area_km2, cell_area_km2 vals, area_km2, per_cell = class_area("海拔分类_beijing_aea.tif") for v, a in zip(vals, area_km2): print(f"class {v}: {a:.2f} km² (单像元 {per_cell*1e6:.0f} m²)")

函数里首先从 GeoTransform 里动态读取 x 方向和 y 方向的分辨率，不硬编码 30 米，防止某些重投影操作后分辨率变化。然后对 Byte 分类栅格做 0/255 排除，用np.unique统计每个类别的有效像元个数。面积单位先转成平方千米，输出时保留两位小数。如果某个类别的面积和全市面积数量级相差很远，优先检查有没有把经纬度栅格误算进来。

有两点值得提醒：第一，统计时以重投影后的 GeoTIFF 为准，不要让原始 WGS84 栅格参与计算面积；第二，Count 字段来自原属性表，如果做过裁剪或重投影，必须重新统计，原 Count 已失去意义。

5. 生产环境验证：精度核对、异常值分析与叠加制图技巧

数据落到项目里之前，最好用 20 分钟做一次精度核对，否则分类图和 DEM 对比时出现系统性位移，前期的缓冲分析都会被带偏。这一章只讲三个最有效的验证手段。

5.1 用 .tfw 核对像元尺寸和左上角坐标

.tfw 是一个六行文本文件，直接查看内容即可验证空间分辨率：

cat landform_北京市.tfw

前两行是 x 方向和 y 方向的像元尺寸，第三、四行是旋转参数，通常为 0。第五、第六行是左上角的 X、Y 坐标。如果数据是 Albers 投影的 .tfw，这里的数值会落在北京区域的坐标量级。重点检查第一行的绝对值是否接近 30，第五行和第六行是否和 gdalinfo 里Upper Left一致。如果 gdalinfo 读出的坐标和 tfw 不一致，说明 tif 内嵌 GeoTransform 与外部世界文件冲突，这种情况通常以 tfw 为准。

5.2 异常值、NoData 与调色板问题

在地貌分类栅格里，0 和 255 往往是背景或无效值。处理方式是在分类统计中把它们排除。其次，注意 aux.xml 里面可能写死了旧的统计信息，如果裁剪或重投影后继续使用原目录的 aux.xml，某些 GIS 软件会直接读取旧的统计结果，导致唯一值列表刷新不出来。安全的做法是每次输出新文件后删除对应目录下的 .aux.xml，让软件重新计算统计值。调色板 TIF 在 QGIS 里需要设置“样式→渲染类型→单波段伪彩色”，在 ArcGIS 里则要在“符号系统→唯一值”下手动指定，否则可能被自动拉伸成连续色带，掩盖分类类型。

5.3 与高分辨率影像叠加检查边界

最直接的精度验证是把 landform_北京市.tif 和天地图影像或高分影像叠加，检查山脊线是否与影像上的地形转折一致。通常设置分类栅格不透明度 60%，把影像作为底图，山地和平原的边界应当与影像上的植被、阴影和纹理变化吻合。如果偏移超过一个像元，大概率是原始投影错误或 tfw 被误替换，这时候换用 WGS84 版本重新对齐。验证通过后再输出 PNG 制图，图例按分类顺序排列，不要按字母排序。全部操作可写成 GDAL 与 Python 脚本组成的批处理，在更换其他城市数据时只改路径和边界文件即可复用。

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