江苏省30m DEM地形处理实战：从RAR解压到坡度分类

简介：江苏省地形地貌最新30m精度数据包，面向地理信息、测绘、国土规划与环境研究从业者，提供统一按省整理的tif栅格数据。内容包括海拔分级、起伏程度分类、陆地地貌类型等图层，并附带WGS84与Albers投影坐标参考，便于直接进行地形分析、制图与模型应用。压缩包共30个文件，以tif栅格、tfw坐标信息、dbf属性表及xml元数据为主，另有使用说明，整体仅4.12MB，结构简洁，可快速在ArcGIS等平台中加载使用。目前已有199人学习下载。数据覆盖低海拔至高海拔、丘陵至极大起伏等分级，以及冲积、洪积、海积、冰碛等成因类型，适用于区域地貌对比、灾害风险评价和资源开发规划等场景。

1. 拿到“江苏省地形地貌最新30m精度.rar”之后，先去体检再去解压

一个写着“地形地貌最新30m精度”的RAR包，听起来像是可直接拖进ArcGIS/QGIS出图的DEM，但实际问题是：RAR是压缩格式，里面装的可能是未经分幅处理的全球SRTM、ASTER或ALOS切片，也可能是一套带着坡度、坡向、地貌分类的SHP和TIFF混合包。江苏省大部分地区海拔不足50米，里下河、太湖周边甚至有大面积负地形，垂直误差在30m数据里如果处理不对，会导致淹没模拟和选址分析的结论完全反过来。这篇内容面向GIS工程师和数据开发，先讲怎么判断压缩包里到底是什么源、什么坐标系，再给出用GDAL与Rasterio从裁剪到地形分类的可复现路径，最后用三个技巧在出图前把伪影和坐标偏差拦下来。

2. 解压RAR前先摸清DEM源、文件结构和投影，这步比跑算法更关键

2.1 30m分辨率DEM的主要来源，以及各版本怎么选

市面上说“30m精度”的全球DEM源，常见有SRTM 1弧秒V3、ASTER GDEM V3、ALOS AW3D30 V3.2、Copernicus GLO-30。它们标称分辨率都是约30米，但垂直误差和数据生产方式差异不小。

如果压缩包文件名里没写源，解压前先看内部文件命名。srtm_*就是由NASA按经纬度分块的1度格网；N032E118_DEM开头的是ALOS；带_DSM_COG_前缀的通常是Copernicus。ASTER GDEM在中国东部容易有云影和残余异常，江苏省水体密集，如果包内出现ASTER，建议优先用另外三套做交叉验证。

2.2 用unrar或7z检查RAR包内容，先别急着全量解压

RAR包里的内容通常包含这几种东西：全球或分省DEM的GeoTIFF、行政边界SHP、PRJ投影文件、可能还有一张数据说明PDF或TXT。先用UnRAR列出内部文件，确保没有缺卷或损坏的压缩块。

unrar l 江苏省地形地貌最新30m精度.rar unrar t 江苏省地形地貌最新30m精度.rar unrar x 江苏省地形地貌最新30m精度.rar ./dem_raw/

unrar l是列出清单，能直接看到路径、原始文件大小和压缩后大小，重点看是否包含.prj或.tif.aux.xml。如果只有IMG或BIL格式，说明是国家级基础测绘派生数据，后续还需要看同目录的*_meta.txt。unrar t是完整性测试，只解压校验和，不落盘；这一步能发现RAR卷中损坏的DEM栅格块。最后unrar x才把文件解压到dem_raw/下，避免直接覆盖当前目录。

Windows环境如果没有UnRAR，用7-Zip的命令行版本也可以：7z t 江苏省地形地貌最新30m精度.rar、7z x -odem_raw 江苏省地形地貌最新30m精度.rar。注意GDAL无法直接读取RAR内的GeoTIFF，必须先把DEM文件解压出来再交给后续工具。

