分享计算道路转弯半径的方法

GISer35

在交通工程、自动驾驶、BIM/CIM 和道路安全评估中，最小转弯半径（Minimum Turning Radius） 是衡量道路几何线形可通行性的核心指标——它直接决定车辆（尤其是大型车）能否安全、不擦碰地通过弯道。

QGIS 原生不提供自动识别道路弯道并计算其曲率半径的一键算法（如 ArcGIS 的 Calculate Curve Radius 或 Civil 3D 的 Alignment Analysis）。但我们可以用 QGIS 现有工具 + 几何原理，同样实现该目标。

什么是“道路转弯半径”？

对矢量线要素，转弯半径并非整条线的属性，而是沿线上局部弯曲段的几何特征。我们采用通用定义：

曲率（Curvature）：描述线在某点处弯曲程度：κ = 1/R，单位 1/m；R 越小 → 弯越急；
最小转弯半径 R_min：矢量线所有局部弯曲段中，对应最大曲率 κmax 的倒数：R_min = 1 / κmax（单位：米）

步骤1：提取道路所有折点

打开 Processing Toolbox → 搜索 extract vertices → 运行。

运行后得到图层包括含字段：

vertex_index: 折点在线上的序号（0=起点，-1=终点）
distance: 沿线累计距离（米）
angle: 折点处切线方位角（度）

步骤2：为每个折点计算曲率半径

原理：对任意折点 P_i，取其前后相邻两点 P{i-1} 和 P{i+1}，构成三点。若三点不共线，则唯一确定一个圆，其半径即为该点处的近似转弯半径。
公式（三点外接圆半径）：
R = (a × b × c) / (4 × Δ)

其中 a,b,c 为三点两两距离，Δ 为三角形面积。

在 QGIS 的 Python 控制台中运行以下代码，计算曲率。
注意：将顶点图层修改为你的点图层名称。


# -*- coding: utf-8 -*-
"""
【QGIS PyQGIS 脚本】为 '顶点' 图层计算每个折点的局部转弯半径（R = d / (2·sin(θ/2))）
✅ 适配 QGIS 3.40.9 + EPSG:32648（米制坐标）
✅ 输入：图层名 '顶点'（必须含字段：vertex_index, distance, angle）
✅ 输出：新增字段 'radius_m'（double），并自动选中最小半径点
⚠️ 注意：请确保 '顶点' 图层已加载
"""

from qgis.core import QgsProject, QgsVectorLayer, QgsField, QgsExpression, QgsFeatureRequest
from qgis.PyQt.QtCore import QVariant
import math

# === 步骤 1：获取 '顶点' 图层 ===
layer_name = "顶点"
layer = QgsProject.instance().mapLayersByName(layer_name)
if not layer:
    raise ValueError(f"❌ 错误：未找到图层 '{layer_name}'。请确认图层名称是否准确（区分大小写/中文）。")
layer = layer[0]

if not isinstance(layer, QgsVectorLayer) or not layer.isValid():
    raise ValueError(f"❌ 错误：'{layer_name}' 不是有效的矢量图层。")

print(f"✅ 正在处理图层：{layer_name}（共 {layer.featureCount()} 个点）")

# === 步骤 2：检查必要字段 ===
required_fields = ["vertex_index", "angle"]
missing = [f for f in required_fields if f not in layer.fields().names()]
if missing:
    raise ValueError(f"❌ 错误：图层缺少必要字段 → {missing}")

# === 步骤 3：添加 radius_m 字段（如不存在）===
field_name = "radius_m"
if layer.fields().lookupField(field_name) == -1:
    print(f"➕ 正在添加字段：{field_name}（Decimal）")
    layer.startEditing()
    layer.addAttribute(QgsField(field_name, QVariant.Double, len=10, prec=2))
    layer.commitChanges()
    print(f"✅ 字段 '{field_name}' 添加成功")
else:
    print(f"ℹ️  字段 '{field_name}' 已存在，将覆盖计算值")

# === 步骤 4：批量计算 radius_m（核心算法）===
layer.startEditing()
feat_count = layer.featureCount()
print(f"⚙️  开始计算曲率半径（共 {feat_count} 个点）...")

# 预加载所有要素（提升性能，避免重复遍历）
all_feats = list(layer.getFeatures())
id_to_feat = {f.id(): f for f in all_feats}

# 缓存 geometry 和属性（关键：避免多次调用 $geometry）
geom_cache = {}
attr_cache = {}
for f in all_feats:
    geom_cache[f.id()] = f.geometry()
    attr_cache[f.id()] = {field.name(): f[field.name()] for field in layer.fields()}

