企业官网建设流程全解析

简介：这份云南省30米分辨率DEM数据面向地理信息、测绘、环境研究与城市规划等领域的从业者和学习者，提供覆盖全省的高程地形基础数据，可用于洪水风险评估、地质灾害分析、地形地貌研究及交通线路设计等场景。压缩包共10个文件，约1MB，以GeoTIFF栅格影像为核心，辅以shp、shx、dbf、prj等矢量边界文件，以及tfw、xml、sbn、sbx等坐标与索引辅助文件，采用WGS84坐标系，便于在QGIS、ArcGIS等软件中直接加载与空间分析。目前已有1057人学习下载。数据基于ASTER GDEM V3版本，精度为30米，读者可据此完成云南省地形可视化、坡度坡向提取、流域分析等操作，并借助配套边界文件快速裁剪研究区，是开展区域地理研究、课程实验与项目建模的实用底图资料。

1. 云南省DEM（30米分辨率）：从数据认知到地形分析的落地路径

拿到一份省级30米DEM，很多人第一反应是“不就是个高程栅格吗”，结果一打开ArcGIS就卡在投影不对、范围对不上、nodata值离谱这些破事上。云南省DEM（30米分辨率）覆盖全省39.4万平方公里，从滇西北海拔6740米的卡瓦格博峰到河口县76米的河谷，高差超过6600米，这种极端地形对数据精度和后续分析流程的要求比平原地区苛刻得多。这份资源适合做水文分析、坡度坡向提取、地形起伏度计算、工程选址、遥感影像正射校正的从业者，也适合需要省级尺度地形底图做空间建模的研究人员。它解决的核心问题是：给你一套可以直接进入GIS工作流的高程数据，省掉从公开源拼接、裁剪、重投影的重复劳动。但能不能用好，取决于你对投影、分辨率、nodata和地形分析参数的理解深度。

2. 30米DEM的技术底座：投影、基准与分辨率选型

2.1 为什么是30米，而不是12.5米或90米

30米分辨率对应的是SRTM（航天飞机雷达地形测绘任务）和ASTER GDEM的经典格网间距，1弧秒约等于30米。这个尺度在省级尺度上是一个平衡点：90米太粗，做小流域水文分析时河道会断线；12.5米更精细，但云南省面积大，全量数据量会膨胀到几十GB，普通工作站处理起来内存吃紧。30米在云南省范围内，单波段16位整型，压缩后大约几百MB到1GB出头，笔记本16GB内存能扛住分块处理。常见做法是：省级宏观分析用30米，重点区域再用更高分辨率数据补充。

选型时还要注意一个坑：不同来源的30米DEM垂直精度差异很大。SRTM v3的绝对高程精度约±16米，ASTER GDEM v3约±17米，但在地形陡峭区误差会放大。云南省滇西北高山峡谷区，DEM高程误差在坡度大于30度的区域可能达到20米以上，做坡度计算时这个误差会被进一步放大。如果你的应用对高程绝对值敏感，比如大坝选址，30米DEM只能做初筛，不能做最终设计依据。

2.2 投影与基准：WGS84地理坐标还是UTM投影

拿到数据第一件事是确认坐标系。云南省DEM常见分发格式是WGS84地理坐标系（EPSG:4326），单位是度，像元大小0.000277778度。这种格式适合存储和分发，但不适合做面积、距离、坡度计算，因为经纬度格网在高纬度地区不是正方形。

做地形分析前，我一般会先投影到适合云南的投影坐标系。云南省跨UTM 47N和48N两个带，中央经线分别是99°E和105°E。更推荐用Albers等面积投影，自定义参数如下：

# GDAL 投影转换：WGS84 地理坐标 -> Albers 等面积投影（云南适用） gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=29 +lat_0=27 +lon_0=102 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" \ -tr 30 30 \ -r bilinear \ -of GTiff \ -co COMPRESS=LZW \ -co TILED=YES \ yunnan_dem_wgs84.tif \ yunnan_dem_albers.tif

这段命令做了三件事：把地理坐标转成Albers等面积投影，重采样到30米×30米像元，用LZW压缩并分块存储。-tr 30 30指定输出分辨率，-r bilinear是双线性插值，适合连续型高程数据。-co TILED=YES让后续按块读取更快，不然GDAL每次都要读整行。注意：如果原始数据已经是投影坐标，不要再转一次，否则重采样会引入额外误差。

2.3 nodata与无效值处理

云南省DEM的nodata值常见是-9999或-32768。如果直接拿去做坡度计算，这些值会被当成真实高程参与运算，结果就是边界出现巨大异常值。处理方式是在投影转换时同时设nodata：

gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=29 +lat_0=27 +lon_0=102 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" \ -tr 30 30 \ -r bilinear \ -srcnodata -9999 \ -dstnodata -9999 \ -of GTiff \ -co COMPRESS=LZW \ yunnan_dem_wgs84.tif \ yunnan_dem_albers.tif

