中国地貌栅格数据处理:从解压到面积统计与叠加分析
2026/9/20 23:30:51 网站建设 项目流程

简介:中国地貌栅格数据是一套覆盖全国的地貌分区栅格产品,按海拔与起伏度划分出平原、台地、丘陵及高中低山等大类,适合GIS、地理学界研究者及高校学生用于地貌分类制图、空间分析与教学示例。压缩包共7个文件,其中主要数据为tif栅格、ovr金字塔与tfw坐标参考文件,配合dbf属性表及docx代码表说明,可方便在ArcGIS等平台直接加载与查询。数据包仅3.92MB,便于下载分享,当前已有1231人学习参考。借助这套数据,读者可快速获取全国尺度的地貌类型空间分布,能够支撑区域对比分析、专题图制作与地貌相关课程案例设计,同时可延伸学习栅格数据组织方式与编码释义。

1. 中国地貌栅格数据.rar:打开这个文件前先弄清的几件事

如果你是从某个科研共享站或同事手里拿到名为「中国地貌栅格数据.rar」的压缩包,第一反应是解压后丢进 ArcGIS 或 QGIS 里出图,那你大概率会在前五分钟内遇到两个问题:一是图层显示出来是全黑或全灰,二是想按行政区统计面积时发现坐标系里的单位并不是米。这份数据本质上是把中国地貌的形态分类结果按网格像元存储的栅格数据集,每个像元的值并不代表高程或坡度,而是一个地貌类型编码,比如 1 表示中山、2 表示低山、5 表示平原。理解这一点之前,任何操作都会偏离方向。

这篇文章不讨论数据从哪来、谁发布的,只从一线处理角度讲清楚:地貌栅格数据内部的结构是什么、如何用 GDAL 把这些像元值变成可分析的信息、以及怎样跟全国城市形态栅格数据集或地下水位等辅助数据做叠加。适合刚接触栅格数据、但已经会点 Python 或命令行的人;如果你手里只有 ArcGIS 的工具栏,某些思路也能平移过去。

2. 地貌栅格数据的组织方式:格式、坐标系与地貌分类体系

2.1 栅格数据的常见封装格式,解压后先看文件尾

压缩包里的内容往往不是单一文件,而是一个文件夹,里面可能是 GeoTIFF、ERDAS IMAGINE 的 .img,或者 ESRI GRID 的 .adf 一组文件。三种格式用 GDAL 都能读,但处理方式略有不同。

格式文件特征适合场景常见问题
GeoTIFF.tif,单个文件,内嵌地理参考分发、Python 处理重压缩后可能丢坐标
IMG.img,可能带同名字典 .rrd.NET 或老项目元数据读取慢
ESRI GRID.adf 多个文件存储于同名文件夹ArcGIS 原生直接删文件夹会破坏数据

我一般解压后第一件事不是打开,而是用gdalinfo看一眼文件路径下的实际内容。.rar里的文件名可能含中文,尤其在地貌数据这种国内发布的数据集中,路径中带「中国地貌」或拼音是常见情况。GDAL 对中文路径支持尚可,但建议先复制到纯英文路径下,避免后面 Python 脚本或命令行工具出现file not found这类诡异问题。

2.2 坐标系与分辨率:先对齐再分析

地貌栅格数据常用的坐标系可能是 Albers 等积投影,也可能是 WGS84 经纬度。为什么强调这一点?因为地貌分类统计要计算面积,如果数据是 WGS84 但像元大小是 0.0083 度,直接按像元数乘 1000 米只会得到错误结果。等积投影适合面积计算,而经纬度适合做简单浏览。

分辨率的常见值有 90 米、500 米、1 公里,这决定了你统计结果的最小空间单元。用 GDAL 读取这些元数据只需一条命令:

gdalinfo -stats 中国地貌.tif

输出中重点关注这几行:

Size is 5360, 3720 Coordinate System is: PROJCRS["Albers_Conic_Equal_Area", ...] Origin = (72.0000000000000000,54.0000000000000000) Pixel Size = (0.0083333333000000,-0.0083333333000000)