2.3 用gdalinfo读元数据，确认坐标系是WGS84还是CGCS2000

解压后用GDAL带上完整路径看一眼栅格元数据，别直接信文件名里的范围。

gdalinfo dem_raw/N032E118_DEM.tif | head -40

输出里重点看四行：Origin、Pixel Size、Coordinate System is、NoData Value。多数全球开源DEM是WGS84经纬度，Pixel Size接近0.0002777778度，也就是1弧秒。如果出现Pixel Size = (30, -30)，那是用投影坐标直接重采样过的30米栅格，坐标系通常是UTM或高斯-克吕格，这时再做坡度计算前要明确使用水平单位还是地理度数。

提示：确认NoData值为负或极大值很重要，比如SRTM用-32768表示水域，Copernicus用0和-9999都出现过，后续裁剪前要把这些值先标记成NaN，否则高程统计会全部被拉偏。

3. 用GDAL与Rasterio对江苏省范围做裁剪、重投影和NoData清洗

3.1 江苏省的经纬度范围与投影选择，避免裁出来是菱形

江苏省大致位于东经116.3°～121.9°、北纬30.7°～35.1°之间，但直接用这个矩形范围裁剪全球DEM会把安徽、山东和浙江的一部分包进来。更稳妥的做法是先准备一份江苏行政边界SHP，用ogr2ogr转换到统一的EPSG后再做gdalwarp -cutline。

为了后续量算坡度、面积和距离，建议把经纬度坐标重投影到CGCS2000 3度高斯-克吕格投影。江苏横跨中央经线120°和121.5°两带，省级制图常见做法是统一用中央经线120°的3度带，EPSG为4521（CGCS2000 / 3-degree Gauss-Kruger zone 39）。实际使用中如果只做省级宏观展示，也可以用EPSG:4526或WGS84 / UTM 50N，但30m地形因子计算建议保留高精度的CGCS2000投影，避免在边界处发生米制接边偏差。

3.2 用gdalwarp按行政边界裁剪为江苏省DEM

先把边界转为与DEM一致的坐标系，再用gdalwarp配合-cutline进行裁剪，同时完成重投影和NoData替换。

ogr2ogr -t_srs EPSG:4521 jiangsu_4521.shp jiangsu_boundary.shp gdalwarp -t_srs EPSG:4521 \ -cutline jiangsu_4521.shp -crop_to_cutline \ -tr 30 30 -r cubic \ -dstnodata -9999 -dstalpha \ dem_raw/N032E118_DEM.tif dem_raw/N032E119_DEM.tif \ jiangsu_dem_30m.tif

-tr 30 30表示输出栅格像素尺寸为30米×30米；-r cubic对DEM重采样用三次卷积，比bilinear平滑且不会像nearest那样出现台阶感，但不要在坡度分析后再插值。-dstnodata -9999统一NoData值，方便后续使用Rasterio处理。-dstalpha生成透明波段，能把江苏边界外的像素全部变透明，避免后期用大范围0值做无效数据。多景输入时gdalwarp会自动按地理位置拼接，但要求所有输入DEM的像素尺度和坐标系一致，所以这里先做了-t_srs统一。

注意：如果RAR包内已经是一份分幅好的江苏省DEM，就不要再叠加cutline，否则需要先合并再裁剪，白白增加计算量。

3.3 用Rasterio做NoData清洗与特殊海拔修正

即使DEM元数据里有NoData值，实际读取时也会遇到三类“脏数据”：填充负值、异常极大值、湖面高程忽高忽低。下面用Rasterio读取STAC边界后的江苏DEM并做处理。

