吉林市30米DEM数据处理全流程：从解压到坡度坡向与填洼分析

简介：吉林省吉林市30米分辨率DEM数字高程数据包，面向GIS学习者、城乡规划、环境分析及地质评估等从业者，可精准反映地表起伏，用于地形建模、遥感配准、灾害研判等场景。压缩包共12个文件，涵盖核心TIFF高程栅格、范围Shapefile及其配套投影参数（.prj）、属性库（.dbf）、空间索引（.sbn/.sbx）与坐标元数据等，整体约81.88MB，解压后可直接在ArcGIS、QGIS中加载分析。当前已有712人学习使用，适合需要高精度市级地形底图、开展空间统计或制作专题地图的用户。内容丰富且结构完整，省去自行拼接边界与坐标转换的步骤，能快速获得吉林市域连续高程数据与行政范围图层，提升地理信息处理效率。

1. 拿到这份DEM压缩包，先别急着解压

吉林市地处松花江中上游，城区沿江而建，东西两侧分别是低山丘陵和松嫩平原边缘，地形高差从接近200米的河谷到一千多米的山脊都有分布。做防洪淹没模拟、地质灾害易发性评价、基站选址或者光伏阴影分析时，30米分辨率的DEM是性价比很高的起点。这份数据之所以是zip格式，是因为里面同时装了两类东西：一类是30m的DEM栅格文件，另一类是吉林市本级的行政区划范围shp。后者看似只是顺手附带的边界框，实际价值并不小，国内公开渠道能直接下载到地级市精确边界的渠道很少，很多项目卡在边界上。本文以这份数据为例，把解压检查、坐标系处理、按shp裁剪以及坡度坡向、填洼、等高线这些后续产出一次梳理到位，配置环境和命令都直接给出，可以在ArcGIS、QGIS和Python环境里对照着用。

2. 解压后先做三件事：验包、看结构、读分辨率

2.1 zip文件完整性检查和文件清单

下载下来的zip先别急着双击，用命令行验一遍包是最稳妥的。Windows下可以用系统自带tar命令，也可以装7-Zip后使用其命令行工具；Linux和macOS直接用unzip即可。

unzip -t 吉林省吉林市DEM数字高程数据30m（含本市级范围shp文件）.zip

这段命令会逐个测试包内文件的CRC32校验值并输出“No errors detected in compressed data”之类的结论。要是网络传输过程中丢包，解压时容易遇到“unexpected end of file”或“CRC failed”，问题不是出在压缩工具上，而是源文件本身不完整，重下比修复快得多。

验包没问题后，列出zip内部结构，确认栅格和矢量文件都齐全再解压：

unzip -l 吉林省吉林市DEM数字高程数据30m（含本市级范围shp文件）.zip

-l参数只列出内容不解压，能看到文件名、原始大小和压缩后大小。一个成熟的DEM压缩包里，栅格部分应该有一个tif或img主文件（30米分辨率可能配一个tfw世界文件或prj投影文件），shp部分则是一组配套文件。如果只有孤零零一个shp而无dbf和shx，说明打包方漏文件了，这类shp在很多GIS软件里根本打不开。

解压建议单独建目录，不要直接撒在桌面或下载文件夹里，后续用Python或ArcGIS处理时，纯英文绝对路径能省掉大量编码问题。

mkdir -p /data/jilin_dem && unzip 吉林省吉林市DEM数字高程数据30m（含本市级范围shp文件）.zip -d /data/jilin_dem

2.2 shp不是单个文件，是一组配套文件的合称

很多第一次接触矢量数据的同学会对着“.shp”找文件，其实一个完整的shapefile至少包含三个基础文件，完整形态更多。这份数据既然标注“含本市级范围shp文件”，理应带齐下面这些。

查看文件清单时，如果.prj存在，先打开看一眼坐标系，这决定后面所有处理步骤。如果缺.prj，可以用ArcGIS的Define Projection工具补上，但前提是你得知道数据原本的坐标系，瞎补一个会导致整个空间分析全军覆没。

2.3 用gdalinfo或rasterio读出真正的分辨率

“30m”是标题里的标签，实际数据到底是什么分辨率，不能光看名字。打开栅格元数据看一眼才算数。

gdalinfo dem.tif

重点看这几行输出：

