简介一份覆盖内蒙古乌兰察布市全域的30米分辨率DEM数字高程数据包除核心高程栅格外还附带本市级行政范围Shapefile边界便于直接裁剪与分析。数据适用于国土规划、生态评估、洪涝模拟、通信基站选址等GIS场景也适合GIS初学者以市级尺度上手地形数据处理。压缩包共12个文件约220.2MB以GeoTIFF高程数据、市级范围矢量文件SHP及配套DBF、SHX、SBN/SBX索引为主另含PRJ投影定义、TFW世界文件与XML元数据可在ArcGIS、QGIS等平台中完整解析空间参考。矢量边界可为后续坡度、坡向、可视域分析提供统一掩膜范围配套投影文件则确保坐标精度与多源数据叠加的一致性。该数据包已有253人学习下载。获取后可快速在GIS软件中加载dem.tif进行地形渲染并结合范围Shp实现批量裁剪、等高线提取、流域分割等二次开发省去自行爬取与拼接大范围DEM的流程。1. 乌兰察布市30米DEM数据拿到手后先别急着叠加把压缩包解压成十几个文件随手把乌兰察布市dem.tif拖进 ArcGIS第一眼看到的多半是一片灰黑相间的色块。这不是数据坏了而是你还没做符号化拉伸。真正影响后续分析的关键问题是这个 DEM 与附带的乌兰察布市范围.shp是否处于同一坐标系以及 DEM 的像元落入坐标能不能精确对齐边界。这篇文不打算给你重讲 DEM 原理而是针对这个市本级 30 米分辨率的数据压缩包把 tif、tfw、ovr、shp 这一组文件的读取顺序、坐标系检查、裁剪与坡度分析整套流程跑通。适合的读者是常年与栅格数据打交道的 GIS 工程师或者接到活儿需要快速用地形数据做评估的人。2. 解压后的文件家族tif、tfw、ovr、shp 各自的角色2.1 从后缀看数据血缘哪个文件缺了会出问题解压后你看到那一堆.sbx、.ovr、.tfw不是冗余垃圾而是 GIS 软件在不同阶段要用到的配套件。我习惯先把它们分成两类一类是给栅格用的一类是给矢量边界用的。栅格侧以乌兰察布市dem.tif为核心tfw是纯文本的世界文件ovr是金字塔概览xml里记录着处理管线矢量侧则是一整套标准 Shapefileshp存几何、dbf存属性、prj存投影shx/sbn/sbx是空间索引。理解每个文件的角色排错时才能少走弯路。文件类型作用缺失后果乌兰察布市dem.tif栅格高程主数据像元值代表海拔米数没有它就没有 DEM乌兰察布市dem.tfw文本仿射变换参数记录像元大小与左上角坐标部分软件无法定位 tif拖入后显示在 0,0 附近乌兰察布市dem.tif.ovr栅格概览重采样金字塔加速缩放显示大范围浏览变慢分析不受影响乌兰察布市dem.tif.xml元数据投影历史、处理时间戳排查坐标异常时会漏掉关键线索乌兰察布市范围.shp矢量面市域行政边界用于裁剪与掩膜无法做市域精确裁剪乌兰察布市范围.dbf属性表行政名称、代码等属性shp 属性信息丢失但几何还在乌兰察布市范围.prj投影文本用 WKT 描述矢量坐标系软件按默认坐标系猜测容易叠加错位乌兰察布市范围.shp.xml元数据可选字段说明一般不影响使用.tfw里每一行都有固定语义第 1、4 行是像元宽高第 5、6 行是左上角像素中心的坐标。很多新手以为把tfw删掉只留tif也能用确实能但代价是 ArcMap 这类传统桌面软件会丢失地理配准信息只能用GeoTIFF头里嵌入的坐标。如果这份 tif 本身已经内置投影信息那tfw更像一个备份。但要留意tfw的坐标约定和 GDAL 的GetGeoTransform返回值之间有时差半像素这是二十年前的历史遗留不算 bug。2.2 用 GDAL 验证 tif 与 shp 的坐标系和覆盖范围我会在动手做任何分析之前先跑两条命令确认数据状态。GDAL 是处理这类问题的首选Windows 下用 OSGeo4W ShellLinux 直接终端macOS 用 Homebrew 装gdal即可。gdalinfo -stats 乌兰察布市dem.tif | grep -E Origin|Pixel Size|Band 1|STATISTICS|Coordinate System-stats会强制计算像元统计值并写入附属文件这样你能在输出里看到Minimum、Maximum、Mean和StdDev。grep -E只是过滤避免刷屏。输出的Coordinate System is:这一行最关键比如它写成WGS 84 / UTM zone 49N而你手里的乌兰察布市范围.shp是CGCS2000 / 3-degree Gauss-Kruger CM 111E那后续所有叠加都必须在同一基准下完成。再看Pixel Size是不是30, -30如果你看到0.0003, -0.0003这种带小数的值说明这个 tif 可能被重投影过不再是原始的 30 米地理分辨率。接着检查 shpogrinfo -so -al 乌兰察布市范围.shp-so表示只显示摘要-al遍历所有图层。输出里有Extent、Geometry、Layer SRS把范围与 tif 的坐标范围做对比就能判断 DEM 是否完整覆盖市域。通常 DEM 会向边界外扩几十个像元这是图幅拼接的正常现象不用意外。如果 DEM 范围反而小于 shp 范围裁剪前必须用gdal_translate重新整编边缘区域。提示tif 与 shp 坐标系不一致时绝对不要直接在 ArcMap 里开启动态投影做“表面上的对齐”。动态投影只管显示实际分析时仍按原始坐标计算可能造成几公里级别的空间错位。2.3 从 tfw 反推像元坐标检查是否被暴力二次压缩有时你拿到的 tif 已经被工具改过tfw可能落后于 tif 内部头或者干脆是错的。更可靠的做法是绕过tfw直接读 GeoTIFF 里的仿射变换。下面这段 Python 用 GDAL 绑定完成from osgeo import gdal ds gdal.Open(乌兰察布市dem.tif) gt ds.GetGeoTransform() w, h ds.RasterXSize, ds.RasterYSize minx gt[0] maxx gt[0] w * gt[1] h * gt[2] maxy gt[3] miny gt[3] w * gt[4] h * gt[5] print(fleft{minx:.2f}, bottom{miny:.2f}, right{maxx:.2f}, top{maxy:.2f}) print(f像元宽{gt[1]:.4f}, 像元高{gt[5]:.4f})gt[0]是左上角 X 坐标gt[1]是水平像元分辨率gt[2]通常为 0 表示不旋转gt[4]也类似。gt[3]是左上角 Y 坐标gt[5]为负值表示栅格自上而下排列。如果你发现gt[5]为正数说明数据来自非标准软件QGIS 还能容忍ArcGIS 画出来可能是上下颠倒的。另外检查像元宽高是否约等于 30 米在投影坐标系下或 0.00027 度在 WGS84 下。如果写的是 0.00007 度那这份 DEM 实际上可能被重采样到约 7.8 米分辨率和“30米”并不相符。3. 在 QGIS 与 ArcMap 中加载符号化先调渲染再谈分析3.1 单波段假彩色与直方图拉伸把乌兰察布市dem.tif直接拖进 QGIS默认用灰度渲染地形起伏几乎看不出来。正确做法是在图层属性的“符号化”里选择“单波段假彩色”配色方案选terrain或RdYlGr然后在“具色渐变”的设置里把最小/最大值裁剪到 2% 到 98%。这样做的原因是 DEM 数据里常存在极端的局部空洞像元比如噪点或云影残留它们会拉宽直方图导致主要地貌的色带被压缩在中间段。ArcMap 则是在图层属性“符号化”里选择“拉伸”拉伸类型用“百分比截断”截断量填 1~2%色带选择自己习惯的。如果你用 QGIS 的 Python 控制台可以用一行脚本快速加载并设置渲染iface.addRasterLayer(乌兰察布市dem.tif, dem)但注意这行代码只负责加载图层并没有设置拉伸参数。更完整的方式是通过QgsSingleBandPseudoColorRenderer创建渲染器但日常使用中手动点几下比写脚本更快除非你有超过几十个图层需要批量处理。软件推荐拉伸方式参数参考适用场景QGIS百分比截断截断量 2% 到 98%快速浏览大范围地形ArcMap标准差拉伸2 到 3 倍标准差保留局部起伏细节命令行 gdaldem颜色映射自定义色带文件出图用不受软件限制3.2 先对齐投影再做掩膜不要依赖动态投影很多文档会教你在 ArcToolbox 里用“按掩膜提取”直接裁剪但那个工具对输入之间的空间参考一致性要求极高。如果你不先对齐工具大概率会弹出“几何坐标不一致”的警告然后给你一个空结果或错位结果。我的工作流里第一步永远是ogr2ogr把边界 shp 转换为与 tif 相同的投影gdalinfo -proj4 乌兰察布市dem.tif这条命令会把投影参数转成 proj4 形式然后取它的 EPSG 编号比如 CGCS2000 UTM 49N。接下来ogr2ogr -t_srs EPSG:xxxx -overwrite 乌兰察布市范围_proj.shp 乌兰察布市范围.shp参数说明-t_srs指定目标投影-overwrite允许覆盖中间文件。转换后最好再跑一次ogrinfo -so -al确认 Extent 已经变化到新坐标系下这样gdalwarp裁剪时不会因为边界和 DEM 范围差异而出现边缘黑带。3.3 用 gdalwarp 一次性完成重投影与裁剪裁剪 DEM 常见做法是执行gdalwarp它同时内置了重投影能力所以投影对齐和裁剪可以一步完成不必像文档流程那样分两步走。我经常对乌兰察布这种市级数据做“投影转换 边界裁剪”的复合操作gdalwarp -t_srs EPSG:4326 \ -cutline 乌兰察布市范围_proj.shp \ -crop_to_cutline \ -dstalpha \ -tr 0.0003 0.0003 \ -r cubic \ 乌兰察布市dem.tif 乌兰察布市dem_4326_masked.tif-t_srs强制输出到 WGS84 经纬度-cutline指定边界矢量-crop_to_cutline让输出范围严格贴合边界外接矩形顺便把边界外的像元全部裁掉-dstalpha增加一个 alpha 通道边界外透明显示-tr是目标像元尺寸0.0003 度约等于 33 米如果仍想要 30 米可以算准确值或使用-te指定输出范围-r cubic用三次卷积重采样在一定程度上平滑锯齿适合 DEM。如果你要保留原始高程统计特性也可以改用-r bilinear。注意-crop_to_cutline裁剪的是整个边界外接矩形而不是按边界形状做非矩形裁剪。边界外区域虽然保留但会被 alpha 通道遮罩。要得到完全按面要素裁剪的栅格后续还需要gdal_rasterize做掩膜。4. 坡度、坡向与地形起伏度Z因子和重采样方法的选择4.1 填洼不是DEM分析的第一站我经常看到有人一上来就做 Fill Sinks认为 DEM 必须一律先填洼。这是个误解。填洼的目的是消除伪洼地从而获得连续水流方向只有水文分析才需要这一步。如果只是做坡度、坡向、地形起伏度或视域分析填洼反而会抹掉真实的地貌细节比如人工梯田、天然坑塘、采石场坑底。特别是在乌兰察布这类草原与农牧交错区微地形对地表径流和植被分布的影响很明显先算坡度坡向更合理。4.2 gdaldem 计算坡度-p 和 -s 参数的正确姿势GDAL 自带gdaldem工具一行就能生成坡度图gdaldem slope 乌兰察布市dem.tif 乌兰察布市slope.tif -p -s 111120-p控制输出格式加上它输出百分数坡度rise over run最高可超几百不加则输出角度制坡度0–90 度。-s是水平距离的缩放因子直接决定结果是否准确它只在输入为地理坐标系度时才有意义。因为 tif 内部高程单位是米而水平单位是度不设置-sGDAL 会默认水平单位等于高程单位导致坡度几乎全部接近 90 度。这里的111120表示 1 度约等于 111120 米这是全球平均估算值。如果你想更精确就要按当地纬度余弦校正乌兰察布大约北纬 41~43 度cos(42°)≈0.743所以更贴近实际的缩放因子是 111120×0.743≈82562。不过要注意这样做会引入纬度相关的形变不推荐在对精度有要求的科研项目里用。最稳妥的做法是先把 DEM 投影到等面积或 UTM 投影然后让水平和垂直单位都为米此时不需要-s。参数作用常见坑-p输出百分数坡度不加则输出角度值-s地理坐标系下的缩放因子忘记加或加错会得到离谱坡度-r重采样方法默认 nearest 在 DEM 上不推荐-z垂直夸张因子仅用于视觉效果4.3 用 Python rasterio 批量统计坡度分区面积用 gdaldem 生成坡度图只是第一步分析时经常需要知道不同坡度级别的面积占比。用 QGIS 的“栅格计算器”也能做但脚本化更适合多期数据对比。下面这段代码用 rasterio 读取坡度图按国标坡度分级统计面积import rasterio import numpy as np import pandas as pd bins [-1, 5, 15, 25, 35, 90] labels [0-5°, 5-15°, 15-25°, 25-35°, 35°] with rasterio.open(乌兰察布市slope.tif) as src: data src.read(1) nodata src.nodata data np.where(data nodata, np.nan, data) cell_area abs(src.transform[0] * src.transform[4]) flat data.flatten() flat flat[~np.isnan(flat)] counts, _ np.histogram(flat, binsbins) df pd.DataFrame({ grade: labels, area_km2: np.round(counts * cell_area / 1e6, 2) }) df.to_csv(slope_report.csv, indexFalse) print(df)参数说明src.transform[0]是像元宽度src.transform[4]通常为负值所以取绝对值得到面积。如果之前做过投影转换这段脚本能直接吃到 30 米或 33 米的像元面积。分箱边界bins可根据自己的目标改成 0-3°、3-8° 等。输出 CSV 可以直接用 Excel 或后续 Python 画图。5. 一条命令与一个瘦脚本把30米DEM变成可交付的裁剪结果5.1 从 shp 到 KML边界预览与检查在做数据交付时对方往往不想打开 QGIS而是想快速在网页或谷歌地图里看一眼边界。最轻量的是从 Shapefile 转换成 KMLogr2ogr -f KML 乌兰察布市范围.kml 乌兰察布市范围.shp如果你已经做了投影转换用_proj.shp转如果没有转换KML 里会保留原始投影信息但 Google Earth 会按 WGS84 勉强显示可能出现边缘不完全贴合。建议先转成 EPSG:4326 再输出。5.2 用 gdal.Warp 批量裁剪并自动识别投影一致性当你手里有多个 DEM 文件需要按同一边界裁剪时GIS 菜单反复点击的效率太低。我一般用一段 Python 脚本直接从 SDL 调用 GDAL 的WarpAPI。这样既避开了命令行拼接文件名的坑也能在脚本里加投影一致性检查。from osgeo import gdal, ogr import glob cutline 乌兰察布市范围_proj.shp tif_list glob.glob(*.tif) # 先读取边界投影 ref ogr.Open(cutline) ref_lsr ref.GetLayer(0).GetSpatialRef() for tif in tif_list: ds gdal.Open(tif) src_sr ds.GetProjection() if not src_sr and ref_lsr is None: print(f{tif}: 缺少投影信息跳过) out tif.replace(.tif, _masked.tif) gdal.Warp(out, tif, cutlineDSNamecutline, cropToCutlineTrue, dstNodata-9999, options[-overwrite]) print(f{tif} - {out})这里gdal.Warp的第一个参数是输出路径第二个是输入cutlineDSName对应命令行-cutlinecropToCutline对应-crop_to_cutlinedstNodata给输出栅格设定无效值。options[-overwrite]允许覆盖已有文件。代码里的投影检查只是最基础一层实际还可以加一句“如果输入与边界的 EPSG 不一致就先用gdal.Warp做一次srcSRS到边界投影的转换”但更干净的做法是在前面章节里先把 tif 统一成 CGCS2000 / 3-degree Gauss-Kruger再跑这个批量裁剪。批量处理完后用gdalinfo -stats抽查一个新输出文件确认Minimum不等于-9999否则说明边界内部可能还有无效像元残留需要补做填补或参考相邻像元处理。本文还有配套的精品资源点击获取 30米DEM数据处理：GDAL坐标校验与坡度分析全流程
简介一份覆盖内蒙古乌兰察布市全域的30米分辨率DEM数字高程数据包除核心高程栅格外还附带本市级行政范围Shapefile边界便于直接裁剪与分析。数据适用于国土规划、生态评估、洪涝模拟、通信基站选址等GIS场景也适合GIS初学者以市级尺度上手地形数据处理。压缩包共12个文件约220.2MB以GeoTIFF高程数据、市级范围矢量文件SHP及配套DBF、SHX、SBN/SBX索引为主另含PRJ投影定义、TFW世界文件与XML元数据可在ArcGIS、QGIS等平台中完整解析空间参考。矢量边界可为后续坡度、坡向、可视域分析提供统一掩膜范围配套投影文件则确保坐标精度与多源数据叠加的一致性。该数据包已有253人学习下载。获取后可快速在GIS软件中加载dem.tif进行地形渲染并结合范围Shp实现批量裁剪、等高线提取、流域分割等二次开发省去自行爬取与拼接大范围DEM的流程。1. 乌兰察布市30米DEM数据拿到手后先别急着叠加把压缩包解压成十几个文件随手把乌兰察布市dem.tif拖进 ArcGIS第一眼看到的多半是一片灰黑相间的色块。这不是数据坏了而是你还没做符号化拉伸。真正影响后续分析的关键问题是这个 DEM 与附带的乌兰察布市范围.shp是否处于同一坐标系以及 DEM 的像元落入坐标能不能精确对齐边界。这篇文不打算给你重讲 DEM 原理而是针对这个市本级 30 米分辨率的数据压缩包把 tif、tfw、ovr、shp 这一组文件的读取顺序、坐标系检查、裁剪与坡度分析整套流程跑通。适合的读者是常年与栅格数据打交道的 GIS 工程师或者接到活儿需要快速用地形数据做评估的人。2. 解压后的文件家族tif、tfw、ovr、shp 各自的角色2.1 从后缀看数据血缘哪个文件缺了会出问题解压后你看到那一堆.sbx、.ovr、.tfw不是冗余垃圾而是 GIS 软件在不同阶段要用到的配套件。我习惯先把它们分成两类一类是给栅格用的一类是给矢量边界用的。栅格侧以乌兰察布市dem.tif为核心tfw是纯文本的世界文件ovr是金字塔概览xml里记录着处理管线矢量侧则是一整套标准 Shapefileshp存几何、dbf存属性、prj存投影shx/sbn/sbx是空间索引。理解每个文件的角色排错时才能少走弯路。文件类型作用缺失后果乌兰察布市dem.tif栅格高程主数据像元值代表海拔米数没有它就没有 DEM乌兰察布市dem.tfw文本仿射变换参数记录像元大小与左上角坐标部分软件无法定位 tif拖入后显示在 0,0 附近乌兰察布市dem.tif.ovr栅格概览重采样金字塔加速缩放显示大范围浏览变慢分析不受影响乌兰察布市dem.tif.xml元数据投影历史、处理时间戳排查坐标异常时会漏掉关键线索乌兰察布市范围.shp矢量面市域行政边界用于裁剪与掩膜无法做市域精确裁剪乌兰察布市范围.dbf属性表行政名称、代码等属性shp 属性信息丢失但几何还在乌兰察布市范围.prj投影文本用 WKT 描述矢量坐标系软件按默认坐标系猜测容易叠加错位乌兰察布市范围.shp.xml元数据可选字段说明一般不影响使用.tfw里每一行都有固定语义第 1、4 行是像元宽高第 5、6 行是左上角像素中心的坐标。很多新手以为把tfw删掉只留tif也能用确实能但代价是 ArcMap 这类传统桌面软件会丢失地理配准信息只能用GeoTIFF头里嵌入的坐标。如果这份 tif 本身已经内置投影信息那tfw更像一个备份。但要留意tfw的坐标约定和 GDAL 的GetGeoTransform返回值之间有时差半像素这是二十年前的历史遗留不算 bug。2.2 用 GDAL 验证 tif 与 shp 的坐标系和覆盖范围我会在动手做任何分析之前先跑两条命令确认数据状态。GDAL 是处理这类问题的首选Windows 下用 OSGeo4W ShellLinux 直接终端macOS 用 Homebrew 装gdal即可。gdalinfo -stats 乌兰察布市dem.tif | grep -E Origin|Pixel Size|Band 1|STATISTICS|Coordinate System-stats会强制计算像元统计值并写入附属文件这样你能在输出里看到Minimum、Maximum、Mean和StdDev。grep -E只是过滤避免刷屏。输出的Coordinate System is:这一行最关键比如它写成WGS 84 / UTM zone 49N而你手里的乌兰察布市范围.shp是CGCS2000 / 3-degree Gauss-Kruger CM 111E那后续所有叠加都必须在同一基准下完成。再看Pixel Size是不是30, -30如果你看到0.0003, -0.0003这种带小数的值说明这个 tif 可能被重投影过不再是原始的 30 米地理分辨率。接着检查 shpogrinfo -so -al 乌兰察布市范围.shp-so表示只显示摘要-al遍历所有图层。输出里有Extent、Geometry、Layer SRS把范围与 tif 的坐标范围做对比就能判断 DEM 是否完整覆盖市域。通常 DEM 会向边界外扩几十个像元这是图幅拼接的正常现象不用意外。如果 DEM 范围反而小于 shp 范围裁剪前必须用gdal_translate重新整编边缘区域。提示tif 与 shp 坐标系不一致时绝对不要直接在 ArcMap 里开启动态投影做“表面上的对齐”。动态投影只管显示实际分析时仍按原始坐标计算可能造成几公里级别的空间错位。2.3 从 tfw 反推像元坐标检查是否被暴力二次压缩有时你拿到的 tif 已经被工具改过tfw可能落后于 tif 内部头或者干脆是错的。更可靠的做法是绕过tfw直接读 GeoTIFF 里的仿射变换。下面这段 Python 用 GDAL 绑定完成from osgeo import gdal ds gdal.Open(乌兰察布市dem.tif) gt ds.GetGeoTransform() w, h ds.RasterXSize, ds.RasterYSize minx gt[0] maxx gt[0] w * gt[1] h * gt[2] maxy gt[3] miny gt[3] w * gt[4] h * gt[5] print(fleft{minx:.2f}, bottom{miny:.2f}, right{maxx:.2f}, top{maxy:.2f}) print(f像元宽{gt[1]:.4f}, 像元高{gt[5]:.4f})gt[0]是左上角 X 坐标gt[1]是水平像元分辨率gt[2]通常为 0 表示不旋转gt[4]也类似。gt[3]是左上角 Y 坐标gt[5]为负值表示栅格自上而下排列。如果你发现gt[5]为正数说明数据来自非标准软件QGIS 还能容忍ArcGIS 画出来可能是上下颠倒的。另外检查像元宽高是否约等于 30 米在投影坐标系下或 0.00027 度在 WGS84 下。如果写的是 0.00007 度那这份 DEM 实际上可能被重采样到约 7.8 米分辨率和“30米”并不相符。3. 在 QGIS 与 ArcMap 中加载符号化先调渲染再谈分析3.1 单波段假彩色与直方图拉伸把乌兰察布市dem.tif直接拖进 QGIS默认用灰度渲染地形起伏几乎看不出来。正确做法是在图层属性的“符号化”里选择“单波段假彩色”配色方案选terrain或RdYlGr然后在“具色渐变”的设置里把最小/最大值裁剪到 2% 到 98%。这样做的原因是 DEM 数据里常存在极端的局部空洞像元比如噪点或云影残留它们会拉宽直方图导致主要地貌的色带被压缩在中间段。ArcMap 则是在图层属性“符号化”里选择“拉伸”拉伸类型用“百分比截断”截断量填 1~2%色带选择自己习惯的。如果你用 QGIS 的 Python 控制台可以用一行脚本快速加载并设置渲染iface.addRasterLayer(乌兰察布市dem.tif, dem)但注意这行代码只负责加载图层并没有设置拉伸参数。更完整的方式是通过QgsSingleBandPseudoColorRenderer创建渲染器但日常使用中手动点几下比写脚本更快除非你有超过几十个图层需要批量处理。软件推荐拉伸方式参数参考适用场景QGIS百分比截断截断量 2% 到 98%快速浏览大范围地形ArcMap标准差拉伸2 到 3 倍标准差保留局部起伏细节命令行 gdaldem颜色映射自定义色带文件出图用不受软件限制3.2 先对齐投影再做掩膜不要依赖动态投影很多文档会教你在 ArcToolbox 里用“按掩膜提取”直接裁剪但那个工具对输入之间的空间参考一致性要求极高。如果你不先对齐工具大概率会弹出“几何坐标不一致”的警告然后给你一个空结果或错位结果。我的工作流里第一步永远是ogr2ogr把边界 shp 转换为与 tif 相同的投影gdalinfo -proj4 乌兰察布市dem.tif这条命令会把投影参数转成 proj4 形式然后取它的 EPSG 编号比如 CGCS2000 UTM 49N。接下来ogr2ogr -t_srs EPSG:xxxx -overwrite 乌兰察布市范围_proj.shp 乌兰察布市范围.shp参数说明-t_srs指定目标投影-overwrite允许覆盖中间文件。转换后最好再跑一次ogrinfo -so -al确认 Extent 已经变化到新坐标系下这样gdalwarp裁剪时不会因为边界和 DEM 范围差异而出现边缘黑带。3.3 用 gdalwarp 一次性完成重投影与裁剪裁剪 DEM 常见做法是执行gdalwarp它同时内置了重投影能力所以投影对齐和裁剪可以一步完成不必像文档流程那样分两步走。我经常对乌兰察布这种市级数据做“投影转换 边界裁剪”的复合操作gdalwarp -t_srs EPSG:4326 \ -cutline 乌兰察布市范围_proj.shp \ -crop_to_cutline \ -dstalpha \ -tr 0.0003 0.0003 \ -r cubic \ 乌兰察布市dem.tif 乌兰察布市dem_4326_masked.tif-t_srs强制输出到 WGS84 经纬度-cutline指定边界矢量-crop_to_cutline让输出范围严格贴合边界外接矩形顺便把边界外的像元全部裁掉-dstalpha增加一个 alpha 通道边界外透明显示-tr是目标像元尺寸0.0003 度约等于 33 米如果仍想要 30 米可以算准确值或使用-te指定输出范围-r cubic用三次卷积重采样在一定程度上平滑锯齿适合 DEM。如果你要保留原始高程统计特性也可以改用-r bilinear。注意-crop_to_cutline裁剪的是整个边界外接矩形而不是按边界形状做非矩形裁剪。边界外区域虽然保留但会被 alpha 通道遮罩。要得到完全按面要素裁剪的栅格后续还需要gdal_rasterize做掩膜。4. 坡度、坡向与地形起伏度Z因子和重采样方法的选择4.1 填洼不是DEM分析的第一站我经常看到有人一上来就做 Fill Sinks认为 DEM 必须一律先填洼。这是个误解。填洼的目的是消除伪洼地从而获得连续水流方向只有水文分析才需要这一步。如果只是做坡度、坡向、地形起伏度或视域分析填洼反而会抹掉真实的地貌细节比如人工梯田、天然坑塘、采石场坑底。特别是在乌兰察布这类草原与农牧交错区微地形对地表径流和植被分布的影响很明显先算坡度坡向更合理。4.2 gdaldem 计算坡度-p 和 -s 参数的正确姿势GDAL 自带gdaldem工具一行就能生成坡度图gdaldem slope 乌兰察布市dem.tif 乌兰察布市slope.tif -p -s 111120-p控制输出格式加上它输出百分数坡度rise over run最高可超几百不加则输出角度制坡度0–90 度。-s是水平距离的缩放因子直接决定结果是否准确它只在输入为地理坐标系度时才有意义。因为 tif 内部高程单位是米而水平单位是度不设置-sGDAL 会默认水平单位等于高程单位导致坡度几乎全部接近 90 度。这里的111120表示 1 度约等于 111120 米这是全球平均估算值。如果你想更精确就要按当地纬度余弦校正乌兰察布大约北纬 41~43 度cos(42°)≈0.743所以更贴近实际的缩放因子是 111120×0.743≈82562。不过要注意这样做会引入纬度相关的形变不推荐在对精度有要求的科研项目里用。最稳妥的做法是先把 DEM 投影到等面积或 UTM 投影然后让水平和垂直单位都为米此时不需要-s。参数作用常见坑-p输出百分数坡度不加则输出角度值-s地理坐标系下的缩放因子忘记加或加错会得到离谱坡度-r重采样方法默认 nearest 在 DEM 上不推荐-z垂直夸张因子仅用于视觉效果4.3 用 Python rasterio 批量统计坡度分区面积用 gdaldem 生成坡度图只是第一步分析时经常需要知道不同坡度级别的面积占比。用 QGIS 的“栅格计算器”也能做但脚本化更适合多期数据对比。下面这段代码用 rasterio 读取坡度图按国标坡度分级统计面积import rasterio import numpy as np import pandas as pd bins [-1, 5, 15, 25, 35, 90] labels [0-5°, 5-15°, 15-25°, 25-35°, 35°] with rasterio.open(乌兰察布市slope.tif) as src: data src.read(1) nodata src.nodata data np.where(data nodata, np.nan, data) cell_area abs(src.transform[0] * src.transform[4]) flat data.flatten() flat flat[~np.isnan(flat)] counts, _ np.histogram(flat, binsbins) df pd.DataFrame({ grade: labels, area_km2: np.round(counts * cell_area / 1e6, 2) }) df.to_csv(slope_report.csv, indexFalse) print(df)参数说明src.transform[0]是像元宽度src.transform[4]通常为负值所以取绝对值得到面积。如果之前做过投影转换这段脚本能直接吃到 30 米或 33 米的像元面积。分箱边界bins可根据自己的目标改成 0-3°、3-8° 等。输出 CSV 可以直接用 Excel 或后续 Python 画图。5. 一条命令与一个瘦脚本把30米DEM变成可交付的裁剪结果5.1 从 shp 到 KML边界预览与检查在做数据交付时对方往往不想打开 QGIS而是想快速在网页或谷歌地图里看一眼边界。最轻量的是从 Shapefile 转换成 KMLogr2ogr -f KML 乌兰察布市范围.kml 乌兰察布市范围.shp如果你已经做了投影转换用_proj.shp转如果没有转换KML 里会保留原始投影信息但 Google Earth 会按 WGS84 勉强显示可能出现边缘不完全贴合。建议先转成 EPSG:4326 再输出。5.2 用 gdal.Warp 批量裁剪并自动识别投影一致性当你手里有多个 DEM 文件需要按同一边界裁剪时GIS 菜单反复点击的效率太低。我一般用一段 Python 脚本直接从 SDL 调用 GDAL 的WarpAPI。这样既避开了命令行拼接文件名的坑也能在脚本里加投影一致性检查。from osgeo import gdal, ogr import glob cutline 乌兰察布市范围_proj.shp tif_list glob.glob(*.tif) # 先读取边界投影 ref ogr.Open(cutline) ref_lsr ref.GetLayer(0).GetSpatialRef() for tif in tif_list: ds gdal.Open(tif) src_sr ds.GetProjection() if not src_sr and ref_lsr is None: print(f{tif}: 缺少投影信息跳过) out tif.replace(.tif, _masked.tif) gdal.Warp(out, tif, cutlineDSNamecutline, cropToCutlineTrue, dstNodata-9999, options[-overwrite]) print(f{tif} - {out})这里gdal.Warp的第一个参数是输出路径第二个是输入cutlineDSName对应命令行-cutlinecropToCutline对应-crop_to_cutlinedstNodata给输出栅格设定无效值。options[-overwrite]允许覆盖已有文件。代码里的投影检查只是最基础一层实际还可以加一句“如果输入与边界的 EPSG 不一致就先用gdal.Warp做一次srcSRS到边界投影的转换”但更干净的做法是在前面章节里先把 tif 统一成 CGCS2000 / 3-degree Gauss-Kruger再跑这个批量裁剪。批量处理完后用gdalinfo -stats抽查一个新输出文件确认Minimum不等于-9999否则说明边界内部可能还有无效像元残留需要补做填补或参考相邻像元处理。本文还有配套的精品资源点击获取