-srcnodata告诉GDAL原始数据里哪个值是无效的，-dstnodata指定输出无效值。如果原始nodata是-32768，就改成-32768。处理完之后用gdalinfo -stats看一眼统计值，确认最小值不是-9999，否则说明nodata没设对。

3. 从DEM到地形因子：坡度、坡向与起伏度提取实操

3.1 坡度与坡向：算法选择与单位陷阱

坡度计算有两种常用算法：Horn算法和Zevenbergen-Thorne算法。GDAL默认用Horn，ArcGIS默认用Zevenbergen-Thorne。两者在平缓地区差异不大，但在陡峭地形区，Zevenbergen-Thorne对噪声更敏感。云南省高山峡谷区，我一般用Horn，结果更平滑。

用GDAL算坡度：

# 坡度计算，输出单位为度 gdaldem slope \ -alg Horn \ -compute_edges \ -of GTiff \ -co COMPRESS=LZW \ yunnan_dem_albers.tif \ yunnan_slope.tif # 坡向计算，输出0-360度，正北为0 gdaldem aspect \ -compute_edges \ -of GTiff \ -co COMPRESS=LZW \ yunnan_dem_albers.tif \ yunnan_aspect.tif

-alg Horn指定算法，-compute_edges让边缘像元也能算出值，不然边界一圈是nodata。坡度输出默认是度，如果你要弧度，加-p。坡向输出是0到360度，0是正北，90是正东。注意：坡向在平坦区域（坡度接近0）没有意义，算出来是随机值，后续分析时要按坡度阈值掩膜掉。

3.2 地形起伏度：窗口大小的选择逻辑

地形起伏度是指定窗口内最大高程与最小高程之差，反映地形破碎程度。窗口大小直接决定结果：窗口太小，结果噪声大；窗口太大，起伏度被平滑，失去局部特征。云南省地形复杂，我一般用3×3窗口做微观起伏，用11×11窗口做中观起伏，用21×21窗口做宏观起伏。

用Python和GDAL计算起伏度：

from osgeo import gdal import numpy as np from scipy.ndimage import maximum_filter, minimum_filter # 打开DEM ds = gdal.Open("yunnan_dem_albers.tif") band = ds.GetRasterBand(1) dem = band.ReadAsArray().astype(np.float32) nodata = band.GetNoDataValue() # 把nodata设为nan，避免参与计算 dem[dem == nodata] = np.nan # 定义窗口大小 window_size = 11 # 计算最大值和最小值滤波 max_dem = maximum_filter(dem, size=window_size) min_dem = minimum_filter(dem, size=window_size) # 起伏度 = 最大值 - 最小值 relief = max_dem - min_dem # 恢复nodata relief[np.isnan(relief)] = nodata # 写出结果 driver = gdal.GetDriverByName("GTiff") out_ds = driver.Create("yunnan_relief_11x11.tif", ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_band = out_ds.GetRasterBand(1) out_band.WriteArray(relief) out_band.SetNoDataValue(nodata) out_band.FlushCache() out_ds = None

这段代码用scipy.ndimage的maximum_filter和minimum_filter做窗口统计，比手写循环快几个数量级。window_size=11对应11×11像元，在30米分辨率下约等于330米×330米的窗口。如果你要算21×21，改成21就行。注意：maximum_filter和minimum_filter对nan的处理是传播nan，所以窗口内只要有一个nan，输出就是nan，这正好符合nodata不参与计算的需求。

3.3 水文分析：填洼与流向计算

DEM做水文分析前必须填洼，否则水流方向会被局部凹陷打断。GDAL没有内置填洼，常用的是richdem或WhiteboxTools。用richdem：

import richdem as rd # 读取DEM dem = rd.LoadGDAL("yunnan_dem_albers.tif") # 填洼 dem_filled = rd.FillDepressions(dem, epsilon=True, in_place=False) # 计算流向（D8算法） flow_dir = rd.FlowDirections(dem_filled, method='D8') # 保存结果 rd.SaveGDAL("yunnan_dem_filled.tif", dem_filled) rd.SaveGDAL("yunnan_flow_dir.tif", flow_dir)

epsilon=True表示填洼时保留微小坡度，避免大面积平坦区导致流向不确定。method='D8'是八方向流向算法，每个像元流向八个邻居中坡度最陡的那个。填洼后的DEM再算流向，河网才连续。注意：richdem对内存要求较高，省级DEM建议分块处理，或者用WhiteboxTools的FillDepressions，它支持分块。

4. 避坑与排查：30米DEM处理中的五个血泪教训

4.1 现象：坡度计算结果出现大面积0值或异常高值

原因：nodata值没设对，或者投影转换时重采样方法选错。如果原始nodata是-32768，你没设，GDAL会把它当真实高程，算出来的坡度要么是0（因为-32768周围都是-32768），要么是巨大值（因为-32768和真实高程之间高差几千米）。

