基于道路缓冲区的滑坡风险分析

GISer35

接着上一篇文章中准备的数据，本文继续讲解如何计算基于道路缓冲区的滑坡风险分析。

任务目标

评估道路周边区域的滑坡易发性：

靠近道路的区域（人为扰动大）
坡度较陡的区域
特定坡向（如南向日照强、干燥或北向潮湿等，视地区而定）

分析思路

我们将采用 加权叠加法（Weighted Overlay） 构建滑坡风险指数图：
1. 距道路距离
数据来源：道路图层 → 缓冲+欧式距离
原理：越靠近道路，风险越高

2. 坡度
数据来源：坡度图层
原理：坡度越大，风险越高

3. 坡向
数据来源：坡度图层
原理：可重分类为日晒/背阴面或直接用于模型

执行PyQGIS代码
以下代码在QGIS的Python控制台中运行，生成道路风险、坡度风险、坡向风险图层。请根据不同情况，微调参数。
注意：请根据运行环境，修改输出目录。

# -*- coding: utf-8 -*-
from qgis.core import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
import processing
import os

# 初始化处理框架
if 'processing' not in globals():
    from processing.core.Processing import Processing
    Processing.initialize()

# 设置输出目录
output_dir = r"d:/output"
os.makedirs(output_dir, exist_ok=True)

# 获取项目实例
project = QgsProject.instance()
root = project.layerTreeRoot()

# === 第一步：获取所有需要的图层 ===
try:
    slope_layer = project.mapLayersByName("坡度")[0]
    distance_to_road = project.mapLayersByName("距道路距离（米）")[0]
    aspect_layer = project.mapLayersByName("坡向")[0] if project.mapLayersByName("坡向") else None
except IndexError as e:
    raise Exception("缺少必要图层，请确保存在【坡度】和【距道路距离（米）】图层。")

# 共享空间参数（来自坡度图层）
cell_size = slope_layer.rasterUnitsPerPixelX()
crs = slope_layer.crs()
extent = slope_layer.extent()  # 使用统一地理范围

print(f"使用分辨率: {cell_size:.6f} 米, CRS: {crs.authid()}")

# === Step 1: 对【坡度】重分类（越高风险越大）===
# 规则（通用山区标准）：
# [0–15°]: 低风险 (1)
# [15–25°]: 中 (2)
# [25–35°]: 高 (3)
# >35°: 极高 (4)

print("正在重分类【坡度】...")
reclassified_slope = os.path.join(output_dir, "risk_slope_reclass.tif")

params_slope = {
    'INPUT_RASTER': slope_layer,
    'RASTER_BAND': 1,
    'TABLE': [0, 15, 1, 15, 25, 2, 25, 35, 3, 35, 90, 4],
    'RANGE_BOUNDARIES': 1,  # min <= value < max
    'NODATA_FOR_MISSING': False,
    'DATA_TYPE': 5,  # Float32
    'OUTPUT': reclassified_slope
}
result_slope = processing.run("native:reclassifybytable", params_slope)
slope_risk = QgsRasterLayer(result_slope['OUTPUT'], "坡度风险")
QgsProject.instance().addMapLayer(slope_risk)

# === Step 2: 对【距道路距离（米）】重分类（越近风险越高）===
# 设定：
# >500m: 低风险 (1)
# 200–500m: 中 (2)
# 100–200m: 高 (3)
# <100m: 极高 (4)

print("正在重分类【距道路距离】...")
reclassified_road = os.path.join(output_dir, "risk_road_reclass.tif")

params_road = {
    'INPUT_RASTER': distance_to_road,
    'RASTER_BAND': 1,
    'TABLE': [
        0, 100, 4,
        100, 200, 3,
        200, 500, 2,
        500, 10000, 1
    ],
    'RANGE_BOUNDARIES': 1,
    'NODATA_FOR_MISSING': False,
    'DATA_TYPE': 5,
    'OUTPUT': reclassified_road
}
result_road = processing.run("native:reclassifybytable", params_road)
road_risk = QgsRasterLayer(result_road['OUTPUT'], "道路风险")
QgsProject.instance().addMapLayer(road_risk)

# === Step 3: （可选）对【坡向】重分类 —— 北坡潮湿更危险（以中国南方为例）===
# 南方地区北坡阴湿 → 易滑坡
# 重分类规则：
# 0–22.5 或 337.5–360 → 北坡 → 风险高 (3)
# 157.5–202.5 → 南坡 → 干燥稳定 → 风险低 (1)
# 其他 → 中等 (2)

aspect_risk = None
if aspect_layer:
    print("正在重分类【坡向】...")
    reclassified_aspect = os.path.join(output_dir, "risk_aspect_reclass.tif")

    params_aspect = {
        'INPUT_RASTER': aspect_layer,
        'RASTER_BAND': 1,
        'TABLE': [
            0, 22.5, 3,
            22.5, 157.5, 2,
            157.5, 202.5, 1,
            202.5, 337.5, 2,
            337.5, 360, 3
        ],
        'RANGE_BOUNDARIES': 2,  # ≤ ≤
        'NODATA_FOR_MISSING': True,
        'NO_DATA': 0,
        'DATA_TYPE': 5,
        'OUTPUT': reclassified_aspect
    }
    result_aspect = processing.run("native:reclassifybytable", params_aspect)
    aspect_risk = QgsRasterLayer(result_aspect['OUTPUT'], "坡向风险")
    QgsProject.instance().addMapLayer(aspect_risk)

