简介内蒙古乌兰察布市30米分辨率DEM数字高程数据包面向GIS开发者、测绘工程人员、城乡规划研究者及高校地学专业学生提供覆盖乌兰察布市行政边界的高精度地形基础数据可直接支撑区域地形分析、坡度提取、洪涝模拟与可视域计算等专业工作。压缩包共12个文件、约220.2MB核心为tif格式高程栅格另含shp市级范围矢量、dbf属性表、prj坐标系统、tfw地理配准信息、ovr金字塔优化文件及元数据xml等已形成一套规范完整的地理空间数据组合解压后即可在ArcGIS或QGIS中直接加载使用。30米分辨率意味着每个像元对应地面30米×30米的区域数据范围按市级行政边界框定可能包含周边少量过渡地带便于边界分析。已有253人学习下载。利用该tif图层可计算坡度坡向、提取等高线、开展流域划分shp边界文件则用于裁剪研究区、叠加统计从而为城市规划、生态环境评价、土地利用调查及灾害风险管控提供可靠的地形数据支撑。1. 拿到 zip 先别急着加载市级 DEM 交付物要过三道关压缩包文件名里写“30m”不代表打开栅格后像元尺寸一定就是 30m。很多 30m 级别 DEM 的原始文件以度为网格单位例如 SRTM 1 弧秒在乌兰察布这样的中高纬度地区一个像元的东西向实际距离会明显小于 30m南北向则接近 30m。真正能做坡度、坡向和断面分析的应该是经过投影、像元落到 30m 的平面网格。文件名里的“含本市级范围 shp 文件”说明数据包预置了裁剪用的行政边界但这只是素材不是结果。它解决的是在指定行政区范围内做高程分析的标准场景先拿到面状边界再用边界约束 DEM。下面这套流程覆盖了拆包、元数据识别、按 SHP 裁剪、派生地形产品以及交付前的完整性验证既能跑通一次也能写成脚本反复复用。2. 解压 zip 后先查 DEM 元数据SHP、投影和栅格三者的关系打开压缩包之前我最常做的是检查 zip 是否完整确认 SHP 是否带投影信息再用 gdalinfo 看 DEM 的真实网格尺寸。这三步不做完就扔进 ArcMap后面大概率会得到一张错位或全黑的图。2.1 用 unzip -t 检查数据完整性先别急着解压对于这种带中文长文件名的 zip直接在文件管理器里双击解压一旦网络传输中断能解压但栅格中间会出现空洞。速度快、可复现的方式是命令行测试cd ~/data/ulaanqab unzip -t 内蒙古乌兰察布市DEM数字高程数据30m含本市级范围shp文件.zip-t参数是 test 模式只校验每个条目的 CRC 和文件结构不解压文件。输出全部为 “No errors detected” 之后继续如果看到bad CRC或truncated file最省事的处理是重新下载不要尝试修复。数据源通常还会提供 md5 或 sha1下载时保存一份校验文件随后用md5sum -c核对md5sum -c checksum.md5匹配后再解压。这一步防止的是“处理到一半才发现 DEM 负值区域异常”这类隐性损坏。2.2 一个完整 SHP 文件不只是 .shp文件名里写“含本市级范围 shp 文件”但 SHP 实际上是一组文件的集合。常见分发 zip 里至少应该有下面这几个成员后缀用途缺了会怎样.shp几何要素没有几何.shx几何索引部分软件能打开但响应极慢.dbf属性表字段全部丢失.prj投影描述GIS 按默认坐标系读取大概率错位.cpg属性字符编码中文地名乱码.sbn/.sbx空间索引可选ArcGIS 有时会提示重建我曾经遇到过只把.shp单个文件从压缩包拖出来的情况ArcGIS 直接报类似failed to copy spatial iop zip的错误。这不是数据坏了而是软件试图从 zip 内访问.prj却被中间层拦截。常见做法是先把整个 zip 解压到本地工作目录再用ogrinfo检查ogrinfo -so -al 本市范围.shp | grep -E Extent|Layer SRS-so表示只输出摘要-al表示读取全部图层。重点看 Extent 的坐标量级如果范围是 1100000 到 1140000 这种七位数它多半是投影坐标系如果只有 110 到 114则是经纬度。这个判断直接影响下一步是否要把 SHP 转换到 DEM 的坐标。2.3 用 gdalinfo 确认 30m 像元尺寸和 NoData对 DEM 栅格执行gdalinfo dem_30m_ulaanqab.tif | sed -n 1,30p重点看两类输出Pixel Size (30.000000, -30.000000)表示 x、y 方向都是 30m这是已经投影好的等距网格栅格文件自带平面坐标参考。Pixel Size (0.00027777778, 0.00027777778)表示原始数据是 1 弧秒地理网格并不是严格意义上的 30m 像元。要进入坡度分析必须先投影并重采样。同一行附近还有Coordinate System is和NoData Value。NoData 常见值是-9999、-32768或-3.40282e38。把 NoData 记下来后续 mask 裁剪时如果不显式传入输出文件可能沿用源文件的 NoData也可能被默认值覆盖导致整片裁出区域显示为空值。3. 用 Python 裁剪 DEM 数据基于本市级范围 SHP 的掩膜方案QGIS 里有 Raster Extraction 工具手动点一次没问题但当你需要换边界、换季节、换分辨率重跑时脚本是更能复用的方案。这里用 geopandas 读取 SHP用 rasterio 的 mask 函数完成按面裁剪。3.1 先只保留外边界线并统一 CRS再进行掩膜拿到“本市范围.shp”后第一步不是直接 mask而是把面状数据统一到 DEM 的坐标系。DEM 的 CRS 可以通过rasterio.open(dem_path).crs读出来不需要手动指定 EPSG。另一方省略这个步骤直接用原始经纬度的 SHP 去裁剪 UTM 的 DEM那输出的范围要么为空要么产生严重形变。import geopandas as gpd import rasterio from rasterio.mask import mask from shapely.geometry import mapping shp_path 本市范围.shp dem_path dem_30m_ulaanqab.tif region_gdf gpd.read_file(shp_path) with rasterio.open(dem_path) as src: region_gdf region_gdf.to_crs(src.crs) # 用融合后的多边形做裁剪避免多部件面的边界重叠产生细缝 combined region_gdf.dissolve() geom combined.geometry.iloc[0] out_img, out_transform mask( src, [mapping(geom)], cropTrue, filledTrue, nodata-32768.0, all_touchedTrue ) out_meta src.meta.copy() out_meta.update({ height: out_img.shape[1], width: out_img.shape[2], transform: out_transform, nodata: -32768.0 }) with rasterio.open(ulaanqab_dem_clip.tif, w, **out_meta) as dst: dst.write(out_img) # 如果需要边界线可在同一对象上生成 LineString gpd.GeoSeries(geom.boundary).to_file(边界线.shp)dissolve()会把该市所有区县面要素合并成一个要素消除邻接面内部边界。这里的逻辑是裁剪只关心全域覆盖不关心内部的行政区划线若把区县面全部直接传给 mask重叠部分会重复计算边界衔接处还可能产生 1 个像元的缝隙。若这个 SHP 本来就是单一市域面dissolve()不会改变结果只是更保险。mapping(geom)把 shapely 几何转成 rasterio mask 函数需要的 GeoJSON 字典格式。不能直接把geom作为 shapes 传入某些旧版本会报元数据序列化错误。3.2 rasterio.mask.mask 参数表和选择依据mask函数名字简单但参数做好合适不太容易。常用参数如下参数建议值作用cropTrue把输出窗口裁剪到 SHP 边界去掉外部大范围冗余像元filledTrue边界外填充为 NoDataFalse 时返回 numpy 掩码数组nodata-32768.0显式指定 NoData避免沿用源文件中的不合理值all_touchedTrue表示与边界相切的像元也保留减少边缘锯齿invertFalse反向裁剪只保留边界以外区域一般用不到all_touchedTrue对河流、山谷这类线状地貌更友好能避免边缘像元被过度剥离但如果后续要统计面积和体积all_touchedFalse会让结果更接近因为只有像元中心落在多边形内才保留。常用于工程量算时使用后者。3.3 裁剪后范围对不上时的三个排查点如果拿到ulaanqab_dem_clip.tif后发现范围明显偏斜或整幅图全是 NoData按下面顺序排查with rasterio.open(dem_path) as src: print(DEM 范围:, src.bounds) print(SHP 范围:, region_gdf.to_crs(src.crs).total_bounds)先看 DEM 范围和 SHP 范围是否在同一个数量级。如果 SHP 范围是110.5, 40.5, 114.5, 43.0而 DEM 范围是410000, 4400000说明to_crs没有生效检查是否在with块外重新对变量赋值。再看 NoData。裁剪结果边缘的黑边本来正常但如果面积的一半以上是黑检查传入的nodata是否与原数据一致。源文件用-9999你传-32768边界外会被填成-32768而源数据里的真实空洞仍是-9999两层 NoData 混在一起。最后看投影。乌兰察布市跨 UTM 49N 的东半区如果 DEM 投影是 CGCS2000 而 SHP 是 WGS84两者的差值可能只有几十米视觉上能对上但坡度计算时会出现边缘错位。这时不以文件名后缀做判断统一以src.crs为准。4. 从 30m DEM 生成坡度、山体阴影与 SHP 转 TXT 提取值裁剪只是把范围限制住了真正要交付的是坡度、山体阴影或者点值表。这里的方向是先在投影坐标系下重采样到 30m再派生产品不要拿一个度坐标系 DEM 直接跑坡度命令。4.1 用 gdalwarp 重投影并用 gdaldem 生成坡度如果第 2 步发现 DEM 仍是经纬度网格要先重投影。乌兰察布市的经度范围大致落在 UTM 49N可以这样对待gdalwarp -t_srs EPSG:32649 -tr 30 30 -r bilinear \ -cutline 本市范围.shp -crop_to_cutline \ dem_30m_ulaanqab.tif ulaanqab_utm.tif-t_srs EPSG:32649指定 WGS 84 / UTM 49N市区高程分析常用。-tr 30 30强制输出像元为 30m但要注意这与原始网格的对齐方式有关如果原始像元是 0.000277 度重采样后边缘会产生插值误差。-r bilinear适合连续高程表面不要用 nearest否则等高线会出现块状台阶。重投影后再处理坡度和山体阴影gdaldem slope ulaanqab_utm.tif slope_deg.tif -p gdaldem hillshade ulaanqab_utm.tif hillshade.tif -z 2.5 -az 315 -alt 45-p表示坡度用度输出不加时是百分比。山体阴影的-z是垂直拉伸因子平原地区用到 2.54 才能看出微地形山区用 1 即可-az 315是光源方位角-alt 45是光源高度角。4.2 用 SHP 点提取高程以及 shp 转 txt 的最简路径如果交付对象不是 GIS 用户而是一个需要把高程写回属性表的分析团队常见做法是把点状 SHP 与 DEM 叠加采样再导出 txtimport geopandas as gpd import rasterio points_gdf gpd.read_file(采样点.shp) with rasterio.open(ulaanqab_utm.tif) as src: points_gdf points_gdf.to_crs(src.crs) coords [(p.x, p.y) for p in points_gdf.geometry] vals list(src.sample(coords, maskedTrue)) points_gdf[dem_value] [v[0] for v in vals] points_gdf.drop(columnsgeometry).to_csv( points_dem.txt, sep\t, indexFalse, encodingutf-8-sig )src.sample接受坐标对列表返回每个位置的像元值maskedTrue会把 NoData 位置标记为掩码而不是输出一个刺眼的负数。to_csv(sep\t)输出制表符分隔文件正好对应各类 shp 转 txt 的需求加入encodingutf-8-sig是为了让 Excel 打开时中文属性不乱码。如果要把整个栅格导出成 xyz 文本使用 gdal_translategdal_translate -of XYZ -srcnodata -32768 -dstnodata NaN \ ulaanqab_utm.tif ulaanqab_dem.xyz-of XYZ输出三列文本X、Y、Z。大范围 DEM 的 XYZ 文件可能上千万行直接用文本编辑器打开会很吃力后续一般用awk或pandas按边界过滤。4.3 DSM 与 DEM 的区别以及 12.5m 数据混用的边界在 OpenTopography 下载 DEM 教程里你会发现同一地区有 DEM 和 DSM 两种产品。DEM 是地表裸面模型去掉了树冠和建筑顶DSM 是表面模型包含屋顶和树冠。ALOS 12.5m 这类数据在官方说明里也明确部分是 DSM不是所有“高程数据”都能当 DEM 用。把 12.5m DSM 重采样成 30m 网格并叠加到现有数据时不要直接用gdalwarp默认 nearest 后再镶嵌否则建筑物屋顶会以离散高亮点形式出现在地形渲染中。先做重采样并对异常高差做低通滤波更稳妥gdalwarp -tr 30 30 -r cubic alos_12_5_dsm.tif alos_12_5_30m.tif-r cubic是平滑插值能抑制边缘锯齿但也会让建筑边界模糊。把重采样后的数据与 30m DEM 做差值差值大于 5m 的位置大概率是建筑或树木需要从最终地形产品中剔除。5. 交付前用渔网分割 SHP 检查空洞再压缩 zip 并验证到这里数据已经裁剪并派生完成但直接压缩交付风险仍然存在比如 DEM 某个瓦片缺失或 SHP 的投影信息在压缩时被漏掉。可以用渔网把整个市域划分成规则块快速定位空洞再做完整压缩检查。5.1 用渔网分割 SHP快速发现 DEM 空洞区域渔网分割 SHP 是处理大面积栅格时常用的检查手段。把本市边界按 5km 边长切分成方格每个格子统计有无有效像元。import geopandas as gpd from shapely.geometry import box boundary gpd.read_file(边界线.shp) xmin, ymin, xmax, ymax boundary.total_bounds size 5000 # 按投影坐标单位UTM 下即 5km grid [] cols int((xmax - xmin) // size) 1 rows int((ymax - ymin) // size) 1 for r in range(rows): for c in range(cols): grid.append(box( xmin c * size, ymin r * size, min(xmin (c 1) * size, xmax), min(ymin (r 1) * size, ymax) )) grid_gdf gpd.GeoDataFrame(geometrygrid, crsboundary.crs) grid_gdf.to_file(fishnet_5km.shp)这个渔网只做辅助检查不上生产。用它和裁剪后的 DEM 做叠加任何一个格子如果没有高程数据就能快速定位该格子的经纬度范围再去查原始数据是缺失还是被 NoData 掩盖。5.2 只保留外边界线并完成 zip 完整性检查如果需要向外输出边界线而不是整个面状范围用几何对象的boundary属性生成线boundary_line boundary.geometry.boundary boundary_line.to_file(外边界线.shp)提交前把派生文件统一压缩ogr2ogr -f KML 外边界线.kml 外边界线.shp zip -r final_dem_package.zip *.tif *.shp* README.txt unzip -t final_dem_package.zip压缩时记得把.prj和.cpg放进包里*.shp*会覆盖.shp/.shx/.dbf/.prj/.cpg这些主要附属文件。README 至少写三行投影 EPSG、NoData 值、像元尺寸。最后用unzip -t对交付包重新做一次完整性测试确保别人拿到手解压后不会看到“文件已损坏”的提示。本文还有配套的精品资源点击获取