地形起伏度RdlsChina1km数据集：从原理到栅格读取与建模应用

简介中国地形起伏度公里网格数据集RdlsChina1km是一套基于SRTM 90米数字高程模型重采样至1公里分辨率并通过专门模型计算得到的全国陆地地形起伏度栅格与统计资料主要面向地理信息系统、地貌学、生态环境与区域规划等领域的研究者和数据分析人员可支持大尺度地形特征提取、分区对比、模型建模及可视化制图等应用。包内既包含中国全国和分省的地形起伏度公里网格空间数据也提供分省、分地区、分县的统计数据兼顾空间分析与表格查询两类需求。资源共522个文件以ARCGIS GRID格式的adf、nit、dat等栅格文件为主辅以log操作日志、xml元数据文件以及一个xlsx格式的统计汇总表压缩包整体约15.34MB结构清晰便于按目录调用与归档。该数据集可直接导入ArcGIS等专业平台进行显示、掩膜提取与栅格计算也可利用Excel快速查看不同行政区域的起伏度数值为相关研究提供基础数据支撑目前已有653人学习下载适合需要获取中国地形起伏度底图或构建区域模型数据的用户参考使用。1. 地形起伏度不能只看高差RdlsChina1km 帮你把“地形破碎度”也算进去同样是高海拔山区很多人第一反应是去查最高点与最低点的海拔差。可到了藏北高原海拔超过四千米但地势开阔平缓到了横断山区的怒江峡谷平均海拔不算极端几十公里内高差却超过两千米地形对人的活动限制完全不是一个量级。地形起伏度Relief Degree of Land SurfaceRDLS正是把“高度差有多大”和“平地占比有多高”同时纳入计算的综合指标比单纯的海拔标准差或坡度更贴近地形对工程建设、灾害易发性和人口分布的约束。RdlsChina1km中国地形起伏度公里网格数据集是一个已经按 1km×1km 公里网格生成好的全境栅格产品下载下来通常是一个 zip 包里面常见的是一张 GeoTIFF 文件。打开后每个像元的值不再是海拔而是该格网内的起伏度数值。它高频出现在人口密度估算、地震地质灾害评价、电网与公路选线、生态红线划定这几类工作里适合 GIS 开发者和做地学建模的工程师直接拿去当特征图层。2. 先理解 RdlsChina1km 的算法逻辑起伏度不是坡度也不是高程方差2.1 国内常用 RDLS 公式里的三个关键变量地形起伏度的计算在不同机构的生产流程里存在多个变体但国内万级格网产品中最常见的是把“最高海拔与最低海拔之差”和“平地面积占比”组合起来。以单个 1km 网格为例常见形式为RDLS [ Max(H) - Min(H) ] × [ 1 - P(A) / A ] / 500其中Max(H) - Min(H)是这个网格内的最大最小海拔差单位是米A是网格总面积P(A)是网格内坡度小于等于 5° 的平地面积。分母 500 的作用是把数值压缩到便于制图和分级的量级典型输出范围在 0 到 200 之间数值越大代表地表越破碎、越不适合大规模工程建设。注意这个公式同时惩罚两种地形高差大的山地和虽然没有高差但几乎没有平坦土地的破碎丘陵后者正是单看海拔方差容易漏掉的部分。2.2 为什么用 1km 网格而不是原始 DEM 分辨率SRTM 或 ASTER GDEM 的原始分辨率是 30m 到 90m直接拿原始像元计算结果会被单条冲沟或陡崖放大形成大量噪声。RdlsChina1km 的做法是先对高精度 DEM 做坡度计算和填洼处理再在目标网格尺度上做统计。1km 格网起到空间平滑作用让一个县域内部的起伏度呈现连续渐变而不是像 30m 数据那样打成碎斑。这里要格外注意一个反直觉结论网格越大RDLS 数值反而越小因为窗口内最高最低点的极差增速远低于平地面积的增速。所以当你拿到数据集后不要用 90m 分辨率产品常引用的 30-70-120 那套分级阈值需要先看这个数据集自身的分位数分布。2.3 与海拔标准差、坡度、地形粗糙度的选型差异指标计算成本对“平地碎部”的敏感度典型使用场景海拔标准差低低宏观地貌分区坡度低高但会把单像元噪声放大坡耕地识别、滑坡面提取地形起伏度 RDLS中中高空间平滑后更稳人口分布建模、灾害易损性、选线地形粗糙度中中依赖窗口大小地表径流模拟、地表粗糙度参数化实际使用中我一般会把坡度当作“局部指标”把 RDLS 当作“区域指标”。如果做省级尺度的地质灾害易发性评价坡度适合做孕灾因子RDLS 适合做承灾体暴露度的修正因子。如果你发现模型的变量共线性过高优先保留 RDLS 而不是同时放进坡度和海拔标准差因为它已经包含了一部分两者的信息。3. 从 zip 到可用栅格用 GDAL 和 rasterio 完成读取与坐标核对3.1 解压后用 gdalinfo 先看元数据别急着写代码拿到RdlsChina1km.zip之后第一步不是直接扔进 Python而是先用命令行工具确认坐标系、位深、无数据值和压缩格式。很多坑在这一步就能提前暴露。mkdir -p rdls_data unzip RdlsChina1km.zip -d rdls_data gdalinfo rdls_data/RdlsChina1km.tif执行后重点看四行输出Size给出栅格宽高应为约 4900×3300 级别Coordinate System确认是经纬度还是投影坐标NoData Value记录空值的具体数值Compression查看是否已经做过 LZW 或 DEFLATE 压缩。如果 zip 里面除了 tif 还有.prj或.tfw文件说明数据附带 ESRI 风格投影定义和世界文件它们能够辅助判断生产方用的是 CGCS2000 还是 WGS84这两者在 1km 格网上的偏差虽然只有几十米但叠加入口普查边界数据时会产生明显的边缘错位。3.2 用 rasterio 按坐标提取任意点的起伏度值读 GeoTIFF 最稳的组合是rasterio加上numpy。下面的代码演示如何用经纬度坐标反算像元行列号再只读取那一个小窗口避免把整张栅格载入内存。全图约上亿像元直接read(1)会占用大几百 MB 内存不优雅。import rasterio tif_path rdls_data/RdlsChina1km.tif with rasterio.open(tif_path) as src: print(CRS:, src.crs) print(NoData:, src.nodata) print(像元尺寸:, src.width, x, src.height) # 给定一个北京附近的经纬度点 lon, lat 116.38, 39.90 # index 方法做的是仿射变换的逆运算返回 (row, col) row, col src.index(lon, lat) print(f行列号: row{row}, col{col}) # window 只读目标像元所在的一格而不是整幅影像 val src.read(1, window((row, row 1), (col, col 1))) print(该点起伏度:, val[0, 0])这里需要说明两个参数的含义。src.crs决定了index()方法内部如何处理坐标如果你的输入点是 WGS84 经纬度而栅格是投影坐标必须先调用src.transform完成变换再取行列号否则取出来的点会偏移几个格网。src.nodata通常是 -9999 或 0如果某个像元落在境外或者原始 DEM 缺失区域读出来的值就是这个标记把它当成真实起伏度参与建模会造成极大误差。3.3 批量提取多个点时用 transform 避免重复开文件如果只是验证几个点每点单独打开文件可以接受但要做抽样验证时我一般会把 shp 里的点一次性读进来用rasterio.features.geometry_mask或直接调src.sample()完成批量采样。import geopandas as gpd import rasterio points gpd.read_file(sample_points.shp) with rasterio.open(tif_path) as src: # sample 接受迭代器内部按行列批量取像元值比 for 循环逐个 index 快得多 coords [(geom.x, geom.y) for geom in points.geometry] values [v[0] for v in src.sample(coords)] points[rdls_value] values points.to_file(sample_points_rdls.shp)sample()方法的一个隐含风险是它不做边界检查落在栅格范围之外的点会返回nodata值所以输出后需要value ! nodata过滤一遍。这样批量处理几千个点基本是毫秒级足够支撑一次快速数据质检。4. 建模型前先做四件事裁剪、重投影、分区统计与分位数重分类4.1 用 gdalwarp 把全国数据裁剪到研究区并统一坐标系绝大多数分析不会直接用全国范围而是先切到省、市或某个流域。这里常见做法是用gdalwarp配合矢量边界完成裁剪与重投影一步到位。gdalwarp -cutline study_area.shp -crop_to_cutline \ -t_srs EPSG:3857 -r bilinear -of GTiff \ rdls_data/RdlsChina1km.tif rdls_study.tif参数说明-cutline指定矢量裁剪边界-crop_to_cutline让输出栅格的范围严格贴住边界而非外接矩形-t_srs把坐标系统一到 Web Mercator 或你后续建模使用的投影这里 EPSG:3857 只是一个示例如果做面积统计更推荐用 Albers 等积投影-r bilinear指定重采样算法RDLS 是连续浮点变量用双线性插值比用最近邻更平滑但如果后续要做类别比较建议改回-r near以保留原始数值。裁剪完成后用gdalinfo再确认一下像元尺寸因为重投影后 1km 格网在不同的纬度带会略微变形。4.2 用 zonal_stats 快速统计县域尺度的平均起伏度接下去最常见的需求是把栅格聚合到行政区边界。rasterstats库的zonal_stats函数是这里最顺手的工具。import geopandas as gpd from rasterstats import zonal_stats counties gpd.read_file(counties.shp) # 一次统计均值、最大值和标准差用于下游特征合并 stats zonal_stats( counties, rdls_study.tif, stats[mean, max, std, count], geojson_outTrue, nodata-9999 ) for feature in stats: feature[properties][rdls_mean] feature[properties].pop(mean) feature[properties][rdls_max] feature[properties].pop(max) gpd.GeoDataFrame.from_features(stats).to_file(counties_rdls.shp)这里的stats参数列表决定了输出的聚合字段count可以用来检查每个县域内有效像元的数量如果某个县的count远小于理论值多半是裁剪时出现了大块 NoData 区域。nodata参数如果不显式传rasterstats会读取 tif 的元数据但一旦你前面重投影时把 NoData 值改了就必须手动指定否则统计均值会被 -9999 严重拉低。另一个提示不要先对整个全国栅格做全量排序再查前几名那种做法既慢又浪费内存用分区统计先聚合再对结果做 pandas 排序就足够了。4.3 用分位数重分类代替固定阈值避免跨数据集误用1km 格网的 RDLS 数值整体会低于 30m 产品直接用学术文献里的固定阈值比如 30、70、120很容易把全国都分进“低起伏度”一类。稳妥做法是先看研究区的分位数再按分位数定义低、中、高、极高四类。import rasterio import numpy as np with rasterio.open(rdls_study.tif) as src: data src.read(1) valid data[data ! src.nodata] # 取四分位数作为切分点 q25, q50, q75 np.percentile(valid, [25, 50, 75]) print(f分位数切分点: {q25:.2f}, {q50:.2f}, {q75:.2f}) classified np.zeros_like(data, dtypenp.uint8) classified[(data 0) (data q25)] 1 classified[(data q25) (data q50)] 2 classified[(data q50) (data q75)] 3 classified[(data q75)] 4 classified[data src.nodata] 0这样分级的好处是自动适应当前研究区的数据分布。如果研究区整体是平原分位数会把细微起伏也区分出来如果是横断山区分位数会让强起伏区内部进一步分层不会出现一个省几乎全落在最高级的情况。分位数切分点参考下表目标用途推荐分级方式说明人口分布建模3 级低/中/高均值聚合后按自然断点地质灾害易发性4 级强调极端区分位数 最大值辅助修正生态红线与保护区评价5 级保留渐变分位数 空间平滑后再分电网选线2 级可建设/不可建设阈值取研究区 90 分位数np.percentile计算时如果数组里含 NaN需要先用np.isnan过滤这里使用valid data[data ! src.nodata]做了脏值隔离是更通用的写法。5. 验证与排错NoData 掩膜、压缩格式和边缘效应怎么处理栅格数据的排错和写代码一样需要“单测”思路固定一组已知地理位置的点对比采样值和真实地形认知再逐步排查系统性问题。第一件事是确认 NoData 值是否被正确识别。某些数据源会把无数据区写成 0这会造成沿海和境外区域出现“平地”假象。验证方法是统计全图的最小值和 0 值像元占比如果最小值为 0 且分布集中在国界外侧建议把 0 统一重映射成 -9999 再继续后续操作。gdal_calc.py -A rdlsChina1km.tif --outfilerdls_nodata.tif \ --calcA*(A0) --NoDataValue-9999第二件要检查的是文件压缩格式。市面上大量辅助工具对 tif 的兼容性取决于压缩算法LZW 是无损压缩且兼容度最高如果解压后发现 tif 是 DEFLATE 压缩而某个建模框架读不了先用gdal_translate -co COMPRESSLZW转换一次。第三件容易被忽略的是边缘效应。以像元为中心的 1km 格网在国界、省界处会跨越边界取到邻域数据导致边界像元值被境外地形拉高或拉低。如果研究区靠近边境对比counties_rdls.shp中边界县和内陆县的均值分布凡是标准差明显偏大的边界县都要单独标记。更细的验证方法是在原图上取一个 5km×5km 的窗口用公式重算一次 RDLS与数据集自带值比对误差超过 10% 说明窗口定义的起点存在偏移。最后建议在输出成果时把 NoData 掩膜单独导出为一张二值 tif而不是靠数值约定传递。多数后续做机器学习的同事不会去读头文件里的 NoData 标记直接把掩膜作为额外输入特征或裁剪条件能省掉下游一大半的排查时间。本文还有配套的精品资源点击获取