# === Step 4: 加权叠加计算综合风险指数 ===
# 权重建议（可根据区域调整）：
# - 坡度：50%
# - 道路距离：40%
# - 坡向：10%

print("正在进行加权叠加分析...")

# 准备参与计算的图层列表
input_layers = [
    (road_risk, 'A'),   # A@1 = 道路风险
    (slope_risk, 'B'),  # B@1 = 坡度风险
]
if aspect_risk:
    input_layers.append((aspect_risk, 'C'))

# 创建 Raster Calculator Entries
entries = []
for layer, prefix in input_layers:
    entry = QgsRasterCalculatorEntry()
    entry.ref = f"{prefix}@1"
    entry.raster = layer
    entry.bandNumber = 1
    entries.append(entry)

# 表达式：0.4*A + 0.5*B [+ 0.1*C]
expression_parts = ["0.4 * A@1", "0.5 * B@1"]
if aspect_risk:
    expression_parts.append("0.1 * C@1")
expression = " + ".join(expression_parts)

output_index = os.path.join(output_dir, "landslide_risk_index.tif")

calc = QgsRasterCalculator(
    expression,
    output_index,
    'GTiff',
    extent,
    cell_size,
    cell_size,
    entries
)
res = calc.processCalculation()
if res != 0:
    raise RuntimeError(f"栅格计算器运行失败，错误码: {res}")

risk_index = QgsRasterLayer(output_index, "滑坡风险指数")
QgsProject.instance().addMapLayer(risk_index)

加权叠加法

权重设定：坡度 50%，道路 40%，坡向 10%。需根据不同情况，调整权重比例。

成功生成坡度风险、道路风险、坡向风险三个图层后，运行以下代码。得到最终滑坡风险指数图。

# -*- coding: utf-8 -*-
from qgis.core import (
    QgsProject, QgsRectangle, QgsCoordinateReferenceSystem, QgsCoordinateTransformContext
)
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
import os

# 设置输出路径
output_dir = r"d:/output"
os.makedirs(output_dir, exist_ok=True)

# 获取项目与图层
project = QgsProject.instance()
layers = project.mapLayersByName

# === 第一步：获取输入图层（使用您已经存在的重分类结果）===
try:
    slope_risk_layer = layers("坡度风险")[0]   # 【坡度风险】
    road_risk_layer = layers("道路风险")[0]     # 【道路风险】
    aspect_risk_layer = layers("坡向风险")[0] if layers("坡向风险") else None  # 可选
except IndexError as e:
    raise Exception(f"缺少必要图层，请检查名称是否准确：{e}")

# 使用统一的空间参数（来自任意图层，比如坡度风险）
extent = slope_risk_layer.extent()        # QgsRectangle
crs = slope_risk_layer.crs()              # EPSG:32648
cell_size = slope_risk_layer.rasterUnitsPerPixelX()

# 计算输出行列数（必须为整数！）
n_cols = int(extent.width() / cell_size)
n_rows = int(extent.height() / cell_size)

print(f"分辨率: {cell_size:.6f}")
print(f"范围宽高: {extent.width():.2f} x {extent.height():.2f}")
print(f"行列数: {n_cols} x {n_rows}")

# === 第二步：准备栅格计算器条目 ===
entries = []

# A@1 → 坡度风险
entry_a = QgsRasterCalculatorEntry()
entry_a.ref = "A@1"
entry_a.raster = slope_risk_layer
entry_a.bandNumber = 1
entries.append(entry_a)

# B@1 → 道路风险
entry_b = QgsRasterCalculatorEntry()
entry_b.ref = "B@1"
entry_b.raster = road_risk_layer
entry_b.bandNumber = 1
entries.append(entry_b)

# C@1 → 坡向风险（如果存在）
if aspect_risk_layer:
    entry_c = QgsRasterCalculatorEntry()
    entry_c.ref = "C@1"
    entry_c.raster = aspect_risk_layer
    entry_c.bandNumber = 1
    entries.append(entry_c)

# === 第三步：定义加权公式（示例权重）===
# 权重建议：坡度 50%，道路 40%，坡向 10%
if aspect_risk_layer:
    formula = "0.5 * A@1 + 0.4 * B@1 + 0.1 * C@1"
else:
    formula = "0.5 * A@1 + 0.5 * B@1"

# 输出文件
output_file = os.path.join(output_dir, "landslide_risk_index.tif")

# === 第四步：创建并运行栅格计算器（关键：确保参数类型正确！）===

# 创建坐标变换上下文（推荐，避免警告）
transform_context = QgsCoordinateTransformContext()

calc = QgsRasterCalculator(
    formula,
    output_file,
    'GTiff',           # 输出格式
    extent,            # 输出范围 (QgsRectangle)
    n_cols,            # 列数 (int) ✅
    n_rows,            # 行数 (int) ✅
    entries,           # 图层条目列表
    transform_context   # 推荐添加
)

# 执行计算
result = calc.processCalculation()
if result == QgsRasterCalculator.Success:
    print("✅ 栅格计算成功完成！")
    # 加载结果到地图
    risk_index = QgsRasterLayer(output_file, "滑坡风险指数")
    project.addMapLayer(risk_index)
else:
    error_msg = calc.lastError()
    print(f"❌ 栅格计算失败，错误码: {result}, 错误信息: {error_msg}")

个人知识有限，内容必有错误，请谨慎参考。欢迎留言交流或联系我们。转载须注明出处。