Size is 12412, 10241 Coordinate System is: PROJCRS["CGCS2000 / 3-degree Gauss-Kruger CM 126E", ... Pixel Size = 25.000000000000000, -25.000000000000000

Size是栅格的像元行列数，Pixel Size是单个像元在地面的大小，单位跟随坐标系定义。如果Pixel Size是0.00027这种带小数的值，说明数据还是经纬度坐标；此时虽然水平间隔大约对应30米，但不同纬度上这个值对应的地面距离并不一致，做面积、距离计算之前必须投影到平面坐标系。

用Python检查也可以，适合要写脚本批量处理多幅数据或加入自动化流程的情况。rasterio是当前最主流的栅格读取库，直接把关键的元数据字段打印出来。

import rasterio with rasterio.open("/data/jilin_dem/dem.tif") as src: print("像元尺寸:", src.res) print("坐标系:", src.crs) print("数据范围:", src.bounds) print("无效值:", src.nodata) print("波段数:", src.count) print("数据类型:", src.dtypes[0])

src.res给出x和y方向的像元大小，如果打印结果是(30.0, 30.0)，说明这确实是真正的30米分辨率，而且大概率已经处于投影坐标系中；如果出现(0.00027, 0.00027)之类的小数，后续需要做投影转换才能进入面积量算。src.nodata这个值也很重要，本数据可能用-9999或-32767表示无数据区域，裁剪和统计分析时如果不把nodata排除在外，算出来的坡度平均值和地形起伏度都会是错的。

3. 坐标系是DEM和shp能不能叠到一起的前提

3.1 用.prj识别数据自带坐标系

DEM和shp能叠到正确位置，前提是两个数据的坐标系一致，或者至少能在软件里动态转换。吉林市地处东经125°40′到127°56′、北纬42°31′到44°40′之间，本地项目常用坐标系集中在两种：一种是CGCS2000或WGS84下的经纬度坐标，另一种是高斯-克吕格投影，中央经线多为126°E或129°E。打开shp的.prj文件，能看到类似“GCS_WGS_1984”或“CGCS2000_3_Degree_GK_CM_126E”的字样，这就直接说明了边界数据的空间参考。

一个很常见的坑是：DEM是WGS84经纬度坐标，shp是CGCS2000高斯投影，两者虽然都是“2000系”，但一个在度上、一个在米上，直接扔进ArcGIS会提示地理坐标系与投影坐标系不匹配。ArcGIS虽能在显示层面实时投影对齐，但做重采样、面积统计这类操作时必须先把两者统一到同一坐标系。

也可以在Python里读一下两个文件的crs，快速判断是否需要转换：

import geopandas as gpd import rasterio shp = gpd.read_file("/data/jilin_dem/吉林市本级.shp") with rasterio.open("/data/jilin_dem/dem.tif") as src: print("shp坐标系:", shp.crs) print("DEM坐标系:", src.crs)

3.2 30米在经纬度和投影坐标系下的差别

为什么说“30m”这个说法在经纬度坐标里是近似值？地球不是正球体，经线在赤道最疏、在两极汇于一点。吉林市所在纬度约北纬43度，1度经度对应的地面距离大约是81公里，而1度纬度约111公里。如果一个DEM标称30米但仍是经纬度坐标，它实际上只在赤道附近是严格30米，到吉林市已经退化成一个长条形网格，每个像元的地面宽度小于30米。直接在这种数据里做坡度计算，x和y方向的距离基准不一致，坡度和坡向结果都会带偏差。

统一投影时，优先选择高斯-克吕格投影或Albers等积投影，而不是Web墨卡托。Web墨卡托在高纬度地区面积变形巨大，吉林市离赤道超过4700公里，用Web墨卡托做面积统计会明显偏大，适合做底图展示，不适合做地形量算。

3.3 用gdalwarp把DEM重投影到shp所在坐标系

判断完两者的坐标系后，以shp所在的平面坐标系为基准统一数据。用gdalwarp最直接，一个命令完成重投影加重采样。EPSG代码需要从shp.prj里查出后替换，吉林市常见的情况是4490（CGCS2000经纬度）转4547或4548（CGCS2000高斯-克吕格3度带）。

gdalwarp -t_srs EPSG:4547 -tr 30 30 -r bilinear -of GTiff dem_wgs84.tif dem_cgcs2000_30m.tif

-t_srs指定目标坐标系，-tr 30 30强制输出像元为30米×30米，-r bilinear用双线性插值重采样。双线性对地形连续表面的DEM是比较合理的，因为它比最邻近法平滑，又不像三次卷积那样可能把高程值拉出异常尖峰。如果后续做水文分析，填洼之前要保留真实地形细节，用cubic可能会制造伪洼地，这个细节后面还会提到。

4. 用shp裁剪DEM并生成坡度、坡向与水文分析数据

4.1 先统一shp与DEM到同一个坐标系

重投影只是把DEM转过去了，shp本身如果也是WGS84或其他坐标系，最好也转换成和DEM完全一样的EPSG，避免后续每次调用都隐式转换。用ogr2ogr完成：

ogr2ogr -t_srs EPSG:4547 吉林市本级_4547.shp 吉林市本级.shp

如果需要了解坐标转换的差异量级，可以给shp加两个字段记录转换前后的质心坐标，但日常场景没这个必要。重点在于，转换后务必在GIS里叠加查看一次，看shp边界是不是和DEM的地形特征对齐——比如河流位置是否落在山谷线上，如果偏移明显，多是被错误投影或错误带号害的。

4.2 用rasterio提取shp范围内的DEM

裁剪方式有两种，一种是直接裁剪成shp外接矩形范围的规则矩形，另一种是精确到shp边界的带透明区域的裁剪。前者适合模型输入，后者适合出图和生产标准成果。

用Python的rasterio实现精确裁剪，代码不长，但nodata处理是核心：

import rasterio from rasterio.mask import mask import geopandas as gpd shp = gpd.read_file("/data/jilin_dem/吉林市本级_4547.shp") geom = [shp.geometry.unary_union] # 全部面要素合并为一个几何 with rasterio.open("/data/jilin_dem/dem_cgcs2000_30m.tif") as src: out_img, out_transform = mask( src, geom, crop=True, nodata=-9999, all_touched=False ) with rasterio.open( "/data/jilin_dem/dem_jilin_clip.tif", "w", driver="GTiff", height=out_img.shape[1], width=out_img.shape[2], count=1, dtype=out_img.dtype, crs=src.crs, transform=out_transform, nodata=-9999, ) as dst: dst.write(out_img, 1)

all_touched=False意味着只有像元中心落在shp范围内才保留，边缘会更符合真实边界，但会导致边界处少量像元被裁掉；all_touched=True则把任何与边界相交的像元都保留，边缘锯齿感低一些，面积会略大于真实值。统计面积时用False，出地形图时用True，这两个场景的取向本来就是矛盾的。

用QGIS或ArcGIS的裁剪工具原理相同，本质都是gdalwarp -cutline。命令行版本更利于批量操作和生成日志：

gdalwarp -cutline 吉林市本级_4547.shp -crop_to_cutline -dstnodata -9999 -of GTiff dem_cgcs2000_30m.tif dem_jilin_clip.tif

-dstnodata -9999把裁剪后落在shp外的区域全部设为-9999，后续在arcpy或numpy统计时需要单独排除此值。如果裁剪后发现边缘出现异常高值或低值，先怀疑nodata没有处理干净，查看直方图就能看出来。

4.3 坡度、坡向、填洼的算力与参数取舍

裁剪完成后，最常用的后续产品是坡度和坡向。ArcGIS的Slope工具和QGIS的r.slope.aspect都底层基于Horn算法，该算法用3×3窗口的八个邻域像元拟合局部平面，相比简单差分法对噪声的敏感度更低。用Python里的richdem库或GDAL DEM工具都可以产出同样规格的结果。

GDAL自带dem命令，一行生成坡度：

gdaldem slope dem_jilin_clip.tif slope.tif -p -s 1.0

-p指定输出坡度为百分比，如果不带此参数则输出为度；-s 1.0表示水平和垂直方向比例因子。这套数据已经是30米平面分辨率且单位是米，比例因子保持1.0即可。如果是经纬度数据，这里要按纬度换算比例因子，换算错了坡度会整体偏大或偏小，这地方经常被人忽略。

坡向在分析日照、积温和植被分布时用得多，同样一个gdaldem命令：

gdaldem aspect dem_jilin_clip.tif aspect.tif

坡向输出的角度是从正北方向顺时针测量的，平地将被编码为-9999，分析时记得把这个值单独处理，不能直接参与平均计算。可用圆形统计量（比如把角度分解为sin和cos再求均值）来提取主导坡向，直接对角度做算术平均会导致355度和5度的“假平均”到180度的严重错误。

遇到高起伏山区，下一步通常是水文分析。填洼是水文分析里最容易引起争议的环节。ArcGIS的Fill工具默认会填掉所有洼地，但吉林市这类有喀斯特地貌或采石坑的区域，很多“洼地”是真实地形，填掉就等于抹掉了真实汇水区。常见做法是设置一个填充限制。在arcpy中手工指定z limit可以控制最大填充深度：

arcpy.gp.Fill_sa("dem_jilin_clip.tif", "fill.tif", "", "10")

第三个位置的参数10就是z limit，表示只填充深度在10米以内的洼地，超过10米的保留为真实洼地。深度阈值应该参考项目区域的地貌特征设定，平原区可以给5米，山区给50米甚至100米也不算不合理，关键取决于分析目标是看区域总体汇水还是单个小流域细节。

填洼后，可以输出流向和流量累积栅格，用于提取河网、划定子流域。在QGIS里用r.watershed或者ArcGIS的Flow Direction + Flow Accumulation两条链路即可完成。需要注意，Flow Accumulation计算出的值表示累积像元数量，而非实际径流量。如果要做洪水模拟，需要将像元数乘以30×30=900平方米，再乘以净雨量系数，才能转化为立方米流量。

5. 验证高程精度与几个少有人提的高阶技巧

5.1 用已知高程点检验DEM的基本可靠性

把DEM裁剪完先别急着分析，用一个简单办法验证高程是否靠谱。找3到5个吉林市范围内的已知高程点，可以从当地测绘成果或导航软件里获得，用Python批量提取对应位置的DEM值。

import rasterio points = [(126.55, 43.84), (126.28, 43.95)] # 经度、纬度，按实际数据调整 with rasterio.open("/data/jilin_dem/dem_jilin_clip.tif") as src: for lon, lat in points: row, col = src.index(lon, lat) val = src.read(1, window=((row, row+1), (col, col+1)))[0][0] print(f"经纬度({lon}, {lat}) 处DEM高程: {val}米")

比对结果如果偏差在正负10米内，对30米分辨率数据来说算正常水平；如果偏差超过50米，要么是投影对错了，要么是数据本身质量有问题，再往下做坡度分析没有意义。

5.2 用无数据区占比和直方图判断裁剪质量

裁剪后的DEM是否有大片空洞，直接看统计量。用rasterio读取数组后，计算nodata值在整幅图中的比例，可以快速判断裁出来的成果是否符合要求。

import numpy as np import rasterio with rasterio.open("/data/jilin_dem/dem_jilin_clip.tif") as src: dem = src.read(1).astype(np.float32) nodata = src.nodata valid = dem[dem != nodata] nodata_ratio = (dem.size - valid.size) / dem.size print(f"无数据区占比: {nodata_ratio:.2%}") print(f"有效高程范围: {valid.min():.1f} ~ {valid.max():.1f}米") print(f"平均高程: {valid.mean():.1f}米")

吉林市城区平均海拔大约200到300米，周边山地可达1200米以上。如果打印出来的高程范围完全背离这个数字，比如出现0或负几千，就要考虑原始数据是否已经做过某种预处理，比如把非陆地区域改成了固定值。30米DEM里出现几处nodata是正常的，占比超过5%则要考虑数据源问题。

5.3 一个容易被忽视的技巧：shp转txt提取边界坐标用于定点分析

很多地表位移分析、基站覆盖软件并不直接接受shp，而是要求输入经纬度格式的边界坐标。把吉林市本级的边界shp转成txt，在行业对接、数据库入库和跨平台交换时非常实用。

import geopandas as gpd gdf = gpd.read_file("/data/jilin_dem/吉林市本级_4547.shp") gdf = gdf.to_crs(epsg=4490) # 转成经纬度 with open("/data/jilin_dem/jilin_boundary.txt", "w", encoding="utf-8") as f: for geom in gdf.geometry: if geom.geom_type == "Polygon": coords = geom.exterior.coords elif geom.geom_type == "MultiPolygon": coords = [p.exterior.coords for p in geom.geoms] else: continue if isinstance(coords, list): for ring in coords: for lon, lat in ring: f.write(f"{lon:.6f},{lat:.6f}\n") else: for lon, lat in coords: f.write(f"{lon:.6f},{lat:.6f}\n")

txt生成后，保留6位小数约等于0.1米精度，作为概略边界完全够用。这份txt还能直接当作散点数据导入各类数值模拟软件和Matlab/Python的路径规划模块，省去再写一次解析接口。

另外一个常用技巧是根据30m数据的精细度做渔网分割。地形分析遇到超大区域时，逐像元计算量大、改参数要反复重跑，把整个吉林市按5km×5km或10km×10km切块分开处理是个成熟套路，用ArcGIS的Create Fishnet或Python的geopandas分块裁剪都行。分割后的每块独立做填洼和流向计算，再拼接回整幅，能显著降低单次运算内存占用，也便于分布式并行处理。

最后建议对坡度结果做一次可视化巡视，检查河流两岸是不是出现了带状异常高值，这是DEM噪点常见的暴露形式，特别是山谷两侧的伪陡坎。发现这类情况，先用焦点统计滤波或中值滤波处理原始DEM再重算，不要直接拿带噪点的数据出成果。

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