简介江苏省地形地貌最新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米但垂直误差和数据生产方式差异不小。数据源发布时间线覆盖范围垂直精度中误差(m)典型文件名SRTM 1 V32015以后重新整理北纬60°南纬56°约4~9srtm_58_05.tifASTER GDEM V32019年发布V383°N83°S约7~14ASTGTM_N32E118.tifALOS AW3D30 V3.22021年后更新82°N82°S约4.5N032E118_DEM.tifCopernicus GLO-302022年公开Global 30m全球约3.5Copernicus_DSM_COG_10_N32_00_E118_00_DEM.tif如果压缩包文件名里没写源解压前先看内部文件命名。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为4521CGCS2000 / 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_distance20) # 将高程异常超过5倍中位数的像素视为粗差 median_elev np.nanmedian(filled) filled[np.abs(filled - median_elev) 5 * np.nanstd(filled)] median_elev profile.update(dtypefloat32, nodatanp.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尤其在沿江建筑和高铁桥梁处和实测水准点比对会偏高38米。这时不能直接修改高程值而是准备一批高精度控制点用gdalwarp做一次带多控制点的多项式拟合。gdalwarp -tps -tr 30 30 -r cubic -dstnodata -9999 \ -co COMPRESSDEFLATE -co PREDICTOR2 \ jiangsu_dem_filled.tif jiangsu_dem_warp.tif4. 从DEM提取坡度、坡向、地形起伏度并生成山体阴影底图4.1 gdaldem参数速查表按需设置水平单位gdaldem生成的坡度坡向是对每个像素与周围8邻域计算二阶差分前提是DEM已经是投影坐标系。江苏省用的EPSG:4521是米制坐标所以坡度公式可以直接使用米制水平距离不需要再做纬度系数修正。下面是常用参数。功能命令片段关键参数坡度gdaldem slope input.tif slope.tif -p -s 111120-p百分数坡度-s比例系数若是米制投影可不设坡向gdaldem aspect input.tif aspect.tif -zero_for_flat平地输出0坡向角度0360山体阴影gdaldem hillshade input.tif hillshade.tif -az 315 -alt 45 -z 1.4光线方向、高度角、垂直夸大系数地形起伏度gdaldem tri input.tif tri.tif计算地形粗糙度单位与高程单位一致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后平缓地区坡度值在05之间分级标准可改成0.5%以内为平地0.55%为缓坡515%为中等坡大于15%为陡坡。输出栅格使用32位浮点后续渲染可 直接拉伸拉伸到0255。-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中按下表分类。高程条件TPI条件坡度条件分类结果任意介于-1和1之间坡度P 0.3%平地/水域任意介于-1和1之间0.3% ≤ P 2%微起伏平原海拔30mTPI 3任意丘陵海拔200mTPI 5P ≥ 10%低山任意TPI -3任意山谷/洼地江苏没有真正的高山山地区域主要用于太湖宜兴山区和连云港云台山周边分类阈值可以参考。生成TIF后可以直接用QGIS的Raster calculator把多个因子叠加gdal_calc.py -A jiangsu_dem_filled.tif -B jiangsu_tpi.tif -C jiangsu_slope_pct.tif \ --outfilejiangsu_landform_class.tif \ --calc(abs(B)1)*(C0.3)*1 (abs(B)1)*((C0.3)*(C2))*2 (A30)*(B3)*3 (A200)*(C10)*44.4 把山体阴影和坡度图叠成地形底图最终出图推荐用QGIS加载jiangsu_hillshade.tif作为灰度底图上面叠加带半透明的坡度图。山体阴影用multiply混合模式坡度图用overlay混合模式。GDAL也可以直接把两种栅格融合成一个RGBpython make_terrain_rgb.py使用Python和Rasterio逐波段归一化后写入三通道GeoTIFF并设置与DEM相同的地理变换这样一套数据就能在OpenLayers或Leaflet里切片发布。5. 检验30m成果的三个实用技巧避免出图后返工5.1 验证原始RAR中的高程值有没有被重新投影破坏最快的检验方式是用gdalinfo -stats看最小值、最大值和标准差。江苏陆域正常海拔应在-20600米之间如果最大值出现在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重新检查。本文还有配套的精品资源点击获取