import numpy as np import rasterio from rasterio.fill import fillnodata from rasterio.warp import reproject, Resampling src_dem = "jiangsu_dem_30m.tif" out_dem = "jiangsu_dem_filled.tif" with rasterio.open(src_dem) as src: dem = src.read(1).astype("float32") profile = src.profile.copy() nodata = src.nodata transform = src.transform # 实际现象：局部水域高程为-9999，需要填为周围陆地高程 mask = np.isclose(dem, nodata) | (dem < -100) | (dem > 1500) dem[mask] = np.nan # 对江苏区域，使用邻域插值填充无效像元 filled = fillnodata(dem, mask=~np.isnan(dem), max_search_distance=20) # 将高程异常超过5倍中位数的像素视为粗差 median_elev = np.nanmedian(filled) filled[np.abs(filled - median_elev) > 5 * np.nanstd(filled)] = median_elev profile.update(dtype="float32", nodata=np.nan) with rasterio.open(out_dem, "w", **profile) as dst: dst.write(filled.astype("float32"), 1)

这段代码先把NoData和异常低值全部转成NaN，再用GDAL内置的fillnodata按最大20像素距离做近邻插值。江苏地面高程整体低于1500米，把大于1500m的像素直接视作异常值处理，主要是压制拔地而起的融合伪影。之后再用中位数替代远离统计分布的粗差点，这比直接置0更安全，也方便后续坡度计算。

参数说明里最容易忽略的是max_search_distance。该参数决定空洞附近搜索有效栅格的最大半径，单位是像素。江苏省密集水网区域可能让30m栅格上的空洞连成片，如果设得太小会保留洞，设太大又会让山体边缘变平滑。建议先用gdal_fillnodata.py配合-md 20做一次，若还是有残留，再看是不是湖心连续空洞超过20像素，需要把max_search_distance提高到50。

3.4 与全国或全球模型不一致时的垂直基准对齐

如果RAR包内数据来自ASTER或Copernicus DSM，地物高度超过DEM，尤其在沿江建筑和高铁桥梁处，和实测水准点比对会偏高3～8米。这时不能直接修改高程值，而是准备一批高精度控制点，用gdalwarp做一次带多控制点的多项式拟合。

gdalwarp -tps -tr 30 30 -r cubic -dstnodata -9999 \ -co COMPRESS=DEFLATE -co PREDICTOR=2 \ jiangsu_dem_filled.tif jiangsu_dem_warp.tif

4. 从DEM提取坡度、坡向、地形起伏度，并生成山体阴影底图

4.1 gdaldem参数速查表，按需设置水平单位

gdaldem生成的坡度坡向是对每个像素与周围8邻域计算二阶差分，前提是DEM已经是投影坐标系。江苏省用的EPSG:4521是米制坐标，所以坡度公式可以直接使用米制水平距离，不需要再做纬度系数修正。下面是常用参数。

4.2 江苏大面积缓坡地区的坡度要使用百分数坡度

江苏省大部分平原坡度小于1度，直接输出弧度过小且无法分级，适合使用-p输出坡度百分比。示例命令：

gdaldem slope jiangsu_dem_filled.tif jiangsu_slope_pct.tif -p -of GTiff gdaldem aspect jiangsu_dem_filled.tif jiangsu_aspect.tif -zero_for_flat gdaldem hillshade jiangsu_dem_filled.tif jiangsu_hillshade.tif -az 315 -alt 45 -z 1.4

用-p后，平缓地区坡度值在0～5之间，分级标准可改成：0.5%以内为平地，0.5～5%为缓坡，5～15%为中等坡，大于15%为陡坡。输出栅格使用32位浮点，后续渲染可 直接拉伸拉伸到0～255。-az 315 -alt 45是常见多方向山体阴影中西北光源的参数，能增强苏南丘陵的地形纹理，但要注意海岸堤防等线性地物会留下条状阴影，并不代表真实的高程变化。

4.3 用地形位置指数把地貌分类成平地、丘陵和山地

单纯海拔阈值无法区分江苏的宁镇山脉和废黄河高地，更常用的方法是计算地形位置指数TPI，也就是某点高程与周围环形邻域平均高程的差。