# 计算函数（三点弦长-偏转角法）
def calc_radius_m(prev_id, curr_id, next_id):
    # 检查 ID 是否合法
    if prev_id not in geom_cache or next_id not in geom_cache:
        return None
    
    g_prev = geom_cache[prev_id]
    g_curr = geom_cache[curr_id]
    g_next = geom_cache[next_id]
    
    if not (g_prev and g_curr and g_next):
        return None
    
    # 获取坐标（EPSG:32648 → 单位为米）
    try:
        p1 = g_prev.asPoint()
        p2 = g_curr.asPoint()
        p3 = g_next.asPoint()
    except:
        return None
    
    # 计算向量 P1→P2 和 P2→P3 的方位角（度）
    try:
        a1 = math.degrees(math.atan2(p2.y() - p1.y(), p2.x() - p1.x()))
        a2 = math.degrees(math.atan2(p3.y() - p2.y(), p3.x() - p2.x()))
    except:
        return None
    
    # 归一化夹角到 [-180, 180]
    delta_a = (a2 - a1 + 180) % 360 - 180
    abs_delta_a = abs(delta_a)
    
    # 若转向角过小（< 0.3°），视为直线 → 设极大半径（避免除零 & 合理物理意义）
    if abs_delta_a < 0.3:
        return 100000.0
    
    # 计算弦长 d = |P1P3|（米）
    d = math.sqrt((p3.x() - p1.x())**2 + (p3.y() - p1.y())**2)
    if d < 0.1:  # 异常重合点
        return None
    
    # R = d / (2 * sin(|Δθ|/2))   ← 弧度制
    try:
        r = d / (2 * math.sin(math.radians(abs_delta_a / 2)))
        return max(r, 5.0)  # 物理下限：车辆最小理论转弯半径约 5m（小轿车）
    except (ZeroDivisionError, ValueError):
        return None

# 批量更新
min_radius = float('inf')
min_fid = None

for i, feat in enumerate(all_feats):
    vid = feat["vertex_index"]
    fid = feat.id()
    
    # 跳过首点（vertex_index == 0）和末点（vertex_index == -1）
    if vid == 0 or vid == -1:
        layer.changeAttributeValue(fid, layer.fields().lookupField(field_name), None)
        continue
    
    # 查找前一点与后一点（按 feature id 顺序 ≈ vertex_index 顺序，但更可靠的是用 vertex_index）
    # ✅ 更鲁棒方式：按 vertex_index 排序后取邻点（因 vertex_index 是严格递增整数序列）
    pass  # 我们改用「按 vertex_index 排序」策略（见下方优化）

# === 优化：按 vertex_index 排序，确保逻辑顺序（解决 id 乱序问题）===
sorted_feats = sorted(
    all_feats,
    key=lambda f: (f["LAYER"] if "LAYER" in f.fields().names() else "", f["ID"] if "ID" in f.fields().names() else 0, f["vertex_index"])
)

# 构建 vertex_index → feature 映射（同一道路内）
from collections import defaultdict
road_groups = defaultdict(list)
for f in sorted_feats:
    key = (f["LAYER"], f["ID"]) if "LAYER" in f.fields().names() and "ID" in f.fields().names() else ("unknown", 0)
    road_groups[key].append(f)

# 逐条道路内计算（保证 vertex_index 连续性）
print("🔄 正在按道路分组计算...")
for road_key, feats_in_road in road_groups.items():
    if len(feats_in_road) < 3:
        continue
    # 按 vertex_index 排序（升序）
    feats_sorted = sorted(feats_in_road, key=lambda f: f["vertex_index"])
    
    for idx, curr_feat in enumerate(feats_sorted):
        vid = curr_feat["vertex_index"]
        if vid == 0 or vid == -1:
            continue
            
        # 取前一点（idx-1）、后一点（idx+1）
        if idx == 0 or idx == len(feats_sorted) - 1:
            continue
            
        prev_feat = feats_sorted[idx-1]
        next_feat = feats_sorted[idx+1]
        
        r = calc_radius_m(prev_feat.id(), curr_feat.id(), next_feat.id())
        if r is not None:
            layer.changeAttributeValue(curr_feat.id(), layer.fields().lookupField(field_name), round(r, 2))
            if r < min_radius:
                min_radius = r
                min_fid = curr_feat.id()

# === 步骤 5：提交编辑并高亮最小半径点 ===
layer.commitChanges()
print(f"✅ 所有半径计算完成！最小转弯半径 = {min_radius:.2f} m")

if min_fid is not None:
    layer.selectByIds([min_fid])
    print(f"📍 最小半径点已选中（feature ID = {min_fid}）")
else:
    print("⚠️  未找到有效最小半径点（可能全为 NULL）")

# === 可选：弹窗提示结果 ===
from qgis.PyQt.QtWidgets import QMessageBox
QMessageBox.information(
    None,
    "✅ 计算完成",
    f"已为 '{layer_name}' 图层计算曲率半径。\n\n🔹 最小转弯半径：{min_radius:.2f} 米\n🔹 位置：已选中对应点\n🔹 字段：'{field_name}'（可在属性表查看）"
)

执行后结果，如下图所示。

更多应用问题，欢迎留言或联系我们。转载须注明出处。