简介吉林省吉林市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_dem2.2 shp不是单个文件是一组配套文件的合称很多第一次接触矢量数据的同学会对着“.shp”找文件其实一个完整的shapefile至少包含三个基础文件完整形态更多。这份数据既然标注“含本市级范围shp文件”理应带齐下面这些。扩展名作用缺失后果.shp要素几何信息点线面坐标无几何基本报废.shx几何索引用于快速定位要素很多软件拒绝打开.dbf属性表如城市名称、行政区代码属性丢失但几何还在.prj坐标系描述WKT或ESRI格式无法判断投影极易叠加错位.cpg属性表字符编码声明中文属性乱码.sbn / .sbx空间索引非必须部分软件打开慢查看文件清单时如果.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.000000000000000Size是栅格的像元行列数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里查出后替换吉林市常见的情况是4490CGCS2000经纬度转4547或4548CGCS2000高斯-克吕格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, cropTrue, nodata-9999, all_touchedFalse ) with rasterio.open( /data/jilin_dem/dem_jilin_clip.tif, w, driverGTiff, heightout_img.shape[1], widthout_img.shape[2], count1, dtypeout_img.dtype, crssrc.crs, transformout_transform, nodata-9999, ) as dst: dst.write(out_img, 1)all_touchedFalse意味着只有像元中心落在shp范围内才保留边缘会更符合真实边界但会导致边界处少量像元被裁掉all_touchedTrue则把任何与边界相交的像元都保留边缘锯齿感低一些面积会略大于真实值。统计面积时用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×30900平方米再乘以净雨量系数才能转化为立方米流量。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, row1), (col, col1)))[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(epsg4490) # 转成经纬度 with open(/data/jilin_dem/jilin_boundary.txt, w, encodingutf-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再重算不要直接拿带噪点的数据出成果。本文还有配套的精品资源点击获取 吉林市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_dem2.2 shp不是单个文件是一组配套文件的合称很多第一次接触矢量数据的同学会对着“.shp”找文件其实一个完整的shapefile至少包含三个基础文件完整形态更多。这份数据既然标注“含本市级范围shp文件”理应带齐下面这些。扩展名作用缺失后果.shp要素几何信息点线面坐标无几何基本报废.shx几何索引用于快速定位要素很多软件拒绝打开.dbf属性表如城市名称、行政区代码属性丢失但几何还在.prj坐标系描述WKT或ESRI格式无法判断投影极易叠加错位.cpg属性表字符编码声明中文属性乱码.sbn / .sbx空间索引非必须部分软件打开慢查看文件清单时如果.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.000000000000000Size是栅格的像元行列数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里查出后替换吉林市常见的情况是4490CGCS2000经纬度转4547或4548CGCS2000高斯-克吕格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, cropTrue, nodata-9999, all_touchedFalse ) with rasterio.open( /data/jilin_dem/dem_jilin_clip.tif, w, driverGTiff, heightout_img.shape[1], widthout_img.shape[2], count1, dtypeout_img.dtype, crssrc.crs, transformout_transform, nodata-9999, ) as dst: dst.write(out_img, 1)all_touchedFalse意味着只有像元中心落在shp范围内才保留边缘会更符合真实边界但会导致边界处少量像元被裁掉all_touchedTrue则把任何与边界相交的像元都保留边缘锯齿感低一些面积会略大于真实值。统计面积时用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×30900平方米再乘以净雨量系数才能转化为立方米流量。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, row1), (col, col1)))[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(epsg4490) # 转成经纬度 with open(/data/jilin_dem/jilin_boundary.txt, w, encodingutf-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再重算不要直接拿带噪点的数据出成果。本文还有配套的精品资源点击获取