如何自动计算合适的等高距并输出等高线

GISer35

在地形制图中，合理的等高距通常取决于：

地形起伏程度（高差）
比例尺（可由分辨率推断）
应用目的（工程、规划、展示）

计算等高距时，常用的经验法则：等高距 ≈ 地形标准差 / 5 ~ 10

以下代码采用QGIS的PyQGIS开发语法，参照经验法则，实现DEM等高距的计算，并输出Shp格式的等高线。

需要替换的变量：
output_dir: 输出目录
layer_name: QGIS中的图层名称

# -*- coding: utf-8 -*-
from qgis.core import (
    QgsRasterLayer,
    QgsVectorFileWriter,
    QgsProcessingContext,
    QgsProject,
    QgsProcessingFeedback,
    QgsApplication
)
import os
import processing  # 确保 Processing 工具箱可用

# 初始化 Processing 框架（如果尚未加载）
if 'processing' not in globals():
    import processing

# 设置输出路径
output_dir = r"D:/output"
os.makedirs(output_dir, exist_ok=True)
output_path = os.path.join(output_dir, "contours.shp")

# 查找名为 "dsm" 的图层
layer_name = "dsm"
dsm_layer = None
for layer in QgsProject.instance().mapLayers().values():
    if layer.name() == layer_name and isinstance(layer, QgsRasterLayer):
        dsm_layer = layer
        break

if not dsm_layer:
    raise Exception(f"未找到名为 '{layer_name}' 的栅格图层！")

print(f"使用图层: {dsm_layer.name()}")

# 获取波段统计信息
stats = dsm_layer.dataProvider().bandStatistics(1)
min_elev = stats.minimumValue
max_elev = stats.maximumValue
std_dev = stats.stdDev

print(f"最小高程: {min_elev:.2f}, 最大高程: {max_elev:.2f}, 标准差: {std_dev:.2f}")

# 自动计算等高距：取标准差的 1/10，并四舍五入到最接近的 5 或 10 的倍数
raw_interval = std_dev / 10
interval = max(round(raw_interval / 5) * 5, 5)  # 四舍五入到最近的 5 的倍数，最小为 5
print(f"建议等高距: {interval} 米")

# 设置起始等高线为大于最小值的第一个 interval 的整数倍
start_level = ((min_elev // interval) + 1) * interval
print(f"起始等高线: {start_level}")

# 准备参数调用 GDAL contour 算法
params = {
    'INPUT': dsm_layer,
    'BAND': 1,
    'INTERVAL': interval,
    'FIELD_NAME': 'ELEV',
    'CREATE_3D': False,
    'USE_NODATA': False,
    'NODATA': None,
    'OFFSET': 0,
    'OUTPUT': output_path
}

context = QgsProcessingContext()
feedback = QgsProcessingFeedback()

# 执行 GDAL contour 算法
result = processing.run("gdal:contour", params, context=context, feedback=feedback)

# 输出结果
if result and os.path.exists(output_path):
    print(f"✅ 等高线已成功生成: {output_path}")
    
    # 加载结果到地图
    contour_layer = iface.addVectorLayer(output_path, "等高线", "ogr")
    if contour_layer.isValid():
        print("🎉 等高线图层已添加到地图中。")
else:
    print("❌ 等高线生成失败，请检查日志。")