这里Pixel Size是度,说明数据是地理坐标系,虽然名字可能写的是 Albers,但实际是伪等积或未投影。-stats会强制计算各像元值的分布,同时也暴露了另一个坑:地貌分类编码通常是整数,但文件里如果混入 NoData 会导致统计值都变成负数,比如 -9999。你要记住文件里的 NoData 数值,后面处理时单独处理。

2.3 地貌分类体系与像元值含义

中国地貌分类体系常见的是把形态与成因结合,编码 1-7 或 1-25 都有。你没有说明文档时,可以用gdalinfo -hist输出直方图,大致看出有哪些整数值,再反查分类标准。常见做法是读取唯一值列表:

gdalinfo -hist 中国地貌.tif | tail -20

如果是 7 类分类,通常 1 为平原、2 为台地、3 为丘陵、4 为低山、5 为中山、6 为高山、7 为极高海拔。但不同数据源的编码可能反过来,需要谨慎。不要用颜色区分,因为渲染默认是灰度,容易被两个相近编码骗过去。

3. 用 GDAL / Python 解包并预处理中国地貌栅格数据

3.1 从 .rar 中解压并校验文件完整性

解压 .rar 在 Windows 上用 WinRAR 最省事,但在 Linux 服务器上,我通常用unrar或 Python 的rarfile库。rarfile依赖系统里有unrar命令,所以优先直接用系统命令:

sudo apt install unrar unrar x 中国地貌栅格数据.rar

这里x保留路径结构,e会把所有文件平铺到当前目录,我一般用x防止多个同名文件覆盖。解压后建议用ls -la查看文件大小,如果 .tif 只有几十 MB,而压缩包本身几百 MB,说明可能只是图层索引或者是分块压缩。有些压缩包内还带说明 PDF 或分类表,一并解出来看。

3.2 读取并检查数据有效性:用 Rasterio 打开像元值

Python 端我更喜欢用rasterio,它比 GDAL 命令更符合现代 Python 风格。先读取元数据和数据切片:

import rasterio import numpy as np with rasterio.open('中国地貌.tif') as src: print(src.crs) # 坐标系 print(src.transform) # 仿射变换参数 data = src.read(1) # 读取第一波段 print(np.unique(data)) # 唯一值列表 print(src.nodata) # NoData 值

这段代码的优势是直接看到数据里到底有哪些值。实际操作中我遇到过一个问题:np.unique(data)输出里出现nan,因为数据里用了浮点型但无效值没有标成 NoData。这时候你要决定是把它当成 NoData 还是当成真实地貌类型。如果压缩包里的说明文件没给,我一般做法是先把 NaN 替换成 -1,后面统计时剔除。

data = np.where(np.isnan(data), -1, data) np.save('tree_cover.npy', data) # 保存中间结果

3.3 重投影与裁剪:与全国城市形态栅格数据集对齐

地貌数据要和全国城市形态栅格数据集或地下水位栅格数据 shp 叠加时,必须先统一坐标系和分辨率。最常见的错误是直接加载两个图层看肉眼重合就以为对齐了,但像元边界可能差几十米,统计时产生误差。所以必须重投影。

rasterioreproject方法,或者直接用gdal.Warp。我更倾向于gdal.Warp因为它参数直观,适合命令行:

gdalwarp -t_srs "EPSG:3857" -r near -tr 1000 1000 -overwrite 中国地貌.tif 中国地貌_webmercator.tif

参数含义:

  • -t_srs目标坐标系,这里用了 Web Mercator,实际科研场景建议用目标数据的坐标系,比如城市形态数据是 UTM 50N 或 Albers,就以它为基准。
  • -r near重采样方法。地貌分类是离散类别,不能使用双线性或三次卷积,会插出 2.5 这类不存在的编码,必须用near
  • -tr 1000 1000目标分辨率 1000 米。如果两个数据分辨率不同,取两者中较大的那个,否则面积统计会有大量混合像元。