解决：用gdalinfo -stats看原始数据的统计值，确认nodata。然后在gdalwarp和gdaldem里都显式指定-srcnodata和-dstnodata。如果已经算错了，重新从投影转换那一步开始，不要试图在坡度结果上修补。

4.2 现象：投影转换后范围偏移，和矢量边界对不上

原因：原始数据的坐标系定义缺失或错误。有些DEM分发时没有嵌入投影信息，GDAL默认按WGS84地理坐标处理，但实际可能是其他基准。

解决：先用gdalinfo看有没有Coordinate System is这一行。如果没有，用gdal_edit.py -a_srs EPSG:4326补上。如果补上后还对不上，用gdalsrsinfo对比矢量边界的坐标系，确认基准是否一致。云南省常见的是WGS84和CGCS2000，两者在30米尺度上差异很小，但严格来说不能混用。

4.3 现象：填洼后DEM出现大面积平坦区，流向计算失败

原因：填洼算法把大片区域填成同一高程，D8算法在平坦区无法确定流向。

解决：用richdem的epsilon=True选项，填洼时保留微小坡度。如果已经填平了，用rd.FlowDirections的method='D8'配合rd.ResolveFlats处理平坦区。或者换用WhiteboxTools的FillDepressions，它默认带坡度保持。

4.4 现象：起伏度计算结果在边界出现异常值

原因：窗口滤波在边界处窗口不完整，maximum_filter和minimum_filter默认用reflect模式填充边界，导致边界值失真。

解决：在计算起伏度前，先把DEM边缘裁剪掉窗口半径的宽度，或者用mode='constant'并指定cval=np.nan，让边界输出nan。我一般用后者，然后在后处理时把边界nan掩膜掉。

4.5 现象：省级DEM处理时内存溢出

原因：30米云南省DEM全量读入内存约几GB，加上中间结果，16GB内存容易爆。

解决：用GDAL的分块读取，每次处理一个block。或者用rasterio的block_windows：

import rasterio from rasterio.windows import Window with rasterio.open("yunnan_dem_albers.tif") as src: for ji, window in src.block_windows(1): data = src.read(1, window=window) # 处理data # 写出到目标文件对应window

block_windows按GeoTIFF内部的分块大小逐块读取，内存占用可控。注意：分块处理时窗口滤波需要额外读取边界像元，不然块与块之间会出现接缝。常见做法是每个块向外扩窗口半径，处理完再裁掉。

5. 进阶技巧：用DEM做地形阴影与三维可视化验证

5.1 地形阴影：光照方位角与高度角的参数选择

地形阴影（hillshade）是验证DEM质量最直观的方式。GDAL的gdaldem hillshade默认方位角315度、高度角45度，这是制图惯例。但云南省地形复杂，默认参数在南北向河谷会显得过暗。我一般用方位角315度、高度角30度，让阴影更柔和，细节更丰富。

gdaldem hillshade \ -az 315 \ -alt 30 \ -z 2 \ -compute_edges \ -of GTiff \ -co COMPRESS=LZW \ yunnan_dem_albers.tif \ yunnan_hillshade.tif

-z 2是垂直 exaggeration 因子，把高程放大2倍，让地形起伏更明显。-az 315是光照方位角，315度是西北方向，符合北半球制图习惯。-alt 30是光照高度角，30度比45度更能突出微地形。注意：-z不要设太大，超过3会让阴影失真，看起来像褶皱。

5.2 三维可视化：用QGIS和Blender快速出图

QGIS里加载DEM和hillshade，把hillshade放在DEM下面，DEM设半透明，就能做出立体地形图。如果要更高质量的三维渲染，导出DEM为OBJ格式，用Blender加光照和材质。

# DEM 转 OBJ gdal_translate -of OBJ \ -co FACES=YES \ -co FORMAT=OBJ \ yunnan_dem_albers.tif \ yunnan_dem.obj

FACES=YES生成三角面片，FORMAT=OBJ指定输出格式。导出的OBJ可以直接拖进Blender，加一个太阳光，调整角度，渲染出来就是三维地形。注意：省级DEM转OBJ文件会很大，建议先降采样或者裁剪重点区域。

5.3 验证DEM质量的三个快速检查

第一，看hillshade有没有明显接缝或条带，有的话说明数据拼接有问题。第二，算坡度，看坡度大于60度的区域占比，如果超过5%，可能是DEM噪声。第三，和已知高程点对比，云南省内选几个GNSS控制点，看DEM高程和实测高程的差值，30米DEM在平坦区差值应在10米以内，陡峭区20米以内。

从那以后我每次拿到新DEM，都强制走一遍“gdalinfo看元数据 → 投影转换 → nodata检查 → hillshade目视 → 坡度统计”这个流程，不跳过任何一步。希望帮到你。

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