尧图网页设计 
简介这份资源是四川省自贡市30米分辨率的DEM数字高程数据包面向地理信息系统学习者、测绘与城市规划从业者及高校师生可用于地形分析、洪水模拟、地质灾害评估与地图制作等教学实践场景。压缩包共12个文件约14.72MB核心为自贡市dem.tif高程栅格辅以ovr金字塔与tfw坐标文件便于快速浏览定位同时提供自贡市范围shp矢量边界配套shx、dbf、prj、sbn、sbx等索引与投影文件以及多个xml元数据保证数据可直接在ArcGIS等平台加载使用。数据覆盖自贡市行政区域并延伸至部分邻近地区便于考虑边界效应或开展更大范围研究。目前已有421人学习下载适合需要完整、开箱即用的市级DEM数据来练习空间分析与地形建模的读者。1. 自贡市 30m DEM 到手之后先搞清楚它能干什么、不能干什么拿到一份「四川省自贡市 DEM 数字高程数据 30m含本市级范围 shp 文件.zip」的时候很多人第一反应是直接拖进 ArcGIS 看效果。我建议先别急先想清楚一件事30m 分辨率的 DEM 到底能支撑什么级别的分析。自贡市总面积大约 4381 平方公里30m 格网意味着每格代表实地 900 平方米全市大概铺满 540 万个栅格像元。这个量级做市域尺度的坡度坡向、汇流分析、地形起伏度完全够用但你要是想拿它做某条村道级别的切剖面设计精度就不够了。这份资源的核心价值在于「配套」——DEM 栅格加本市级行政边界 shp两样东西一起给。做过 GIS 项目的人都知道拿到 DEM 之后第一件事往往就是按行政边界裁剪如果边界文件要另外去找、坐标系还对不上光对齐这一步就能耗掉半天。这份数据把这两个文件打包在一起省掉的就是这个匹配成本。它适合做市域地形分析、水文初步建模、选址适宜性评价、三维地形底图这类工作也适合教学场景下让学生直接上手练裁剪和重采样。需要提前说清楚一个边界30m 分辨率不等于 1:1 万地形图精度。它来自公开的全球 DEM 数据集常见的是 SRTM、ASTER GDEM 或 ALOS 这类源垂直精度在平坦地区大概 ±5 到 ±10 米山区会更差一些。你要是做工程土方计算这个精度不够但做宏观地形格局判断、坡度分级、流域划分它完全能打。后面几章我会按「先验证数据 → 再裁剪对齐 → 然后做典型分析 → 最后避坑」的顺序把这份资源从头到尾跑一遍。2. 拆包先验证坐标系、NoData 与边界对齐的三步检查2.1 为什么不能直接开始分析很多人拿到 DEM 压缩包解压后看到 tif 文件和 shp 文件直接往 ArcGIS 或 QGIS 里一拖就开始做坡度分析。这么做大概率会在后面某个环节翻车——要么裁剪出来是空白的要么面积算出来差了好几倍要么坡度值明显不对。根因通常出在三个地方坐标系不一致、NoData 值没处理、边界文件和 DEM 的实际覆盖范围对不上。先说坐标系。国内常见的 DEM 数据可能带三种坐标信息WGS84 地理坐标系经纬度单位是度、CGCS2000 地理坐标系、或者投影后的平面坐标系比如 UTM 或高斯-克吕格单位是米。自贡市位于东经 104° 到 105° 之间如果数据是地理坐标系直接算坡度会得到一个荒谬的结果——因为坡度计算需要以米为单位的平面距离用度做单位算出来的坡度值没有物理意义。所以第一步必须确认坐标系必要时做投影转换。再说 NoData。DEM 数据在边界外或者水体区域通常用一个大负数比如 -9999 或 -32768表示无效值。如果你不处理做坡度分析时这些值会被当成真实高程参与计算产生一堆异常值。最后说边界对齐shp 文件的行政边界和 DEM 的实际覆盖范围如果不完全重叠裁剪时就会出现边缘缺失或者多余区域。2.2 用 Python 做三步验证我一般会用 Python 的 rasterio 和 geopandas 快速做一遍体检比在 GUI 里点来点去更可控。下面这段代码依次检查 DEM 的坐标系、NoData 值、范围和边界文件的范围。import rasterio import geopandas as gpd import numpy as np # 打开 DEM 文件 dem_path zigong_dem_30m.tif with rasterio.open(dem_path) as src: print( DEM 基本信息 ) print(f坐标系: {src.crs}) print(f影像尺寸: {src.width} x {src.height}) print(f波段数: {src.count}) print(fNoData 值: {src.nodata}) print(f地理范围: {src.bounds}) print(f分辨率: {src.res}) # 读取高程数据做统计 elev src.read(1) valid elev[elev ! src.nodata] if src.nodata else elev print(f高程范围: {valid.min():.1f} ~ {valid.max():.1f} 米) print(f高程均值: {valid.mean():.1f} 米) print(f无效值占比: {(1 - valid.size / elev.size) * 100:.2f}%) # 打开行政边界 shp shp_path zigong_boundary.shp gdf gpd.read_file(shp_path) print(\n 边界文件信息 ) print(f坐标系: {gdf.crs}) print(f要素数量: {len(gdf)}) print(f边界范围: {gdf.total_bounds}) print(f面积(原始单位): {gdf.geometry.area.sum():.4f})这段代码的逻辑很直接先看 DEM 的元数据确认坐标系和 NoData 值再读高程数组看有效值的分布范围是否合理自贡市海拔大概在 250 到 900 米之间如果算出来是负数或者上万说明 NoData 没排除干净最后看边界文件的坐标系和范围和 DEM 的 bounds 做对比。参数方面有几个关键点。src.nodata如果返回 None说明这个文件没有显式声明 NoData 值你需要手动检查数据中的异常值。src.res返回的是像元分辨率如果是地理坐标系单位是度30m 约等于 0.000277 度如果是投影坐标系单位是米。gdf.total_bounds返回的是边界文件的四至范围格式是 (minx, miny, maxx, maxy)你可以直接和 DEM 的 bounds 做对比看边界是否完全落在 DEM 覆盖范围内。2.3 坐标系不一致时的处理如果检查发现 DEM 是地理坐标系EPSG:4326而边界文件是投影坐标系比如 EPSG:32648 或 CGCS2000 高斯投影你需要统一到同一个坐标系再做裁剪。常见做法是统一到投影坐标系因为后续做坡度、汇流这些分析都需要以米为单位。# 假设 DEM 是 EPSG:4326边界是 EPSG:4544CGCS2000 3度带 # 统一将边界转到 DEM 的坐标系或者反过来 # 推荐都转到投影坐标系 EPSG:4544 from rasterio.warp import calculate_default_transform, reproject, Resampling dst_crs EPSG:4544 with rasterio.open(dem_path) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds ) kwargs src.meta.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height }) with rasterio.open(zigong_dem_30m_projected.tif, w, **kwargs) as dst: for i in range(1, src.count 1): reproject( sourcerasterio.band(src, i), destinationrasterio.band(dst, i), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.bilinear # 高程数据用双线性插值 )这里有个细节值得注意重投影时重采样方法的选择。高程数据推荐用双线性插值bilinear因为它能保持高程的连续变化趋势不要用最近邻nearest那会产生阶梯状的高程突变。但如果是分类数据比如土地利用类型就必须用最近邻否则会出现不存在的类别值。这个区别在做数据融合时特别容易搞混。3. 按行政边界裁剪 DEM掩膜提取与面积校验3.1 裁剪的两种思路按行政边界裁剪 DEM本质上是用矢量多边形做掩膜提取。常见做法有两种一种是用 ArcGIS 的「按掩膜提取」工具另一种是用 Python 的 rasterio.mask 模块。两种方法结果一样但 Python 方式更容易批量化和版本控制。裁剪之前有一个容易被忽略的问题边界文件可能包含多个多边形要素比如自贡市下辖的自流井区、贡井区、大安区、沿滩区、荣县、富顺县等你需要确认是想要整个市域合并后的边界还是每个区县单独裁剪。如果是前者先用 dissolve 合并如果是后者就按属性字段逐要素裁剪。3.2 用 rasterio.mask 做裁剪import rasterio from rasterio.mask import mask import geopandas as gpd import json # 读取边界并合并为单个多边形 gdf gpd.read_file(zigong_boundary.shp) gdf_union gdf.dissolve() # 合并所有区县为市域整体 # 确保坐标系一致 dem_path zigong_dem_30m_projected.tif with rasterio.open(dem_path) as src: if gdf_union.crs ! src.crs: gdf_union gdf_union.to_crs(src.crs) # 提取几何体 geoms [json.loads(gdf_union.to_json())[features][0][geometry]] # 执行掩膜裁剪 out_image, out_transform mask(src, geoms, cropTrue, nodata-9999) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: -9999 }) # 写出裁剪结果 with rasterio.open(zigong_dem_clipped.tif, w, **out_meta) as dst: dst.write(out_image) print(f裁剪后尺寸: {out_image.shape[1]} x {out_image.shape[2]}) print(f裁剪后范围: {out_transform})这段代码的关键参数有三个。cropTrue表示裁剪后自动缩小影像范围到边界的最小外接矩形去掉多余空白如果设为 False输出影像尺寸和原始 DEM 一样大只是边界外被设为 NoData。nodata-9999指定裁剪后边界外的填充值这个值要和后续分析工具兼容。gdf_union.to_json()把 GeoDataFrame 转成 GeoJSON 格式再提取 geometry 字段这是 rasterio.mask 要求的输入格式。3.3 裁剪后的面积校验裁剪完成后必须做一次面积校验。方法很简单用裁剪后的 DEM 有效像元数乘以单个像元面积和边界文件的实际面积做对比。如果差异超过 5%说明裁剪过程中有问题。import rasterio import numpy as np import geopandas as gpd # 计算 DEM 有效像元面积 with rasterio.open(zigong_dem_clipped.tif) as src: data src.read(1) valid_mask data ! src.nodata pixel_area src.res[0] * src.res[1] # 像元面积平方米 dem_area valid_mask.sum() * pixel_area / 1e6 # 转为平方公里 print(fDEM 有效像元面积: {dem_area:.2f} 平方公里) # 计算边界文件面积 gdf gpd.read_file(zigong_boundary.shp) if gdf.crs and gdf.crs.is_geographic: gdf gdf.to_crs(EPSG:4544) # 转投影坐标系后才能算面积 shp_area gdf.geometry.area.sum() / 1e6 print(f边界文件面积: {shp_area:.2f} 平方公里) print(f面积差异: {abs(dem_area - shp_area) / shp_area * 100:.2f}%)自贡市的官方面积大约 4381 平方公里如果你算出来在这个数值附近±5% 以内说明裁剪和坐标系处理都是对的。如果差得远优先检查坐标系是否统一、边界文件是否包含了所有区县。注意地理坐标系下直接算面积得到的是平方度没有实际意义。必须转到投影坐标系后再计算面积这是新手最容易踩的坑之一。4. 从 DEM 到地形因子坡度、坡向与汇流分析的参数设置4.1 坡度与坡向计算裁剪完成、面积校验通过之后就可以做地形因子提取了。最基础的是坡度和坡向。坡度反映地表倾斜程度坡向反映坡面朝向这两个因子在选址、农业、水文分析中都是基础输入。用 Python 的 richdem 或 GDAL 都可以算我这里用 GDAL 的命令行方式因为它稳定且参数清晰。# 计算坡度单位度 gdaldem slope zigong_dem_clipped.tif zigong_slope.tif -of GTiff -b 1 -s 1.0 # 计算坡向 gdaldem aspect zigong_dem_clipped.tif zigong_aspect.tif -of GTiff -b 1 # 计算地形起伏度用焦点统计窗口 3x3 gdaldem roughness zigong_dem_clipped.tif zigong_roughness.tif -of GTiffgdaldem slope的-s 1.0参数是比例因子。当 DEM 的 XY 单位是米、Z 单位也是米时比例因子设为 1.0如果 Z 单位是英尺而 XY 是米就需要设成 0.3048。这个参数设错了坡度值会差好几倍。-b 1指定用第一个波段因为裁剪后的 DEM 只有一个波段。坡向输出的值范围是 0 到 360 度0 表示正北90 表示正东180 表示正南270 表示正西。平坦区域坡度接近 0的坡向值通常设为 -1 或 9999具体取决于工具实现后续分析时要记得排除。4.2 汇流累积量计算汇流分析是水文建模的核心步骤目的是算出每个像元上游有多少个像元的水会流经它。这个值越大说明该位置越可能是河道或汇水区。完整流程是填洼 → 流向 → 汇流累积。# 第一步填洼消除 DEM 中的伪洼地 gdal_fillnodata.py -md 10 zigong_dem_clipped.tif zigong_dem_filled.tif # 第二步计算流向D8 算法 gdaldem TRI zigong_dem_filled.tif zigong_tri.tif # 注意GDAL 没有直接的流向命令实际常用 richdem 或 TauDEM这里要说明一下GDAL 本身没有内置 D8 流向和汇流累积的命令实际工作中我一般用 richdem 这个 Python 库来做它的 API 简洁且算法实现可靠。import richdem as rd import numpy as np # 读取填洼后的 DEM dem rd.LoadGDAL(zigong_dem_filled.tif) # 填洼如果上一步没做 dem_filled rd.FillDepressions(dem, epsilonTrue, in_placeFalse) # 计算 D8 流向 flow_dir rd.FlowDirections(dem_filled, methodD8) # 计算汇流累积量 accum rd.FlowAccumulation(flow_dir, methodD8) # 保存结果 rd.SaveGDAL(zigong_flow_accum.tif, accum) print(f汇流累积最大值: {np.max(accum)}) print(f汇流累积均值: {np.mean(accum):.1f})FillDepressions的epsilonTrue参数表示在填洼时给平坦区域一个微小梯度避免后续流向计算出现死循环。FlowDirections的methodD8指定用八方向算法这是最经典的流向算法每个像元的水流方向指向周围八个邻居中坡度最陡的那个。FlowAccumulation输出的值表示每个像元上游汇水像元的总数最大值通常出现在流域出口位置。4.3 参数选择与结果解读坡度计算中比例因子是最关键的参数。汇流分析中填洼的阈值和流向算法是核心。D8 算法简单快速但在平坦区域会产生平行流线这是它的固有缺陷。如果你对水文精度要求高可以考虑 D-infinity 或 MFD 算法但计算量会大很多。汇流累积量的解读需要结合阈值。通常认为汇流累积量超过某个阈值的像元就是河道。这个阈值没有固定标准取决于流域面积和 DEM 分辨率。对于 30m 分辨率的自贡市数据我一般先用 1000 作为初始阈值试一下然后根据提取出的河网密度做调整。阈值太小河网过密到处都是沟阈值太大只留下主干河道支流全丢了。提示填洼步骤会改变原始高程值如果你后续需要报告实际高程记得保留一份未填洼的原始数据。填洼只用于水文分析不要用它替代原始 DEM 做其他分析。5. 避坑与排查30m DEM 处理中最容易翻车的五个地方5.1 裁剪结果全空白或只有边缘一圈现象用边界文件裁剪 DEM 后打开结果发现大部分区域是 NoData只有边缘少量像元有值。原因最常见的情况是边界文件和 DEM 的坐标系不一致但两者看起来「差不多」导致裁剪时几何体位置偏移。另一种可能是边界文件的多边形是「洞」而不是「实体」或者几何体本身无效自相交、未闭合。解决先用gdf.is_valid检查几何有效性无效的用gdf.buffer(0)修复。然后强制统一坐标系不要凭肉眼判断。最后用gdf.total_bounds和src.bounds做数值对比确认范围确实重叠。5.2 坡度值异常偏大或偏小现象算出来的坡度最大值超过 80 度或者全市坡度都在 5 度以下明显不符合实际地形。原因比例因子设错了。如果 DEM 是地理坐标系度而你没有做投影转换就直接算坡度GDAL 会默认用度作为 XY 单位算出来的坡度值完全错误。另一种可能是 Z 单位不是米比如是英尺但比例因子没调整。解决算坡度之前确认 DEM 已经转到投影坐标系XY 和 Z 单位都是米。用gdalsrsinfo查看坐标系信息用gdalinfo -stats查看高程统计值确认数值范围合理。5.3 汇流分析出现大量平行流线现象汇流累积结果中平坦区域出现大量平行的直线流线明显不符合自然水系形态。原因这是 D8 算法的固有缺陷。在平坦区域D8 算法无法确定唯一的最陡方向导致水流沿着固定方向平行流动。解决填洼时使用epsilonTrue给平坦区域一个微小梯度。如果问题仍然严重考虑改用 D-infinity 算法它会把水流按角度分配到多个下游像元能有效消除平行流线。代价是计算量增加且结果不再是整数。5.4 边界文件属性表混乱导致裁剪错误现象裁剪时报错「几何体无效」或裁剪结果只包含部分区县。原因边界文件的属性表可能包含多个图层或重复要素dissolve()合并时没有指定正确的分组字段。另外有些 shp 文件的几何类型是 MultiPolygon直接提取 geometry 时格式可能不对。解决先print(gdf.head())看属性表结构确认要素数量和几何类型。用gdf.dissolve(bycity_name)按字段合并而不是无参数合并。提取几何体时用gdf.geometry.unary_union而不是手动转 JSON。5.5 NoData 值被当成真实高程参与统计现象高程统计结果中出现 -9999 或 -32768 这样的异常值导致均值、最大值完全失真。原因读取数据时没有排除 NoData 像元。rasterio 读取的数组是原始值NoData 不会被自动过滤。解决读取后立即用data[data ! src.nodata]过滤或者用np.ma.masked_equal(data, src.nodata)创建掩膜数组。如果src.nodata返回 None需要手动检查数据中的异常值常见的有 -9999、-32768、0 等。6. 进阶技巧用 GDAL 命令行批量处理与结果验证6.1 批量裁剪多个区县实际项目中经常需要按区县分别裁剪 DEM而不是只裁一个市域整体。手动一个个裁太慢用 GDAL 命令行配合 shell 脚本可以批量完成。#!/bin/bash # 按区县批量裁剪 DEM DEMzigong_dem_30m_projected.tif SHPzigong_boundary.shp OUTPUT_DIRcounty_dem mkdir -p $OUTPUT_DIR # 获取所有区县名称 COUNTIES$(ogrinfo -al -so $SHP | grep 区\|县 | awk {print $NF}) for county in $COUNTIES; do echo 正在裁剪: $county # 用 ogr2ogr 提取单个区县边界 ogr2ogr -where NAME$county ${OUTPUT_DIR}/${county}.shp $SHP # 用 gdalwarp 按边界裁剪 gdalwarp -cutline ${OUTPUT_DIR}/${county}.shp \ -crop_to_cutline \ -dstnodata -9999 \ $DEM ${OUTPUT_DIR}/${county}_dem.tif done echo 批量裁剪完成这个脚本的核心是gdalwarp的-cutline和-crop_to_cutline参数。-cutline指定裁剪用的矢量边界-crop_to_cutline表示裁剪后自动缩小影像范围。-dstnodata -9999指定边界外的填充值。相比 Python 的 rasterio.maskgdalwarp 在处理大量文件时更稳定而且支持并行。6.2 用 gdalinfo 做结果验证每个裁剪结果出来后用gdalinfo快速检查一下基本信息。# 查看裁剪结果的统计信息 gdalinfo -stats county_dem/自流井区_dem.tif | grep -E Size|Coordinate|Minimum|Maximum|Mean|NoData # 对比原始 DEM 和裁剪后的范围 gdalinfo -stats zigong_dem_30m_projected.tif | grep Upper Left\|Lower Right-stats参数会强制计算统计值包括最小值、最大值、均值、标准差。你可以快速判断高程范围是否合理。如果某个区县的 DEM 最大值突然变成 9999 或 -9999说明 NoData 处理有问题。6.3 一个我踩过的坑有一次我拿到一份 DEM 数据坐标系显示是 CGCS2000但实际数据是 WGS84 的坐标值只是元数据标错了。结果裁剪出来的边界和 DEM 偏移了大概 50 米肉眼看不出来但做汇流分析时河道位置全错了。从那以后我每次拿到新数据都强制走一遍「坐标系验证」流程用已知地物的坐标做交叉验证比如取一个明显的山峰或河流交汇点在 DEM 和边界文件里分别查坐标对比是否一致。具体做法是在 DEM 上找一个最高点记下它的行列号和坐标然后在 Google Earth 或天地图上找到同一个位置对比坐标差异。如果差异在几十米以内说明坐标系没问题如果差了几百米甚至更多说明坐标系标错了需要手动校正。这个验证步骤花不了五分钟但能避免后面几天的返工。希望帮到你。本文还有配套的精品资源点击获取 自贡市30m DEM数据处理全流程：从验证、裁剪到地形分析避坑指南
简介这份资源是四川省自贡市30米分辨率的DEM数字高程数据包面向地理信息系统学习者、测绘与城市规划从业者及高校师生可用于地形分析、洪水模拟、地质灾害评估与地图制作等教学实践场景。压缩包共12个文件约14.72MB核心为自贡市dem.tif高程栅格辅以ovr金字塔与tfw坐标文件便于快速浏览定位同时提供自贡市范围shp矢量边界配套shx、dbf、prj、sbn、sbx等索引与投影文件以及多个xml元数据保证数据可直接在ArcGIS等平台加载使用。数据覆盖自贡市行政区域并延伸至部分邻近地区便于考虑边界效应或开展更大范围研究。目前已有421人学习下载适合需要完整、开箱即用的市级DEM数据来练习空间分析与地形建模的读者。1. 自贡市 30m DEM 到手之后先搞清楚它能干什么、不能干什么拿到一份「四川省自贡市 DEM 数字高程数据 30m含本市级范围 shp 文件.zip」的时候很多人第一反应是直接拖进 ArcGIS 看效果。我建议先别急先想清楚一件事30m 分辨率的 DEM 到底能支撑什么级别的分析。自贡市总面积大约 4381 平方公里30m 格网意味着每格代表实地 900 平方米全市大概铺满 540 万个栅格像元。这个量级做市域尺度的坡度坡向、汇流分析、地形起伏度完全够用但你要是想拿它做某条村道级别的切剖面设计精度就不够了。这份资源的核心价值在于「配套」——DEM 栅格加本市级行政边界 shp两样东西一起给。做过 GIS 项目的人都知道拿到 DEM 之后第一件事往往就是按行政边界裁剪如果边界文件要另外去找、坐标系还对不上光对齐这一步就能耗掉半天。这份数据把这两个文件打包在一起省掉的就是这个匹配成本。它适合做市域地形分析、水文初步建模、选址适宜性评价、三维地形底图这类工作也适合教学场景下让学生直接上手练裁剪和重采样。需要提前说清楚一个边界30m 分辨率不等于 1:1 万地形图精度。它来自公开的全球 DEM 数据集常见的是 SRTM、ASTER GDEM 或 ALOS 这类源垂直精度在平坦地区大概 ±5 到 ±10 米山区会更差一些。你要是做工程土方计算这个精度不够但做宏观地形格局判断、坡度分级、流域划分它完全能打。后面几章我会按「先验证数据 → 再裁剪对齐 → 然后做典型分析 → 最后避坑」的顺序把这份资源从头到尾跑一遍。2. 拆包先验证坐标系、NoData 与边界对齐的三步检查2.1 为什么不能直接开始分析很多人拿到 DEM 压缩包解压后看到 tif 文件和 shp 文件直接往 ArcGIS 或 QGIS 里一拖就开始做坡度分析。这么做大概率会在后面某个环节翻车——要么裁剪出来是空白的要么面积算出来差了好几倍要么坡度值明显不对。根因通常出在三个地方坐标系不一致、NoData 值没处理、边界文件和 DEM 的实际覆盖范围对不上。先说坐标系。国内常见的 DEM 数据可能带三种坐标信息WGS84 地理坐标系经纬度单位是度、CGCS2000 地理坐标系、或者投影后的平面坐标系比如 UTM 或高斯-克吕格单位是米。自贡市位于东经 104° 到 105° 之间如果数据是地理坐标系直接算坡度会得到一个荒谬的结果——因为坡度计算需要以米为单位的平面距离用度做单位算出来的坡度值没有物理意义。所以第一步必须确认坐标系必要时做投影转换。再说 NoData。DEM 数据在边界外或者水体区域通常用一个大负数比如 -9999 或 -32768表示无效值。如果你不处理做坡度分析时这些值会被当成真实高程参与计算产生一堆异常值。最后说边界对齐shp 文件的行政边界和 DEM 的实际覆盖范围如果不完全重叠裁剪时就会出现边缘缺失或者多余区域。2.2 用 Python 做三步验证我一般会用 Python 的 rasterio 和 geopandas 快速做一遍体检比在 GUI 里点来点去更可控。下面这段代码依次检查 DEM 的坐标系、NoData 值、范围和边界文件的范围。import rasterio import geopandas as gpd import numpy as np # 打开 DEM 文件 dem_path zigong_dem_30m.tif with rasterio.open(dem_path) as src: print( DEM 基本信息 ) print(f坐标系: {src.crs}) print(f影像尺寸: {src.width} x {src.height}) print(f波段数: {src.count}) print(fNoData 值: {src.nodata}) print(f地理范围: {src.bounds}) print(f分辨率: {src.res}) # 读取高程数据做统计 elev src.read(1) valid elev[elev ! src.nodata] if src.nodata else elev print(f高程范围: {valid.min():.1f} ~ {valid.max():.1f} 米) print(f高程均值: {valid.mean():.1f} 米) print(f无效值占比: {(1 - valid.size / elev.size) * 100:.2f}%) # 打开行政边界 shp shp_path zigong_boundary.shp gdf gpd.read_file(shp_path) print(\n 边界文件信息 ) print(f坐标系: {gdf.crs}) print(f要素数量: {len(gdf)}) print(f边界范围: {gdf.total_bounds}) print(f面积(原始单位): {gdf.geometry.area.sum():.4f})这段代码的逻辑很直接先看 DEM 的元数据确认坐标系和 NoData 值再读高程数组看有效值的分布范围是否合理自贡市海拔大概在 250 到 900 米之间如果算出来是负数或者上万说明 NoData 没排除干净最后看边界文件的坐标系和范围和 DEM 的 bounds 做对比。参数方面有几个关键点。src.nodata如果返回 None说明这个文件没有显式声明 NoData 值你需要手动检查数据中的异常值。src.res返回的是像元分辨率如果是地理坐标系单位是度30m 约等于 0.000277 度如果是投影坐标系单位是米。gdf.total_bounds返回的是边界文件的四至范围格式是 (minx, miny, maxx, maxy)你可以直接和 DEM 的 bounds 做对比看边界是否完全落在 DEM 覆盖范围内。2.3 坐标系不一致时的处理如果检查发现 DEM 是地理坐标系EPSG:4326而边界文件是投影坐标系比如 EPSG:32648 或 CGCS2000 高斯投影你需要统一到同一个坐标系再做裁剪。常见做法是统一到投影坐标系因为后续做坡度、汇流这些分析都需要以米为单位。# 假设 DEM 是 EPSG:4326边界是 EPSG:4544CGCS2000 3度带 # 统一将边界转到 DEM 的坐标系或者反过来 # 推荐都转到投影坐标系 EPSG:4544 from rasterio.warp import calculate_default_transform, reproject, Resampling dst_crs EPSG:4544 with rasterio.open(dem_path) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds ) kwargs src.meta.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height }) with rasterio.open(zigong_dem_30m_projected.tif, w, **kwargs) as dst: for i in range(1, src.count 1): reproject( sourcerasterio.band(src, i), destinationrasterio.band(dst, i), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.bilinear # 高程数据用双线性插值 )这里有个细节值得注意重投影时重采样方法的选择。高程数据推荐用双线性插值bilinear因为它能保持高程的连续变化趋势不要用最近邻nearest那会产生阶梯状的高程突变。但如果是分类数据比如土地利用类型就必须用最近邻否则会出现不存在的类别值。这个区别在做数据融合时特别容易搞混。3. 按行政边界裁剪 DEM掩膜提取与面积校验3.1 裁剪的两种思路按行政边界裁剪 DEM本质上是用矢量多边形做掩膜提取。常见做法有两种一种是用 ArcGIS 的「按掩膜提取」工具另一种是用 Python 的 rasterio.mask 模块。两种方法结果一样但 Python 方式更容易批量化和版本控制。裁剪之前有一个容易被忽略的问题边界文件可能包含多个多边形要素比如自贡市下辖的自流井区、贡井区、大安区、沿滩区、荣县、富顺县等你需要确认是想要整个市域合并后的边界还是每个区县单独裁剪。如果是前者先用 dissolve 合并如果是后者就按属性字段逐要素裁剪。3.2 用 rasterio.mask 做裁剪import rasterio from rasterio.mask import mask import geopandas as gpd import json # 读取边界并合并为单个多边形 gdf gpd.read_file(zigong_boundary.shp) gdf_union gdf.dissolve() # 合并所有区县为市域整体 # 确保坐标系一致 dem_path zigong_dem_30m_projected.tif with rasterio.open(dem_path) as src: if gdf_union.crs ! src.crs: gdf_union gdf_union.to_crs(src.crs) # 提取几何体 geoms [json.loads(gdf_union.to_json())[features][0][geometry]] # 执行掩膜裁剪 out_image, out_transform mask(src, geoms, cropTrue, nodata-9999) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: -9999 }) # 写出裁剪结果 with rasterio.open(zigong_dem_clipped.tif, w, **out_meta) as dst: dst.write(out_image) print(f裁剪后尺寸: {out_image.shape[1]} x {out_image.shape[2]}) print(f裁剪后范围: {out_transform})这段代码的关键参数有三个。cropTrue表示裁剪后自动缩小影像范围到边界的最小外接矩形去掉多余空白如果设为 False输出影像尺寸和原始 DEM 一样大只是边界外被设为 NoData。nodata-9999指定裁剪后边界外的填充值这个值要和后续分析工具兼容。gdf_union.to_json()把 GeoDataFrame 转成 GeoJSON 格式再提取 geometry 字段这是 rasterio.mask 要求的输入格式。3.3 裁剪后的面积校验裁剪完成后必须做一次面积校验。方法很简单用裁剪后的 DEM 有效像元数乘以单个像元面积和边界文件的实际面积做对比。如果差异超过 5%说明裁剪过程中有问题。import rasterio import numpy as np import geopandas as gpd # 计算 DEM 有效像元面积 with rasterio.open(zigong_dem_clipped.tif) as src: data src.read(1) valid_mask data ! src.nodata pixel_area src.res[0] * src.res[1] # 像元面积平方米 dem_area valid_mask.sum() * pixel_area / 1e6 # 转为平方公里 print(fDEM 有效像元面积: {dem_area:.2f} 平方公里) # 计算边界文件面积 gdf gpd.read_file(zigong_boundary.shp) if gdf.crs and gdf.crs.is_geographic: gdf gdf.to_crs(EPSG:4544) # 转投影坐标系后才能算面积 shp_area gdf.geometry.area.sum() / 1e6 print(f边界文件面积: {shp_area:.2f} 平方公里) print(f面积差异: {abs(dem_area - shp_area) / shp_area * 100:.2f}%)自贡市的官方面积大约 4381 平方公里如果你算出来在这个数值附近±5% 以内说明裁剪和坐标系处理都是对的。如果差得远优先检查坐标系是否统一、边界文件是否包含了所有区县。注意地理坐标系下直接算面积得到的是平方度没有实际意义。必须转到投影坐标系后再计算面积这是新手最容易踩的坑之一。4. 从 DEM 到地形因子坡度、坡向与汇流分析的参数设置4.1 坡度与坡向计算裁剪完成、面积校验通过之后就可以做地形因子提取了。最基础的是坡度和坡向。坡度反映地表倾斜程度坡向反映坡面朝向这两个因子在选址、农业、水文分析中都是基础输入。用 Python 的 richdem 或 GDAL 都可以算我这里用 GDAL 的命令行方式因为它稳定且参数清晰。# 计算坡度单位度 gdaldem slope zigong_dem_clipped.tif zigong_slope.tif -of GTiff -b 1 -s 1.0 # 计算坡向 gdaldem aspect zigong_dem_clipped.tif zigong_aspect.tif -of GTiff -b 1 # 计算地形起伏度用焦点统计窗口 3x3 gdaldem roughness zigong_dem_clipped.tif zigong_roughness.tif -of GTiffgdaldem slope的-s 1.0参数是比例因子。当 DEM 的 XY 单位是米、Z 单位也是米时比例因子设为 1.0如果 Z 单位是英尺而 XY 是米就需要设成 0.3048。这个参数设错了坡度值会差好几倍。-b 1指定用第一个波段因为裁剪后的 DEM 只有一个波段。坡向输出的值范围是 0 到 360 度0 表示正北90 表示正东180 表示正南270 表示正西。平坦区域坡度接近 0的坡向值通常设为 -1 或 9999具体取决于工具实现后续分析时要记得排除。4.2 汇流累积量计算汇流分析是水文建模的核心步骤目的是算出每个像元上游有多少个像元的水会流经它。这个值越大说明该位置越可能是河道或汇水区。完整流程是填洼 → 流向 → 汇流累积。# 第一步填洼消除 DEM 中的伪洼地 gdal_fillnodata.py -md 10 zigong_dem_clipped.tif zigong_dem_filled.tif # 第二步计算流向D8 算法 gdaldem TRI zigong_dem_filled.tif zigong_tri.tif # 注意GDAL 没有直接的流向命令实际常用 richdem 或 TauDEM这里要说明一下GDAL 本身没有内置 D8 流向和汇流累积的命令实际工作中我一般用 richdem 这个 Python 库来做它的 API 简洁且算法实现可靠。import richdem as rd import numpy as np # 读取填洼后的 DEM dem rd.LoadGDAL(zigong_dem_filled.tif) # 填洼如果上一步没做 dem_filled rd.FillDepressions(dem, epsilonTrue, in_placeFalse) # 计算 D8 流向 flow_dir rd.FlowDirections(dem_filled, methodD8) # 计算汇流累积量 accum rd.FlowAccumulation(flow_dir, methodD8) # 保存结果 rd.SaveGDAL(zigong_flow_accum.tif, accum) print(f汇流累积最大值: {np.max(accum)}) print(f汇流累积均值: {np.mean(accum):.1f})FillDepressions的epsilonTrue参数表示在填洼时给平坦区域一个微小梯度避免后续流向计算出现死循环。FlowDirections的methodD8指定用八方向算法这是最经典的流向算法每个像元的水流方向指向周围八个邻居中坡度最陡的那个。FlowAccumulation输出的值表示每个像元上游汇水像元的总数最大值通常出现在流域出口位置。4.3 参数选择与结果解读坡度计算中比例因子是最关键的参数。汇流分析中填洼的阈值和流向算法是核心。D8 算法简单快速但在平坦区域会产生平行流线这是它的固有缺陷。如果你对水文精度要求高可以考虑 D-infinity 或 MFD 算法但计算量会大很多。汇流累积量的解读需要结合阈值。通常认为汇流累积量超过某个阈值的像元就是河道。这个阈值没有固定标准取决于流域面积和 DEM 分辨率。对于 30m 分辨率的自贡市数据我一般先用 1000 作为初始阈值试一下然后根据提取出的河网密度做调整。阈值太小河网过密到处都是沟阈值太大只留下主干河道支流全丢了。提示填洼步骤会改变原始高程值如果你后续需要报告实际高程记得保留一份未填洼的原始数据。填洼只用于水文分析不要用它替代原始 DEM 做其他分析。5. 避坑与排查30m DEM 处理中最容易翻车的五个地方5.1 裁剪结果全空白或只有边缘一圈现象用边界文件裁剪 DEM 后打开结果发现大部分区域是 NoData只有边缘少量像元有值。原因最常见的情况是边界文件和 DEM 的坐标系不一致但两者看起来「差不多」导致裁剪时几何体位置偏移。另一种可能是边界文件的多边形是「洞」而不是「实体」或者几何体本身无效自相交、未闭合。解决先用gdf.is_valid检查几何有效性无效的用gdf.buffer(0)修复。然后强制统一坐标系不要凭肉眼判断。最后用gdf.total_bounds和src.bounds做数值对比确认范围确实重叠。5.2 坡度值异常偏大或偏小现象算出来的坡度最大值超过 80 度或者全市坡度都在 5 度以下明显不符合实际地形。原因比例因子设错了。如果 DEM 是地理坐标系度而你没有做投影转换就直接算坡度GDAL 会默认用度作为 XY 单位算出来的坡度值完全错误。另一种可能是 Z 单位不是米比如是英尺但比例因子没调整。解决算坡度之前确认 DEM 已经转到投影坐标系XY 和 Z 单位都是米。用gdalsrsinfo查看坐标系信息用gdalinfo -stats查看高程统计值确认数值范围合理。5.3 汇流分析出现大量平行流线现象汇流累积结果中平坦区域出现大量平行的直线流线明显不符合自然水系形态。原因这是 D8 算法的固有缺陷。在平坦区域D8 算法无法确定唯一的最陡方向导致水流沿着固定方向平行流动。解决填洼时使用epsilonTrue给平坦区域一个微小梯度。如果问题仍然严重考虑改用 D-infinity 算法它会把水流按角度分配到多个下游像元能有效消除平行流线。代价是计算量增加且结果不再是整数。5.4 边界文件属性表混乱导致裁剪错误现象裁剪时报错「几何体无效」或裁剪结果只包含部分区县。原因边界文件的属性表可能包含多个图层或重复要素dissolve()合并时没有指定正确的分组字段。另外有些 shp 文件的几何类型是 MultiPolygon直接提取 geometry 时格式可能不对。解决先print(gdf.head())看属性表结构确认要素数量和几何类型。用gdf.dissolve(bycity_name)按字段合并而不是无参数合并。提取几何体时用gdf.geometry.unary_union而不是手动转 JSON。5.5 NoData 值被当成真实高程参与统计现象高程统计结果中出现 -9999 或 -32768 这样的异常值导致均值、最大值完全失真。原因读取数据时没有排除 NoData 像元。rasterio 读取的数组是原始值NoData 不会被自动过滤。解决读取后立即用data[data ! src.nodata]过滤或者用np.ma.masked_equal(data, src.nodata)创建掩膜数组。如果src.nodata返回 None需要手动检查数据中的异常值常见的有 -9999、-32768、0 等。6. 进阶技巧用 GDAL 命令行批量处理与结果验证6.1 批量裁剪多个区县实际项目中经常需要按区县分别裁剪 DEM而不是只裁一个市域整体。手动一个个裁太慢用 GDAL 命令行配合 shell 脚本可以批量完成。#!/bin/bash # 按区县批量裁剪 DEM DEMzigong_dem_30m_projected.tif SHPzigong_boundary.shp OUTPUT_DIRcounty_dem mkdir -p $OUTPUT_DIR # 获取所有区县名称 COUNTIES$(ogrinfo -al -so $SHP | grep 区\|县 | awk {print $NF}) for county in $COUNTIES; do echo 正在裁剪: $county # 用 ogr2ogr 提取单个区县边界 ogr2ogr -where NAME$county ${OUTPUT_DIR}/${county}.shp $SHP # 用 gdalwarp 按边界裁剪 gdalwarp -cutline ${OUTPUT_DIR}/${county}.shp \ -crop_to_cutline \ -dstnodata -9999 \ $DEM ${OUTPUT_DIR}/${county}_dem.tif done echo 批量裁剪完成这个脚本的核心是gdalwarp的-cutline和-crop_to_cutline参数。-cutline指定裁剪用的矢量边界-crop_to_cutline表示裁剪后自动缩小影像范围。-dstnodata -9999指定边界外的填充值。相比 Python 的 rasterio.maskgdalwarp 在处理大量文件时更稳定而且支持并行。6.2 用 gdalinfo 做结果验证每个裁剪结果出来后用gdalinfo快速检查一下基本信息。# 查看裁剪结果的统计信息 gdalinfo -stats county_dem/自流井区_dem.tif | grep -E Size|Coordinate|Minimum|Maximum|Mean|NoData # 对比原始 DEM 和裁剪后的范围 gdalinfo -stats zigong_dem_30m_projected.tif | grep Upper Left\|Lower Right-stats参数会强制计算统计值包括最小值、最大值、均值、标准差。你可以快速判断高程范围是否合理。如果某个区县的 DEM 最大值突然变成 9999 或 -9999说明 NoData 处理有问题。6.3 一个我踩过的坑有一次我拿到一份 DEM 数据坐标系显示是 CGCS2000但实际数据是 WGS84 的坐标值只是元数据标错了。结果裁剪出来的边界和 DEM 偏移了大概 50 米肉眼看不出来但做汇流分析时河道位置全错了。从那以后我每次拿到新数据都强制走一遍「坐标系验证」流程用已知地物的坐标做交叉验证比如取一个明显的山峰或河流交汇点在 DEM 和边界文件里分别查坐标对比是否一致。具体做法是在 DEM 上找一个最高点记下它的行列号和坐标然后在 Google Earth 或天地图上找到同一个位置对比坐标差异。如果差异在几十米以内说明坐标系没问题如果差了几百米甚至更多说明坐标系标错了需要手动校正。这个验证步骤花不了五分钟但能避免后面几天的返工。希望帮到你。本文还有配套的精品资源点击获取