裁剪到行政区范围时,优先用省级行政区划 shp 做掩模:

gdalwarp -cutline 省级边界.shp -crop_to_cutline -dstnodata -1 -of GTiff 中国地貌_webmercator.tif 中国地貌_华东.tif

这里-cutline后的 shp 必须是有效多边形,且坐标系不需要预先转换,GDAL 会自动处理。-crop_to_cutline会把输出范围精确到 shp 的边界,而不是 shp 外接矩形,能节省磁盘空间。

4. 地貌栅格数据的统计分析与可视化

4.1 按地貌类型统计面积与占比

重投影后,如果坐标系是等积投影,像元面积就是x_resolution * y_resolution的固定值。用numpy统计各类别像元数:

import numpy as np import rasterio with rasterio.open('中国地貌_华东.tif') as src: data = src.read(1) pixel_area = src.transform[0] * -src.transform[4] # 像元宽 * 高,单位 map 单位 valid_mask = (data > 0) # 假设 0 不是地貌类型 labels, counts = np.unique(data[valid_mask], return_counts=True) categories = {1: '平原', 2: '台地', 3: '丘陵', 4: '低山', 5: '中山', 6: '高山', 7: '极高海拔'} for lab, cnt in zip(labels, counts): area_km2 = cnt * pixel_area / 1e6 print(f'{categories.get(lab, f"未知{lab}")}: {area_km2:.2f} km²')

这段代码里src.transform[0]是像元宽度,src.transform[4]通常是负的像元高度,乘出来就是单个像元的面积。如果坐标系不是等积投影,比如前面的 WGS84,就不能这么算。你需要用rasteriotransform.xy逐个像元计算经纬度实际距离,或者改用geopandas做栅格转面来计算。常见的替代方案是让gdalwarp先把数据投影到 Albers,再做统计。

4.2 制作专题图时的配色与标注坑

地貌栅格数据可视化最头疼的是类别多,且部分类别面积小,默认色带会把平原画成深绿色、山地画成浅黄色,但小图例很容易混淆。我一般使用matplotlibrasterio.plot,并显式指定离散色条:

import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap cmap = ListedColormap(['#d4e157', '#8bc34a', '#ffeb3b', '#ff9800', '#795548', '#9e9e9e', '#eceff1']) with rasterio.open('中国地貌_华东.tif') as src: data = src.read(1, masked=True) fig, ax = plt.subplots(1, 1, figsize=(8, 6)) im = ax.imshow(data, cmap=cmap, vmin=1, vmax=7) cbar = plt.colorbar(im, ax=ax, orientation='horizontal') cbar.set_ticks([1, 2, 3, 4, 5, 6, 7]) cbar.set_ticklabels(['平原', '台地', '丘陵', '低山', '中山', '高山', '极高海拔']) ax.set_title('华东地区地貌类型分布') plt.tight_layout() plt.savefig('华东地貌.png', dpi=300)

这里masked=True会自动把 NoData 变成透明,避免图上出现黑色色块。vmin/vmax限制色条范围为 1 到 7,否则imshow会根据数据实际最小值最大值拉伸。另一个容易踩的坑是:如果数据里存在 NoData 值 -9999,ListedColormap会把 -9999 也映射到第一个颜色,导致图上大范围出现绿色。所以一定要用 mask 而不是只设置 vmin。

4.3 直方图与数据质量检查

分析之前做一次直方图检查,能发现异常极值。rasterio可以快速生成统计:

python -c "import rasterio; src=rasterio.open('中国地貌_华东.tif'); print(src.statistics(1))"

如果看到最大值是 32767 或类似数值,而不是 7,说明数据里有异常值。这类异常值一般来自边界填充或分块拼接错误。我在处理全国城市形态栅格数据集时碰到过类似问题,拼接后边缘出现 0 值和 255 值混合,导致按 0-255 为有效范围统计时全部错误。解决办法是先做一次条件替换,把超出合法范围的值统一换成 NoData。

