企业官网建设流程全解析

简介：内蒙古乌兰察布市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单个文件从压缩包拖出来的情况，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.40282e+38。把 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)], crop=True, filled=True, nodata=-32768.0, all_touched=True ) 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函数名字简单，但参数做好合适不太容易。常用参数如下：

all_touched=True对河流、山谷这类线状地貌更友好，能避免边缘像元被过度剥离；但如果后续要统计面积和体积，all_touched=False会让结果更接近，因为只有像元中心落在多边形内才保留。常用于工程量算时使用后者。

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.5～4 才能看出微地形，山区用 1 即可；-az 315是光源方位角，-alt 45是光源高度角。

4.2 用 SHP 点提取高程，以及 shp 转 txt 的最简路径

如果交付对象不是 GIS 用户，而是一个需要把高程写回属性表的分析团队，常见做法是把点状 SHP 与 DEM 叠加采样，再导出 txt：

import 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, masked=True)) points_gdf["dem_value"] = [v[0] for v in vals] points_gdf.drop(columns="geometry").to_csv( "points_dem.txt", sep="\t", index=False, encoding="utf-8-sig" )

src.sample接受坐标对列表，返回每个位置的像元值，masked=True会把 NoData 位置标记为掩码，而不是输出一个刺眼的负数。to_csv(sep="\t")输出制表符分隔文件，正好对应各类 shp 转 txt 的需求；加入encoding="utf-8-sig"是为了让 Excel 打开时中文属性不乱码。

如果要把整个栅格导出成 xyz 文本，使用 gdal_translate：

gdal_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(geometry=grid, crs=boundary.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对交付包重新做一次完整性测试，确保别人拿到手解压后不会看到“文件已损坏”的提示。

本文还有配套的精品资源，点击获取