gdaldem TPI jiangsu_dem_filled.tif jiangsu_tpi.tif -p 5 5

参数-p 5 5定义了一个外半径5像素、内半径0像素的矩形环，大致对应水平距离150米的分析窗口。TPI等于0表示该点与周边地面等势，大于0是山顶或山脊，小于0为谷底。随后在QGIS或Python中按下表分类。

江苏没有真正的高山，山地区域主要用于太湖宜兴山区和连云港云台山周边，分类阈值可以参考。生成TIF后可以直接用QGIS的Raster calculator把多个因子叠加：

gdal_calc.py -A jiangsu_dem_filled.tif -B jiangsu_tpi.tif -C jiangsu_slope_pct.tif \ --outfile=jiangsu_landform_class.tif \ --calc="(abs(B)<1)*(C<0.3)*1 + (abs(B)<1)*((C>=0.3)*(C<2))*2 + (A>30)*(B>3)*3 + (A>200)*(C>=10)*4"

4.4 把山体阴影和坡度图叠成地形底图

最终出图推荐用QGIS加载jiangsu_hillshade.tif作为灰度底图，上面叠加带半透明的坡度图。山体阴影用multiply混合模式，坡度图用overlay混合模式。GDAL也可以直接把两种栅格融合成一个RGB：

python make_terrain_rgb.py

使用Python和Rasterio逐波段归一化后写入三通道GeoTIFF，并设置与DEM相同的地理变换，这样一套数据就能在OpenLayers或Leaflet里切片发布。

5. 检验30m成果的三个实用技巧，避免出图后返工

5.1 验证原始RAR中的高程值有没有被重新投影破坏

最快的检验方式是用gdalinfo -stats看最小值、最大值和标准差。江苏陆域正常海拔应在-20～600米之间，如果最大值出现在2000米以上，说明存在未清洗的异常像元；如果最小值是-32768，说明NoData没被当成无效值参与统计。再用gdal_calc.py统计异常值比例，超过总像元数0.5%时就该回头重做第3章的空洞填充。

5.2 用一条剖面线对比SRTM、ALOS与Copernicus的高程折线

在QGIS横断面插件或Python中沿长江南通—南京段画一条剖面线，对三套DEM分别取样。

import rasterio import numpy as np files = ["copernicus_30m.tif", "alos_30m.tif", "srtm_30m.tif"] x = np.linspace(120.1, 118.7, 500) y = np.linspace(32.0, 32.1, 500) for fp in files: with rasterio.open(fp) as src: coords = [(xx, yy) for xx, yy in zip(x, y)] vals = [v[0] for v in src.sample(coords)] print(fp, np.nanmedian(vals), np.nanmin(vals), np.nanmax(vals))

如果Copernicus在长江大桥附近明显高出SRTM数米，那是DSM与DEM的差异，不是错误。相反，如果ALOS在低洼圩区出现锯齿状阶梯，说明该源数据在平坦区域回波噪声过大，建议最终成果以Copernicus辅以SRTM填充为准。

5.3 用行政边界做最后的蒙版并导出成金字塔GeoTIFF

出图前把江苏省边界再次作为cutline，使用gdalwarp -dstalpha生成带透明边的成果，再用gdaladdo构建金字塔：

gdalwarp -cutline jiangsu_boundary.shp -crop_to_cutline -dstalpha \ jiangsu_dem_filled.tif jiangsu_dem_release.tif gdaladdo -r average -ro jiangsu_dem_release.tif 2 4 8 16

带-ro标识的内部金字塔能让QGIS和GeoServer加载大范围地形图时快速缩放，同时不会影响原始像素值。最后对比gdalinfo -stats中的像元数量与江苏总面积，若有效像元数乘以900平方米（30m×30m）得到的面积偏离江苏陆域面积超过5%，则说明裁剪或重投影过程中出现了像元错位，需要从第3章的-tr 30 30和-t_srs重新检查。

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