简介这份资源是四川省自贡市30米分辨率的DEM数字高程数据包面向地理信息系统学习者、测绘与城市规划从业者以及开展地形分析的研究人员可用于地形建模、洪水模拟、地质灾害评估和城市空间分析等教学与实践场景。压缩包共12个文件约14.72MB核心为自贡市dem.tif高程栅格辅以ovr金字塔与tfw坐标文件便于快速浏览和配准自贡市范围.shp、shx、dbf、prj、sbn、sbx构成完整Shapefile边界数据另有多个xml元数据说明数据结构与属性。数据覆盖自贡市行政区域并延伸至部分邻近地区便于处理边界效应或扩展研究范围。目前已有421人学习下载适合需要真实区域高程底图、希望快速搭建GIS分析环境的初中级用户参考使用。1. 自贡市 30m DEM 到手之后先别急着往 GIS 里拖拿到「四川省自贡市DEM数字高程数据30m含本市级范围shp文件.zip」这类数据包很多人的第一反应是解压、拖进 QGIS 或 ArcGIS然后直接出图。我见过太多人卡在这一步要么坐标对不上要么范围裁歪了要么高程值全是异常。自贡这个地方地形很有代表性——西南丘陵区海拔落差从两百多米到一千米出头沱江穿城而过低山、丘陵、平坝交错。30m 分辨率的 DEM 在这个尺度上刚好够用既能看清丘陵起伏又不会让文件大到跑不动。这个数据包的核心价值在于「本市级范围 shp 文件」——它给了你一个现成的行政边界省去了自己去抠边界的麻烦。但边界文件和高程栅格能不能严丝合缝地对上取决于坐标系、裁剪方式和 NoData 处理。这篇笔记就按我实际处理这类数据的顺序把从解压到出可用成果的整条链路讲清楚适合做水文分析、选址评估、地形可视化的从业者照着复现。2. 解压后先看清三样东西栅格、边界、坐标系2.1 压缩包里通常有什么怎么快速判断能不能直接用这类数据包解压后一般包含三类文件DEM 栅格文件常见 .tif 或 .img 格式、本市级行政边界 shp 文件.shp/.shx/.dbf/.prj 一套、可能还有一个说明文档或元数据文件。先别打开 GIS 软件用命令行快速扫一遍文件结构和大小能省不少事。# 解压后进入目录先看文件清单和大小 unzip -l 四川省自贡市DEM数字高程数据30m含本市级范围shp文件.zip # 解压到指定目录 unzip 四川省自贡市DEM数字高程数据30m含本市级范围shp文件.zip -d ./zigong_dem # 查看解压后的文件结构 find ./zigong_dem -type f | head -50 ls -lh ./zigong_dem/逻辑说明unzip -l只列出压缩包内容不解压用来确认里面有没有你需要的 shp 和 tif。find加head是为了防止文件太多刷屏。ls -lh看文件大小30m 分辨率下自贡全市范围约 4381 平方公里的 DEM 栅格通常在几十到一百多 MB 之间如果只有几 MB可能是被压缩过或者范围不全。参数说明-d指定解压目录避免污染当前工作目录。如果压缩包内文件名含中文Linux 下可能需要加-O GBK或-O CP936参数解决乱码Windows 下用 7-Zip 或 Bandizip 更省心。2.2 坐标系不核对后面全白干这是血泪经验里排第一的坑。DEM 栅格和 shp 边界文件如果坐标系不一致你在 GIS 里看到的要么是边界飘到几千公里外要么是栅格和边界错位。自贡市常用的坐标系有两种地理坐标系 CGCS2000经纬度单位度和投影坐标系 CGCS2000 3 度带高斯-克吕格投影单位米自贡大概在 104°E 到 105°E 之间对应 3 度带第 35 带中央经线 105°E。# 用 gdalinfo 查看 DEM 栅格的坐标系和范围 gdalinfo ./zigong_dem/zigong_dem_30m.tif | grep -E Coordinate|Origin|Pixel|Upper|Lower # 用 ogrinfo 查看 shp 文件的坐标系 ogrinfo -al -so ./zigong_dem/zigong_boundary.shp | grep -E Extent|Layer|Geometry|Feature逻辑说明gdalinfo输出的Coordinate System字段告诉你栅格用的什么坐标系Origin和Pixel Size告诉你左上角坐标和像元大小。ogrinfo -al -so只输出摘要信息不打印所有要素Extent字段显示边界的范围。把两者的范围对比一下如果数量级差了几十倍基本可以确定一个是经纬度一个是投影坐标。参数说明-so是 summary only 的意思数据量大时必加否则终端会被要素信息刷爆。如果gdalinfo输出里Coordinate System显示GCS_China_Geodetic_Coordinate_System_2000说明是地理坐标系如果显示CGCS2000_3_Degree_GK_Zone_35或类似带GK字样的就是投影坐标系。提示如果两者坐标系不一致不要直接改栅格或 shp 的坐标系定义那是自欺欺人要用gdalwarp或ogr2ogr做真正的坐标转换。2.3 用 gdalwarp 统一坐标系的最小命令假设 DEM 是地理坐标系shp 是投影坐标系或者反过来统一到投影坐标系对后续面积、距离计算更友好。# 将 DEM 从地理坐标系转换到 CGCS2000 3 度带 35 带投影 gdalwarp -t_srs EPSG:4544 -r bilinear -of GTiff \ ./zigong_dem/zigong_dem_30m.tif \ ./zigong_dem/zigong_dem_30m_proj.tif # 如果 shp 是地理坐标系同样转换 ogr2ogr -f ESRI Shapefile -t_srs EPSG:4544 \ ./zigong_dem/zigong_boundary_proj.shp \ ./zigong_dem/zigong_boundary.shp逻辑说明gdalwarp的-t_srs指定目标坐标系-r bilinear是重采样方法DEM 连续数据用双线性插值比最近邻更平滑。ogr2ogr的-t_srs对矢量做投影转换。EPSG:4544 对应 CGCS2000 3 度带 35 带中央经线 105°E自贡市全域都落在这个带内。参数说明-r可选near最近邻适合分类数据、bilinear双线性适合连续数据、cubic三次卷积更平滑但可能过冲。DEM 用bilinear是常见做法。如果转换后栅格出现条带或异常值检查源数据是否有 NoData 值没设对。3. 按行政边界裁剪 DEMgdalwarp 和按掩膜提取怎么选3.1 两种裁剪方式的本质区别裁剪 DEM 到自贡市范围常见做法有两种一是用gdalwarp -cutline直接按 shp 边界裁二是先用gdal_translate裁矩形范围再用gdalwarp按掩膜提取。前者一步到位但边界外像元会被设为 NoData后者多一步但可控性更强。我一般会先用-cutline试一次如果边界复杂自贡下辖四区两县边界有飞地或狭长区域再考虑分步处理。关键参数是-crop_to_cutline不加这个参数的话输出栅格范围还是原始范围只是边界外变成 NoData文件大小没减。# 方式一一步裁剪到自贡市边界 gdalwarp -cutline ./zigong_dem/zigong_boundary_proj.shp \ -crop_to_cutline -dstnodata -9999 \ -r bilinear -of GTiff \ ./zigong_dem/zigong_dem_30m_proj.tif \ ./zigong_dem/zigong_dem_30m_clip.tif # 方式二先裁矩形再按掩膜提取适合边界复杂或需要分县处理 gdal_translate -projwin 104.5 29.5 105.5 28.8 \ -of GTiff \ ./zigong_dem/zigong_dem_30m_proj.tif \ ./zigong_dem/zigong_dem_30m_rect.tif gdalwarp -cutline ./zigong_dem/zigong_boundary_proj.shp \ -crop_to_cutline -dstnodata -9999 \ ./zigong_dem/zigong_dem_30m_rect.tif \ ./zigong_dem/zigong_dem_30m_clip.tif逻辑说明-cutline指定裁剪边界 shp-crop_to_cutline让输出范围贴合边界外接矩形-dstnodata -9999把边界外像元设为 -9999方便后续识别。方式二的-projwin参数顺序是xmin ymax xmax ymin注意是「左上右下」不是「左下右上」这个顺序搞反了裁出来是空的。参数说明-dstnodata的值要和后续分析工具兼容ArcGIS 里常用 -9999QGIS 里也可以用 -9999 或 0但 0 在 DEM 里可能是真实高程自贡最低点约 240m所以别用 0。-projwin的坐标要跟栅格坐标系一致投影坐标下单位是米地理坐标下是度。3.2 裁剪后必做的三项检查裁完不是就完事了至少检查三样范围对不对、NoData 设没设对、高程值有没有异常。# 检查裁剪后栅格的范围和 NoData 值 gdalinfo ./zigong_dem/zigong_dem_30m_clip.tif | grep -E Upper|Lower|NoData|Size # 统计高程值分布排除 NoData gdalinfo -stats ./zigong_dem/zigong_dem_30m_clip.tif | grep -A 20 STATISTICS # 用 Python 快速检查异常值 python3 -c from osgeo import gdal import numpy as np ds gdal.Open(./zigong_dem/zigong_dem_30m_clip.tif) band ds.GetRasterBand(1) arr band.ReadAsArray() nodata band.GetNoDataValue() valid arr[arr ! nodata] print(f有效像元数: {valid.size}) print(f高程范围: {valid.min():.1f} ~ {valid.max():.1f} 米) print(f均值: {valid.mean():.1f} 米) print(fNoData 像元数: {arr.size - valid.size}) 逻辑说明gdalinfo -stats会计算栅格的统计信息包括最小值、最大值、均值、标准差。自贡市高程范围大概在 240m 到 1000m 出头如果统计出来最小值是 -9999 或最大值是 9999说明 NoData 没设对或者有异常值。Python 脚本用numpy过滤 NoData 后统计更直观。参数说明GetNoDataValue()返回栅格设置的 NoData 值如果返回None说明没设需要手动处理。arr ! nodata做布尔索引注意如果 nodata 是浮点数直接用!比较可能有精度问题稳妥做法是用np.isclose或设一个容差范围。注意如果裁剪后有效像元数远小于预期自贡全市约 4381 平方公里30m 像元约 487 万个检查是不是-crop_to_cutline没加或者 shp 边界本身有问题比如坐标系不对导致边界跑到别处。4. 避坑与排查自贡 DEM 处理中最容易翻车的五个地方4.1 坑一shp 边界和 DEM 范围对不上差了几十公里现象在 QGIS 里同时加载 DEM 和 shp发现边界飘在栅格外面或者只覆盖了栅格的一个角。原因最常见的是坐标系不一致。DEM 是地理坐标系经纬度shp 是投影坐标系米或者反过来。另一个可能是 shp 文件本身的范围就是错的比如从某个在线地图下载的边界坐标系标的是 WGS84 但实际是 GCJ02 偏移过的。解决先用gdalinfo和ogrinfo分别看两者的坐标系和范围。如果坐标系不一致用gdalwarp和ogr2ogr统一。如果坐标系标称一致但范围还是对不上用 QGIS 的「缩放到图层」功能分别看两个图层的实际位置确认是不是数据本身有偏移。自贡地区如果用到从某些在线地图获取的边界注意 GCJ02 偏移问题需要做坐标纠偏。4.2 坑二裁剪后栅格全是 NoData 或者只有一条边有数据现象执行gdalwarp -cutline后输出栅格大部分是 NoData只有边缘一小条有数据。原因-cutline的 shp 坐标系和输入栅格坐标系不一致gdalwarp 按 shp 的坐标去裁栅格但两者不在一个空间参考下导致裁剪区域错位。另一个可能是 shp 的几何类型有问题比如是线而不是面或者面有自相交。解决确认 shp 和栅格坐标系一致后再裁。用ogrinfo -al -so看 shp 的Geometry字段必须是Polygon或MultiPolygon。如果是线需要用ogr2ogr或 QGIS 转成面。面有自相交的话用 QGIS 的「修复几何」工具处理。4.3 坑三高程值出现负值或异常大值现象统计高程时发现最小值是 -9999 或 -32768或者最大值是 9999。原因-9999 通常是 NoData 值没被正确识别-32768 是某些格式如 ERDAS IMG的默认 NoData。9999 可能是原始数据里的填充值或错误值。解决用gdalwarp的-dstnodata重新指定 NoData 值或者在 Python 里用numpy过滤。如果原始数据本身就有异常值需要用gdal_calc或 Python 做条件替换。# 用 Python 将异常值替换为 NoData from osgeo import gdal import numpy as np ds gdal.Open(./zigong_dem/zigong_dem_30m_clip.tif, gdal.GA_Update) band ds.GetRasterBand(1) arr band.ReadAsArray() # 将小于 0 或大于 2000 的值设为 NoData arr[(arr 0) | (arr 2000)] -9999 band.SetNoDataValue(-9999) band.WriteArray(arr) ds None逻辑说明自贡市真实高程不会低于 0 米也不会高于 2000 米用这个范围做过滤是合理的。GA_Update以可写模式打开SetNoDataValue设置 NoData 值WriteArray写回。最后ds None关闭数据集确保数据落盘。参数说明阈值 0 和 2000 是根据自贡实际地形定的如果换到其他地区要调整。写回前建议先备份原始文件避免误操作覆盖。4.4 坑四裁剪后文件太大跑不动现象自贡全市 30m DEM 裁剪后文件还有几百 MB在 QGIS 里缩放卡顿。原因GeoTIFF 默认不压缩30m 分辨率下像元多文件自然大。另外如果 NoData 区域没有用掩膜或压缩存储效率低。解决用gdal_translate加压缩参数重新输出。# 用 LZW 压缩和金字塔重采样减小文件 gdal_translate -co COMPRESSLZW -co PREDICTOR2 \ -co TILEDYES -co BIGTIFFIF_SAFER \ ./zigong_dem/zigong_dem_30m_clip.tif \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif # 添加金字塔加速缩放显示 gdaladdo -r average ./zigong_dem/zigong_dem_30m_clip_compressed.tif 2 4 8 16逻辑说明COMPRESSLZW是无损压缩PREDICTOR2对连续数据如 DEM压缩率更好TILEDYES分块存储加速读取BIGTIFFIF_SAFER在文件可能超过 4GB 时自动用 BigTIFF 格式。gdaladdo添加金字塔QGIS 和 ArcGIS 缩放时会自动调用显示更流畅。参数说明PREDICTOR可选 1不预测、2水平差分、3浮点预测DEM 用 2 或 3 都行。金字塔层级2 4 8 16表示 1/2、1/4、1/8、1/16 分辨率一般加到 16 或 32 就够。4.5 坑五用 ArcGIS 裁剪时结果和 gdalwarp 不一致现象同样的 shp 和 DEMArcGIS 的「按掩膜提取」和 gdalwarp 裁出来边界处像元值不一样。原因两者的重采样默认方法和 NoData 处理逻辑不同。ArcGIS 默认可能用最近邻gdalwarp 默认用最近邻但可以指定双线性。边界处像元如果跨在裁剪线上不同方法取值不同。解决统一重采样方法。gdalwarp 加-r bilinearArcGIS 里在环境设置中把「重采样技术」改为「双线性」或「三次卷积」。另外确认两者的 NoData 值设置一致。如果做定量分析建议全程用同一套工具链别混用。5. 从 DEM 到可用成果坡度坡向提取与水文分析的参数怎么设5.1 坡度坡向提取gdal 和 ArcGIS 的参数差异DEM 最常用的衍生成果是坡度和坡向。自贡丘陵区坡度分析对农业选址、水土保持很有价值。用gdaldem命令行提取坡度参数设置直接影响结果。# 提取坡度度为单位 gdaldem slope -of GTiff -compute_edges \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif \ ./zigong_dem/zigong_slope.tif # 提取坡向 gdaldem aspect -of GTiff -compute_edges \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif \ ./zigong_dem/zigong_aspect.tif # 提取山体阴影可视化用 gdaldem hillshade -of GTiff -z 2 -az 315 -alt 45 \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif \ ./zigong_dem/zigong_hillshade.tif逻辑说明gdaldem slope默认输出度加-p输出百分比坡度。-compute_edges让边缘像元也参与计算不加的话边缘一圈是 NoData。hillshade的-z 2是垂直夸张系数自贡地形起伏不大用 2 到 3 比较合适-az 315是光源方位角西北方向-alt 45是光源高度角。参数说明坡度提取的算法基于 Horn 方法3x3 窗口ArcGIS 的「坡度」工具默认也是 Horn 方法但 ArcGIS 会先做边缘填充。如果两者结果在边缘处不一致是正常现象。坡向输出 0-360 度0 为正北90 为正东ArcGIS 的坡向输出范围是 -1 到 360-1 表示平地。5.2 水文分析填洼、流向、流量累积的关键参数自贡有沱江及其支流做水文分析前必须填洼否则流向计算会断。# 用 gdal_fillnodata 填洼简单场景 gdal_fillnodata.py -md 10 -si 0 \ ./zigong_dem/zigong_dem_30m_clip_compressed.tif \ ./zigong_dem/zigong_dem_filled.tif # 用 Python RichDEM 做更专业的填洼和流向分析 python3 -c import richdem as rd dem rd.LoadGDAL(./zigong_dem/zigong_dem_30m_clip_compressed.tif) dem_filled rd.FillDepressions(dem, epsilonTrue, in_placeFalse) flow_accum rd.FlowAccumulation(dem_filled, methodD8) rd.SaveGDAL(./zigong_dem/zigong_flow_accum.tif, flow_accum) 逻辑说明gdal_fillnodata.py的-md 10是最大搜索距离像元数-si 0是搜索步长。RichDEM 的FillDepressions用epsilonTrue做微填洼避免大范围平坦区域导致流向不确定。FlowAccumulation用 D8 算法八方向适合丘陵区。参数说明填洼的-md参数根据洼地大小调自贡丘陵区一般 10 到 20 够用。RichDEM 的epsilon填洼会在平坦区加微小梯度保证流向唯一。流量累积结果中高值对应河道可以用来提取河网阈值一般设 1000 到 5000 个像元30m 分辨率下约 0.9 到 4.5 平方公里汇水面积。提示如果做正式水文分析建议用 ArcGIS 的 Hydrology 工具箱或 WhiteboxToolsRichDEM 适合快速验证。不同工具的填洼算法有差异结果会有细微不同选一套用到底就行。5.3 用自贡 shp 边界做分区统计的实操拿到坡度、坡向后常需要按行政区统计。用zonal工具或 Python 的rasterstats库。# 用 rasterstats 按自贡市边界统计坡度 from rasterstats import zonal_stats import geopandas as gpd boundary gpd.read_file(./zigong_dem/zigong_boundary_proj.shp) stats zonal_stats(boundary, ./zigong_dem/zigong_slope.tif, stats[min, max, mean, median, std], nodata-9999) print(stats)逻辑说明zonal_stats接受矢量边界和栅格路径stats指定要计算的统计量。nodata参数要和栅格实际 NoData 值一致否则会把 NoData 算进去。输出是一个列表每个要素对应一个字典。参数说明如果边界有多个要素比如自贡下辖的区县结果会按要素顺序返回。all_touchedTrue参数会让边界接触到的所有像元都参与统计默认只统计中心点在边界内的像元。做精确面积统计时用all_touchedTrue更合理。6. 一个容易被忽略的技巧用 DEM 做自贡市地形起伏度分级地形起伏度是描述地表切割程度的指标定义为特定窗口内高程最大值与最小值之差。自贡丘陵区用这个指标做地貌分级比单纯看坡度更直观。我一般用 Python 加scipy做滑动窗口计算窗口大小取 3x3 到 15x15 像元对应 90m 到 450m 的空间尺度。from osgeo import gdal import numpy as np from scipy.ndimage import maximum_filter, minimum_filter # 读取 DEM ds gdal.Open(./zigong_dem/zigong_dem_30m_clip_compressed.tif) band ds.GetRasterBand(1) dem band.ReadAsArray().astype(np.float32) nodata band.GetNoDataValue() # 将 NoData 设为 NaN 避免影响计算 dem[dem nodata] np.nan # 计算 5x5 窗口150m 尺度的地形起伏度 window_size 5 max_dem maximum_filter(dem, sizewindow_size, modenearest) min_dem minimum_filter(dem, sizewindow_size, modenearest) relief max_dem - min_dem # 按起伏度分级参考自贡实际地形 # 30m 平坝30-70m 浅丘70-150m 中丘 150m 深丘/低山 classified np.zeros_like(relief, dtypenp.uint8) classified[relief 30] 1 classified[(relief 30) (relief 70)] 2 classified[(relief 70) (relief 150)] 3 classified[relief 150] 4 classified[np.isnan(relief)] 0 # 输出分级结果 driver gdal.GetDriverByName(GTiff) out_ds driver.Create(./zigong_dem/zigong_relief_class.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_band out_ds.GetRasterBand(1) out_band.WriteArray(classified) out_band.SetNoDataValue(0) out_ds None逻辑说明maximum_filter和minimum_filter是scipy.ndimage提供的滑动窗口极值滤波size5表示 5x5 窗口modenearest让边缘像元用最近的有效值填充。relief就是起伏度。分级阈值 30m、70m、150m 是根据自贡实际地形定的——平坝区起伏度通常小于 30m浅丘 30 到 70m中丘 70 到 150m深丘和低山大于 150m。参数说明窗口大小决定分析尺度5x5 对应 150m 见方适合村级规划15x15 对应 450m 见方适合乡镇级。分级阈值可以根据具体项目调整但建议先统计relief的分位数比如 25%、50%、75% 分位再定阈值比拍脑袋更靠谱。输出用GDT_Byte节省空间NoData 设为 0。这个方法的实用价值在于自贡很多区域坡度看着不大但起伏度很高说明是破碎丘陵修路和选址成本比平坝高得多。把起伏度分级图和 shp 边界叠加能快速识别哪些乡镇以平坝为主、哪些以深丘为主。我一般会把这个结果和坡度图交叉坡度大于 15 度且起伏度大于 70m 的区域标记为「建设难度高」给规划部门做参考。做这类分析最大的教训是别一上来就追求花哨的算法先把坐标系、NoData、裁剪范围这三样核对清楚后面所有分析都顺。我见过太多人卡在坐标系上折腾一整天最后发现只是 shp 的 .prj 文件缺失。拿到数据先gdalinfo和ogrinfo扫一遍花不了五分钟能省几小时。希望帮到你。本文还有配套的精品资源点击获取