5. 进阶:把地貌栅格数据与城市形态栅格数据集叠加分析

5.1 提取不同地貌下城市形态指标

如果你手里还有全国城市形态栅格数据集,其中通常包含建筑密度、容积率、形态分维等连续值指标,可以按地貌类型分区做 zonal statistics。核心思路是:先把地貌数据重采样到与城市形态相同的网格,然后分组统计城市形态值。

gdalwarp -tr 100 100 -r near -t_srs EPSG:3857 中国地貌.tif 地貌_100m.tif

再用 Python 做分组聚合:

import rasterio import numpy as np import pandas as pd with rasterio.open('地貌_100m.tif') as geo_src, rasterio.open('城市形态_100m.tif') as city_src: geo = geo_src.read(1) city = city_src.read(1) valid = (geo > 0) & (city > 0) & (~np.isnan(city)) df = pd.DataFrame({'地貌': geo[valid], '建筑密度': city[valid]}) result = df.groupby('地貌').agg(平均密度=('建筑密度', 'mean'), 面积占比=('建筑密度', lambda x: (x > 0.5).mean())) print(result)

参数说明:-tr 100 100将地貌重采样到 100 米分辨率,如果城市形态数据原本就是 100 米,可以去掉这一步。geo > 0过滤掉无效地貌编码,city > 0过滤掉无效城市形态值。城市形态里的 0 往往代表非建设区,但有时代表数据缺失,需要结合说明判断。这种分组统计不适合用 ArcGIS 的 Zonal Statistics Table,因为地貌类型是离散值,直接用图斑提取更快。

5.2 验证地貌数据与地下水位栅格数据的关系

如果你后续想分析地貌类型与地下水位的关系,常见的栅格数据是以 shp 点文件或插值面形式出现的。比如我国地下水位栅格数据 shp 是矢量,你不能直接与地貌栅格逐像元计算。正确做法是先把地下水位面栅格化,或者用点数据在地貌栅格上采样。但要注意:地下水位是连续变量,地貌是分类变量,做相关性分析前必须检验空间自相关性。命令级别可以这样采样:

gdallocationinfo -geoloc -valonly 中国地貌.tif 东经 北纬

gdallocationinfo-geoloc表示输入的是地理坐标,输出对应像元值。对于大量点,建议用 Python 的rasterio.sample

from rasterio.sample import sample coords = [(120.1, 31.2), (121.3, 30.8)] # 地下水位监测点的坐标 with rasterio.open('中国地貌.tif') as src: values = [v[0] for v in sample(src, coords)] print(values)

注意sample中的坐标顺序是 (x, y) 即 (经度, 纬度),别写成纬度在前,否则采样位置完全错误。如果地貌数据是 Albers 投影而点是 WGS84,需要先通过pyproj转换。

5.3 存储与分享建议:压缩包里的数据该导出成什么格式

当你把加工好的地貌栅格数据分享给同事时,不建议再重新打包成 .rar。因为 RAR 是商业专利格式,很多开源 GIS 环境无法解压。更通用的做法是导出为 GeoTIFF 并使用 LZW 或 DEFLATE 压缩,这样文件可能比原 rar 小,同时不丢失坐标信息:

gdal_translate -of GTiff -co COMPRESS=LZW -co TILED=YES -co BIGTIFF=IF_SAFER 中国地貌.tif 中国地貌_分发.tif

参数说明:TILED=YES让数据按块存取,Web 地图和云服务加载更快;BIGTIFF=IF_SAFER在文件可能超 4GB 时自动切换为 BigTIFF。输出后务必用gdalinfo再验证一次坐标范围,防止.rar里的原始文件本身就没有坐标,压缩包外你需要额外的世界文件.tfw才能定位。没有世界文件时,你只能用gdal_translate -a_srs EPSG:4326 -a_ullr 72 54 135 18手动指定四至范围,这一步如果做错了,后面所有叠加分析都会偏离几公里。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询