北京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。文件后缀类型作用.tifGeoTIFF分类栅格本体携带内嵌 GeoTransform.tfwASCII 世界文件记录像元尺寸、旋转和左上角坐标.vat.dbfdBASE 表分类码到地类名称的属性表.aux.xmlXML缓存统计信息、调色板与坐标描述下面用命令解压并检查文件类型mkdir -p beijing_landform 7z x 北京市地形地貌最新30m精度.rar -obeijing_landform cd beijing_landform file landform_北京市.tif landform_北京市.tfw landform_北京市.tif.vat.dbfWindows 下用 7-Zip 或 WinRAR 直接解压Linux 下用 7z 即可。中文文件名在部分 Linux 环境里解压后可能显示为乱码可以用convmv -f GBK -t UTF-8 --notest *批量转码。解开后先用file确认 tif 是 GeoTIFFtfw 是 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_AreaProj4 里会出现projaea同时带lat_1、lat_2、lon_0等参数。这里有一个常见的坑这类分类栅格经常是调色板文件ColorInterp显示为Palette所以 gdalinfo 看到的Band 1 Block... TypeByte是正常现象不要误判成普通 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 把北京市海拔分成低海拔、中海拔、中高海拔、高海拔、极高海拔五档。严格的分级界值在各学科的用法不完全一致这类产品常见做法是等级海拔范围米典型分布低海拔小于 1000平原与山间盆地中海拔1000 - 2000低山过渡带中高海拔2000 - 3000中山山地高海拔3000 - 5000高山区域极高海拔大于 5000极高山北京的地势整体西北高、东南低西部和北部是西山、军都山所以市区和平原区基本落在低海拔档门头沟、延庆部分区域进入中海拔以上。使用这份分类图的关键是把分类码映射成可读名称而不是直接去读 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, encodingutf-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_countsTrue) 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 projaea lat_125 lat_247 lon_0105 datumWGS84 \ -r near \ -tr 30 30 \ -overwrite \ 海拔分类_北京市.tif 海拔分类_beijing_aea.tif-t_srs指定目标投影-r near指定重采样算法为最邻近分类栅格必须用 near不能使用 bilinear 或 cubic否则会在类别边界插值出不存在的新分类码-tr 30 30表示输出像元分辨率 30 米。如果数据包内有 .prj 文件投影参数要优先以它为准。gdalwarp 参数作用分类栅格建议-r near最邻近重采样必须使用-tr 30 30输出像元大小与源数据一致-overwrite覆盖已存在文件避免残留旧结果-dstalpha输出 Alpha 波段裁剪时推荐4.2 按区界或图幅裁剪北京区域在某些分析里只需中心城区六区或者按街道边界裁剪。使用矢量边界文件裁剪栅格时注意先把边界转到同样的投影避免坐标对齐问题ogr2ogr -t_srs projaea lat_125 lat_247 lon_0105 datumWGS84 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_width30.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_countsTrue) 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(fclass {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 脚本组成的批处理在更换其他城市数据时只改路径和边界文件即可复用。本文还有配套的精品资源点击获取