尧图网页设计 
简介这份资源面向地理信息、城乡规划、测绘及环境研究等领域的从业者与学习者提供湖南省湘西土家族苗族自治州30米分辨率的DEM数字高程数据并附带本市级行政范围矢量边界可用于地形分析、坡度计算、洪水模拟与制图等GIS任务。压缩包共12个文件约50.18MB主要包含tif格式的高程栅格数据、shp格式的行政边界矢量文件以及prj坐标系统、tfw坐标配准、dbf属性表、shx与sbn/sbx索引、ovr金字塔缩略图和xml元数据等配套文件构成一套可直接在ArcGIS、QGIS中加载的完整空间数据集。目前已有301人学习下载。数据覆盖范围完整、分辨率适中既能满足区域尺度地形建模需求也便于与行政边界叠加开展空间统计适合作为课程实验、项目底图或规划分析的基础数据使用。1. 湘西州 30m DEM 数据到手后先搞清楚它能干什么、不能干什么如果你正在做湘西土家族苗族自治州范围内的地形分析、水文建模、选址评估或者三维可视化大概率绕不开一份 30m 分辨率的 DEM 数字高程数据。这次拆的这份资源是覆盖湘西州本市级范围的 30m DEM 栅格外加配套的 shp 边界文件。很多人拿到这类压缩包的第一反应是直接拖进 ArcMap 或者 QGIS 里看效果但真正决定后面顺不顺的是你在打开之前有没有想清楚三件事这份数据的精度够不够支撑你的分析尺度、shp 边界和栅格范围对不对得上、以及你打算用它做提取还是做裁剪。30m 分辨率意味着每个像元代表地面约 900 平方米的范围这个尺度适合流域级、县域级的坡度坡向分析、汇流计算和地形起伏度统计但拿它去做地块级的土方量计算或者精细的场地平整设计精度就不够了。湘西州地处武陵山区地形起伏大、沟谷密集30m 数据在描述山脊线和沟谷走向时会有一定的平滑效应这一点在做水文分析时要特别留意。配套的 shp 文件通常是行政边界或者本市级范围面图层它的作用不是装饰而是帮你把分析范围卡死在目标区域内避免跑出边界算出无意义的结果。适合谁用做区域规划、自然资源调查、遥感预处理、教学演示的从业者以及需要快速拿到一份可用地形底图的新手。不适合谁做厘米级工程测量或者城市级精细建模的团队这个精度撑不住。2. 把 DEM 和 shp 正确装进工作流从坐标系检查到范围对齐2.1 先验坐标系别等裁剪完才发现偏了几百米拿到栅格和矢量数据第一步永远不是裁剪而是确认两者的坐标系是否一致。DEM 常见的是地理坐标系 WGS84 或者 CGCS2000单位是度shp 边界可能是投影坐标系单位是米。如果直接拿投影坐标的 shp 去裁地理坐标的栅格ArcMap 会弹警告QGIS 可能直接给你裁出一片空白或者错位的结果。我一般会先在 ArcMap 里右键图层看属性或者用 QGIS 的图层属性面板确认 CRS。如果两者不一致优先把 shp 投影到和 DEM 相同的坐标系而不是反过来动栅格因为重采样栅格会引入额外的高程误差。# 用 rasterio 和 geopandas 快速检查 DEM 与 shp 的坐标系 import rasterio import geopandas as gpd dem_path xiangxi_dem_30m.tif shp_path xiangxi_boundary.shp with rasterio.open(dem_path) as src: print(DEM CRS:, src.crs) print(DEM 范围:, src.bounds) print(DEM 分辨率:, src.res) gdf gpd.read_file(shp_path) print(SHP CRS:, gdf.crs) print(SHP 范围:, gdf.total_bounds)这段代码的逻辑很直接rasterio 读栅格的 CRS、边界和分辨率geopandas 读矢量的 CRS 和总范围。参数上重点看src.res30m 数据在地理坐标系下大约是 0.000277 度左右如果显示的是 30 或者 0.000277 都算正常但单位不同后续处理方式完全不同。total_bounds返回的是 shp 的四至范围拿它和 DEM 的 bounds 对比如果 shp 范围明显大于 DEM说明边界文件可能包含了州外区域裁剪时要以 DEM 实际覆盖为准。2.2 用面图层裁剪 DEMExtract by Mask 和 Clip 的区别要分清这是搜索热词里反复出现的问题ArcMap 中依靠面图层裁剪 DEM 栅格和依靠面图层掩码提取到底有什么区别。简单说Clip工具是几何裁剪输出栅格的像元值不变只是把边界外的像元去掉边界上的像元按面积比例处理Extract by Mask是掩码提取边界外的像元被设为 NoData边界内的像元值原样保留。两者在边界规整、范围完全包含的情况下结果几乎一样但当 shp 边界和栅格像元边界不重合时Clip 可能会在边缘产生半个像元的问题而 Extract by Mask 更干净。# 用 rasterio 的 mask 功能实现按面裁剪 DEM import rasterio from rasterio.mask import mask import geopandas as gpd with rasterio.open(xiangxi_dem_30m.tif) as src: gdf gpd.read_file(xiangxi_boundary.shp) # 确保 shp 和 DEM 坐标系一致 if gdf.crs ! src.crs: gdf gdf.to_crs(src.crs) geoms [geom for geom in gdf.geometry] out_image, out_transform mask(src, geoms, cropTrue) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(xiangxi_dem_clipped.tif, w, **out_meta) as dest: dest.write(out_image)这段代码的关键参数是cropTrue它会把输出栅格的范围收缩到 shp 的实际边界而不是保留原始 DEM 的完整范围。out_meta继承原栅格的元数据后更新了行列数和仿射变换保证输出文件的地理定位正确。如果你在 ArcMap 里操作对应的是 Spatial Analyst 工具条下的 Extract by Mask输入栅格选 DEM掩码选 shp输出就是裁剪后的结果。注意一点如果 shp 有多个面要素mask 函数会取所有面的并集作为掩码范围这和 ArcMap 里逐个要素裁剪的逻辑不同需要分开处理时得先拆分 shp。2.3 裁剪后必做的三项检查范围、NoData 和统计值裁剪完不是就结束了我见过太多人裁完直接拿去算坡度结果发现边缘一圈全是异常值。第一项检查是范围把裁剪后的栅格和 shp 叠在一起看边界应该严丝合缝不能出现栅格超出 shp 或者 shp 内有大片空白的情况。第二项检查是 NoData 值Extract by Mask 会把边界外设为 NoData但有些 DEM 原始数据本身在河流或者建筑区域就有 NoData裁剪后这些值会混在一起需要在属性里确认 NoData 的数值。第三项是统计值用栅格统计工具看最小、最大和均值湘西州的高程范围大致在 200 到 1700 米之间如果最小值出现负数或者最大值超过 2000大概率是数据本身有问题或者裁剪时混入了异常像元。# 裁剪后快速统计高程分布 import rasterio import numpy as np with rasterio.open(xiangxi_dem_clipped.tif) as src: data src.read(1) nodata src.nodata if nodata is not None: valid data[data ! nodata] else: valid data.flatten() print(最小值:, np.min(valid)) print(最大值:, np.max(valid)) print(均值:, np.mean(valid)) print(有效像元数:, valid.size)这段代码先读取第一个波段然后根据 nodata 值过滤掉无效像元再算统计量。参数上注意src.nodata可能是 None这时候需要手动指定一个明显不合理的值来过滤比如 -9999。湘西州的地形统计如果均值在 600 到 800 米之间、标准差在 200 米左右说明数据分布合理如果均值异常低或者标准差极小可能是裁剪时把大片 NoData 当成了有效值。3. 从 DEM 到 shp等高线、坡度分级和流域边界的提取方法3.1 等高线提取参数设错就是一堆碎线从 DEM 提取等高线是常见需求ArcMap 里用 Contour 工具QGIS 里用等值线提取。核心参数是等高距30m 数据建议等高距设在 10 到 20 米之间设太小会生成大量碎线设太大又丢失地形细节。湘西州山地起伏大我一般用 10 米等高距做整体展示局部区域用 5 米。提取出来的等高线 shp 要检查是否闭合、是否有自相交这些在后续做地形分析时都会导致拓扑错误。# 用 GDAL 从 DEM 提取等高线 from osgeo import gdal, ogr dem gdal.Open(xiangxi_dem_clipped.tif) band dem.GetRasterBand(1) driver ogr.GetDriverByName(ESRI Shapefile) out_ds driver.CreateDataSource(contour_10m.shp) out_layer out_ds.CreateLayer(contour, geom_typeogr.wkbLineString) out_layer.CreateField(ogr.FieldDefn(elev, ogr.OFTReal)) gdal.ContourGenerate(band, 10, 0, [], 0, 0, out_layer, 0, 0) out_ds NoneContourGenerate的参数依次是波段、等高距、基准高程、忽略值列表、是否用 NoData、是否用固定间隔、输出图层、字段索引等。这里等高距设为 10基准高程设为 0表示从 0 开始每 10 米生成一条线。生成的 shp 里 elev 字段记录了每条等高线的高程值方便后续按高程筛选或者做标注。3.2 坡度分级和坡向分析别直接拿原始坡度做规划图DEM 算坡度是基础操作但直接输出的坡度栅格是连续值做规划图或者统计报告时需要分级。湘西州的山地坡度分级一般按 0-5 度、5-15 度、15-25 度、25-35 度、35 度以上来分对应平缓、较缓、中等、陡峭和极陡。ArcMap 里用 Reclassify 工具QGIS 里用栅格重分类。坡向分析同理输出的是 0 到 360 度的方向值需要按北、东北、东等八个方向重分类。注意坡度计算前要确认 DEM 的 Z 单位地理坐标系下 Z 单位是米而 XY 单位是度直接算坡度会得到错误结果需要先用投影坐标系或者设置 Z 因子。# 用 GDAL 计算坡度并重分类 from osgeo import gdal import numpy as np dem gdal.Open(xiangxi_dem_clipped.tif) gt dem.GetGeoTransform() # 地理坐标系下需要设置 Z 因子约等于 1/111320 z_factor 1.0 / 111320.0 slope_ds gdal.DEMProcessing(slope.tif, dem, slope, zFactorz_factor) slope_data slope_ds.GetRasterBand(1).ReadAsArray() # 按 5 度间隔分级 bins [0, 5, 15, 25, 35, 90] classified np.digitize(slope_data, bins) print(各级别像元数:, [np.sum(classified i) for i in range(1, len(bins))])DEMProcessing的zFactor参数是关键地理坐标系下不设置这个值坡度会被严重放大。分级用np.digitize按边界值把连续坡度映射到离散类别输出的是 1 到 5 的整数栅格。实际项目中我会把这个分级结果再转成 shp 或者做掩码统计方便出报告。3.3 流域边界和河网提取填洼是绕不过去的一步用 DEM 做水文分析填洼是第一步也是最多人翻车的地方。原始 DEM 里存在洼地不填洼直接算流向会出现断流或者内流区。ArcMap 里用 Fill 工具QGIS 里用 Wang and Liu 填洼算法。填洼之后算流向、流量累积设定阈值提取河网再转成 shp。湘西州的喀斯特地貌区域洼地特别多填洼的阈值要适当调大否则会把大量真实洼地填平改变水文响应特征。# 用 RichDEM 做填洼和流量累积 import richdem as rd dem rd.LoadGDAL(xiangxi_dem_clipped.tif) dem_filled rd.FillDepressions(dem, epsilonTrue, in_placeFalse) flow_accum rd.FlowAccumulation(dem_filled, methodD8) rd.SaveGDAL(flow_accum.tif, flow_accum)FillDepressions的epsilonTrue表示在填洼时加入微小坡度避免平坦区域流向不确定。FlowAccumulation用 D8 算法每个像元的流向指向最陡的下坡方向流量累积值表示上游汇水像元数。阈值一般设在 1000 到 5000 之间具体看流域面积和河网密度湘西州建议从 2000 开始试。4. 避坑与排查30m DEM 处理中最容易翻车的五个地方4.1 裁剪后栅格范围对不上 shp 边界现象裁剪后的 DEM 边缘和 shp 边界有明显偏移或者 shp 内出现大片 NoData。原因通常是坐标系不一致或者 shp 本身有拓扑错误比如自相交、重叠面。解决方法是先统一坐标系再用拓扑检查工具修复 shp最后重新裁剪。如果 shp 是多个面要素确认是否需要合并成一个面再裁剪。4.2 坡度计算结果明显偏大现象算出来的坡度动辄七八十度明显不符合实际地形。原因是在地理坐标系下直接计算坡度没有设置 Z 因子。解决方法是先把 DEM 投影到投影坐标系或者在坡度工具里手动设置 Z 因子为 1/111320 左右。ArcMap 的 Slope 工具里有 Z factor 参数QGIS 的坡度工具也有类似设置。4.3 填洼后河网位置偏移现象填洼后提取的河网和实际河流走向对不上或者河网密度异常。原因是填洼阈值设置不当把真实洼地填平了或者流量累积阈值太小导致碎河网过多。解决方法是先用小阈值填洼对比填洼前后的高程差异确认没有大面积改变地形流量累积阈值从大到小试找到河网密度合理的值。4.4 等高线提取出现大量碎线现象生成的等高线 shp 里有很多短小、不闭合的线段。原因是 DEM 本身有噪声或者等高距设得太小。解决方法是对 DEM 先做一次平滑滤波或者把等高距调大。如果只是展示用可以在提取后按长度过滤掉短于某个阈值的线段。4.5 裁剪后统计值异常现象裁剪后的 DEM 最小值出现负数或者最大值远超实际高程。原因是 NoData 值没有被正确识别或者裁剪时混入了原始 DEM 的异常像元。解决方法是在裁剪前先检查原始 DEM 的 NoData 定义裁剪后用栅格统计工具确认有效像元范围必要时手动设置 NoData 值再重新统计。5. 进阶用法用这份 DEM 做地形起伏度和剖面分析地形起伏度是描述区域地形特征的常用指标计算方式是在一个滑动窗口内取高程最大值和最小值的差。30m 数据做起伏度分析窗口大小一般选 3×3 到 9×9对应 90 米到 270 米的邻域范围。湘西州的地形起伏度整体偏高武陵山区核心地带的起伏度可以达到 300 米以上这个值在做选址和规划时很有参考意义。# 用 scipy 计算地形起伏度 import rasterio import numpy as np from scipy.ndimage import maximum_filter, minimum_filter with rasterio.open(xiangxi_dem_clipped.tif) as src: dem src.read(1) nodata src.nodata dem_valid np.where(dem nodata, np.nan, dem) # 3x3 窗口的起伏度 window_size 3 max_dem maximum_filter(dem_valid, sizewindow_size) min_dem minimum_filter(dem_valid, sizewindow_size) relief max_dem - min_dem # 保存结果 with rasterio.open(relief_3x3.tif, w, **src.meta) as dst: dst.write(relief.astype(np.float32), 1)maximum_filter和minimum_filter分别取窗口内的最大值和最小值相减得到起伏度。窗口大小根据分析尺度调整做县域级分析用 3×3 就够了做流域级可以用 9×9。注意 NoData 区域在滤波后会产生边缘效应需要在结果里重新标记 NoData。剖面分析是另一个实用技巧沿着某条线路提取高程剖面可以直观看到地形起伏。在 QGIS 里用 Profile Tool 插件ArcMap 里用 3D Analyst 的 Profile Graph。我一般会沿着规划线路或者地质剖面线提取输出高程和距离的对应关系用来判断线路的爬升和下降情况。湘西州的山地线路剖面往往呈现剧烈的锯齿状这时候要结合起伏度数据一起看判断哪些路段需要重点处理。从那以后我每次拿到新的 DEM 数据都会先跑一遍坐标系检查、范围对齐和统计值验证这三步确认没问题再进入分析流程。这套习惯帮我省掉了大量返工时间也希望帮到你。本文还有配套的精品资源点击获取 湘西州30m DEM数据处理全流程：从坐标系检查到地形分析实战
简介这份资源面向地理信息、城乡规划、测绘及环境研究等领域的从业者与学习者提供湖南省湘西土家族苗族自治州30米分辨率的DEM数字高程数据并附带本市级行政范围矢量边界可用于地形分析、坡度计算、洪水模拟与制图等GIS任务。压缩包共12个文件约50.18MB主要包含tif格式的高程栅格数据、shp格式的行政边界矢量文件以及prj坐标系统、tfw坐标配准、dbf属性表、shx与sbn/sbx索引、ovr金字塔缩略图和xml元数据等配套文件构成一套可直接在ArcGIS、QGIS中加载的完整空间数据集。目前已有301人学习下载。数据覆盖范围完整、分辨率适中既能满足区域尺度地形建模需求也便于与行政边界叠加开展空间统计适合作为课程实验、项目底图或规划分析的基础数据使用。1. 湘西州 30m DEM 数据到手后先搞清楚它能干什么、不能干什么如果你正在做湘西土家族苗族自治州范围内的地形分析、水文建模、选址评估或者三维可视化大概率绕不开一份 30m 分辨率的 DEM 数字高程数据。这次拆的这份资源是覆盖湘西州本市级范围的 30m DEM 栅格外加配套的 shp 边界文件。很多人拿到这类压缩包的第一反应是直接拖进 ArcMap 或者 QGIS 里看效果但真正决定后面顺不顺的是你在打开之前有没有想清楚三件事这份数据的精度够不够支撑你的分析尺度、shp 边界和栅格范围对不对得上、以及你打算用它做提取还是做裁剪。30m 分辨率意味着每个像元代表地面约 900 平方米的范围这个尺度适合流域级、县域级的坡度坡向分析、汇流计算和地形起伏度统计但拿它去做地块级的土方量计算或者精细的场地平整设计精度就不够了。湘西州地处武陵山区地形起伏大、沟谷密集30m 数据在描述山脊线和沟谷走向时会有一定的平滑效应这一点在做水文分析时要特别留意。配套的 shp 文件通常是行政边界或者本市级范围面图层它的作用不是装饰而是帮你把分析范围卡死在目标区域内避免跑出边界算出无意义的结果。适合谁用做区域规划、自然资源调查、遥感预处理、教学演示的从业者以及需要快速拿到一份可用地形底图的新手。不适合谁做厘米级工程测量或者城市级精细建模的团队这个精度撑不住。2. 把 DEM 和 shp 正确装进工作流从坐标系检查到范围对齐2.1 先验坐标系别等裁剪完才发现偏了几百米拿到栅格和矢量数据第一步永远不是裁剪而是确认两者的坐标系是否一致。DEM 常见的是地理坐标系 WGS84 或者 CGCS2000单位是度shp 边界可能是投影坐标系单位是米。如果直接拿投影坐标的 shp 去裁地理坐标的栅格ArcMap 会弹警告QGIS 可能直接给你裁出一片空白或者错位的结果。我一般会先在 ArcMap 里右键图层看属性或者用 QGIS 的图层属性面板确认 CRS。如果两者不一致优先把 shp 投影到和 DEM 相同的坐标系而不是反过来动栅格因为重采样栅格会引入额外的高程误差。# 用 rasterio 和 geopandas 快速检查 DEM 与 shp 的坐标系 import rasterio import geopandas as gpd dem_path xiangxi_dem_30m.tif shp_path xiangxi_boundary.shp with rasterio.open(dem_path) as src: print(DEM CRS:, src.crs) print(DEM 范围:, src.bounds) print(DEM 分辨率:, src.res) gdf gpd.read_file(shp_path) print(SHP CRS:, gdf.crs) print(SHP 范围:, gdf.total_bounds)这段代码的逻辑很直接rasterio 读栅格的 CRS、边界和分辨率geopandas 读矢量的 CRS 和总范围。参数上重点看src.res30m 数据在地理坐标系下大约是 0.000277 度左右如果显示的是 30 或者 0.000277 都算正常但单位不同后续处理方式完全不同。total_bounds返回的是 shp 的四至范围拿它和 DEM 的 bounds 对比如果 shp 范围明显大于 DEM说明边界文件可能包含了州外区域裁剪时要以 DEM 实际覆盖为准。2.2 用面图层裁剪 DEMExtract by Mask 和 Clip 的区别要分清这是搜索热词里反复出现的问题ArcMap 中依靠面图层裁剪 DEM 栅格和依靠面图层掩码提取到底有什么区别。简单说Clip工具是几何裁剪输出栅格的像元值不变只是把边界外的像元去掉边界上的像元按面积比例处理Extract by Mask是掩码提取边界外的像元被设为 NoData边界内的像元值原样保留。两者在边界规整、范围完全包含的情况下结果几乎一样但当 shp 边界和栅格像元边界不重合时Clip 可能会在边缘产生半个像元的问题而 Extract by Mask 更干净。# 用 rasterio 的 mask 功能实现按面裁剪 DEM import rasterio from rasterio.mask import mask import geopandas as gpd with rasterio.open(xiangxi_dem_30m.tif) as src: gdf gpd.read_file(xiangxi_boundary.shp) # 确保 shp 和 DEM 坐标系一致 if gdf.crs ! src.crs: gdf gdf.to_crs(src.crs) geoms [geom for geom in gdf.geometry] out_image, out_transform mask(src, geoms, cropTrue) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(xiangxi_dem_clipped.tif, w, **out_meta) as dest: dest.write(out_image)这段代码的关键参数是cropTrue它会把输出栅格的范围收缩到 shp 的实际边界而不是保留原始 DEM 的完整范围。out_meta继承原栅格的元数据后更新了行列数和仿射变换保证输出文件的地理定位正确。如果你在 ArcMap 里操作对应的是 Spatial Analyst 工具条下的 Extract by Mask输入栅格选 DEM掩码选 shp输出就是裁剪后的结果。注意一点如果 shp 有多个面要素mask 函数会取所有面的并集作为掩码范围这和 ArcMap 里逐个要素裁剪的逻辑不同需要分开处理时得先拆分 shp。2.3 裁剪后必做的三项检查范围、NoData 和统计值裁剪完不是就结束了我见过太多人裁完直接拿去算坡度结果发现边缘一圈全是异常值。第一项检查是范围把裁剪后的栅格和 shp 叠在一起看边界应该严丝合缝不能出现栅格超出 shp 或者 shp 内有大片空白的情况。第二项检查是 NoData 值Extract by Mask 会把边界外设为 NoData但有些 DEM 原始数据本身在河流或者建筑区域就有 NoData裁剪后这些值会混在一起需要在属性里确认 NoData 的数值。第三项是统计值用栅格统计工具看最小、最大和均值湘西州的高程范围大致在 200 到 1700 米之间如果最小值出现负数或者最大值超过 2000大概率是数据本身有问题或者裁剪时混入了异常像元。# 裁剪后快速统计高程分布 import rasterio import numpy as np with rasterio.open(xiangxi_dem_clipped.tif) as src: data src.read(1) nodata src.nodata if nodata is not None: valid data[data ! nodata] else: valid data.flatten() print(最小值:, np.min(valid)) print(最大值:, np.max(valid)) print(均值:, np.mean(valid)) print(有效像元数:, valid.size)这段代码先读取第一个波段然后根据 nodata 值过滤掉无效像元再算统计量。参数上注意src.nodata可能是 None这时候需要手动指定一个明显不合理的值来过滤比如 -9999。湘西州的地形统计如果均值在 600 到 800 米之间、标准差在 200 米左右说明数据分布合理如果均值异常低或者标准差极小可能是裁剪时把大片 NoData 当成了有效值。3. 从 DEM 到 shp等高线、坡度分级和流域边界的提取方法3.1 等高线提取参数设错就是一堆碎线从 DEM 提取等高线是常见需求ArcMap 里用 Contour 工具QGIS 里用等值线提取。核心参数是等高距30m 数据建议等高距设在 10 到 20 米之间设太小会生成大量碎线设太大又丢失地形细节。湘西州山地起伏大我一般用 10 米等高距做整体展示局部区域用 5 米。提取出来的等高线 shp 要检查是否闭合、是否有自相交这些在后续做地形分析时都会导致拓扑错误。# 用 GDAL 从 DEM 提取等高线 from osgeo import gdal, ogr dem gdal.Open(xiangxi_dem_clipped.tif) band dem.GetRasterBand(1) driver ogr.GetDriverByName(ESRI Shapefile) out_ds driver.CreateDataSource(contour_10m.shp) out_layer out_ds.CreateLayer(contour, geom_typeogr.wkbLineString) out_layer.CreateField(ogr.FieldDefn(elev, ogr.OFTReal)) gdal.ContourGenerate(band, 10, 0, [], 0, 0, out_layer, 0, 0) out_ds NoneContourGenerate的参数依次是波段、等高距、基准高程、忽略值列表、是否用 NoData、是否用固定间隔、输出图层、字段索引等。这里等高距设为 10基准高程设为 0表示从 0 开始每 10 米生成一条线。生成的 shp 里 elev 字段记录了每条等高线的高程值方便后续按高程筛选或者做标注。3.2 坡度分级和坡向分析别直接拿原始坡度做规划图DEM 算坡度是基础操作但直接输出的坡度栅格是连续值做规划图或者统计报告时需要分级。湘西州的山地坡度分级一般按 0-5 度、5-15 度、15-25 度、25-35 度、35 度以上来分对应平缓、较缓、中等、陡峭和极陡。ArcMap 里用 Reclassify 工具QGIS 里用栅格重分类。坡向分析同理输出的是 0 到 360 度的方向值需要按北、东北、东等八个方向重分类。注意坡度计算前要确认 DEM 的 Z 单位地理坐标系下 Z 单位是米而 XY 单位是度直接算坡度会得到错误结果需要先用投影坐标系或者设置 Z 因子。# 用 GDAL 计算坡度并重分类 from osgeo import gdal import numpy as np dem gdal.Open(xiangxi_dem_clipped.tif) gt dem.GetGeoTransform() # 地理坐标系下需要设置 Z 因子约等于 1/111320 z_factor 1.0 / 111320.0 slope_ds gdal.DEMProcessing(slope.tif, dem, slope, zFactorz_factor) slope_data slope_ds.GetRasterBand(1).ReadAsArray() # 按 5 度间隔分级 bins [0, 5, 15, 25, 35, 90] classified np.digitize(slope_data, bins) print(各级别像元数:, [np.sum(classified i) for i in range(1, len(bins))])DEMProcessing的zFactor参数是关键地理坐标系下不设置这个值坡度会被严重放大。分级用np.digitize按边界值把连续坡度映射到离散类别输出的是 1 到 5 的整数栅格。实际项目中我会把这个分级结果再转成 shp 或者做掩码统计方便出报告。3.3 流域边界和河网提取填洼是绕不过去的一步用 DEM 做水文分析填洼是第一步也是最多人翻车的地方。原始 DEM 里存在洼地不填洼直接算流向会出现断流或者内流区。ArcMap 里用 Fill 工具QGIS 里用 Wang and Liu 填洼算法。填洼之后算流向、流量累积设定阈值提取河网再转成 shp。湘西州的喀斯特地貌区域洼地特别多填洼的阈值要适当调大否则会把大量真实洼地填平改变水文响应特征。# 用 RichDEM 做填洼和流量累积 import richdem as rd dem rd.LoadGDAL(xiangxi_dem_clipped.tif) dem_filled rd.FillDepressions(dem, epsilonTrue, in_placeFalse) flow_accum rd.FlowAccumulation(dem_filled, methodD8) rd.SaveGDAL(flow_accum.tif, flow_accum)FillDepressions的epsilonTrue表示在填洼时加入微小坡度避免平坦区域流向不确定。FlowAccumulation用 D8 算法每个像元的流向指向最陡的下坡方向流量累积值表示上游汇水像元数。阈值一般设在 1000 到 5000 之间具体看流域面积和河网密度湘西州建议从 2000 开始试。4. 避坑与排查30m DEM 处理中最容易翻车的五个地方4.1 裁剪后栅格范围对不上 shp 边界现象裁剪后的 DEM 边缘和 shp 边界有明显偏移或者 shp 内出现大片 NoData。原因通常是坐标系不一致或者 shp 本身有拓扑错误比如自相交、重叠面。解决方法是先统一坐标系再用拓扑检查工具修复 shp最后重新裁剪。如果 shp 是多个面要素确认是否需要合并成一个面再裁剪。4.2 坡度计算结果明显偏大现象算出来的坡度动辄七八十度明显不符合实际地形。原因是在地理坐标系下直接计算坡度没有设置 Z 因子。解决方法是先把 DEM 投影到投影坐标系或者在坡度工具里手动设置 Z 因子为 1/111320 左右。ArcMap 的 Slope 工具里有 Z factor 参数QGIS 的坡度工具也有类似设置。4.3 填洼后河网位置偏移现象填洼后提取的河网和实际河流走向对不上或者河网密度异常。原因是填洼阈值设置不当把真实洼地填平了或者流量累积阈值太小导致碎河网过多。解决方法是先用小阈值填洼对比填洼前后的高程差异确认没有大面积改变地形流量累积阈值从大到小试找到河网密度合理的值。4.4 等高线提取出现大量碎线现象生成的等高线 shp 里有很多短小、不闭合的线段。原因是 DEM 本身有噪声或者等高距设得太小。解决方法是对 DEM 先做一次平滑滤波或者把等高距调大。如果只是展示用可以在提取后按长度过滤掉短于某个阈值的线段。4.5 裁剪后统计值异常现象裁剪后的 DEM 最小值出现负数或者最大值远超实际高程。原因是 NoData 值没有被正确识别或者裁剪时混入了原始 DEM 的异常像元。解决方法是在裁剪前先检查原始 DEM 的 NoData 定义裁剪后用栅格统计工具确认有效像元范围必要时手动设置 NoData 值再重新统计。5. 进阶用法用这份 DEM 做地形起伏度和剖面分析地形起伏度是描述区域地形特征的常用指标计算方式是在一个滑动窗口内取高程最大值和最小值的差。30m 数据做起伏度分析窗口大小一般选 3×3 到 9×9对应 90 米到 270 米的邻域范围。湘西州的地形起伏度整体偏高武陵山区核心地带的起伏度可以达到 300 米以上这个值在做选址和规划时很有参考意义。# 用 scipy 计算地形起伏度 import rasterio import numpy as np from scipy.ndimage import maximum_filter, minimum_filter with rasterio.open(xiangxi_dem_clipped.tif) as src: dem src.read(1) nodata src.nodata dem_valid np.where(dem nodata, np.nan, dem) # 3x3 窗口的起伏度 window_size 3 max_dem maximum_filter(dem_valid, sizewindow_size) min_dem minimum_filter(dem_valid, sizewindow_size) relief max_dem - min_dem # 保存结果 with rasterio.open(relief_3x3.tif, w, **src.meta) as dst: dst.write(relief.astype(np.float32), 1)maximum_filter和minimum_filter分别取窗口内的最大值和最小值相减得到起伏度。窗口大小根据分析尺度调整做县域级分析用 3×3 就够了做流域级可以用 9×9。注意 NoData 区域在滤波后会产生边缘效应需要在结果里重新标记 NoData。剖面分析是另一个实用技巧沿着某条线路提取高程剖面可以直观看到地形起伏。在 QGIS 里用 Profile Tool 插件ArcMap 里用 3D Analyst 的 Profile Graph。我一般会沿着规划线路或者地质剖面线提取输出高程和距离的对应关系用来判断线路的爬升和下降情况。湘西州的山地线路剖面往往呈现剧烈的锯齿状这时候要结合起伏度数据一起看判断哪些路段需要重点处理。从那以后我每次拿到新的 DEM 数据都会先跑一遍坐标系检查、范围对齐和统计值验证这三步确认没问题再进入分析流程。这套习惯帮我省掉了大量返工时间也希望帮到你。本文还有配套的精品资源点击获取