土石流区域設定支援システムと同等な横断測線、流下方向（縦断測線）をQGISで再現

はじめに

土石流区域設定支援システム（dem.mdb）のCSVデータより、横断測線、縦断測線をQGIS上で可視化・活用するためのPythonスクリプトを作成しました。

システム独自の「カンマ区切り座標」や「mm単位」を自動変換し、以下の2つのレイヤを自動生成します。

  1. 流下方向（縦断）：勾配に応じて色分けし、始点・終点の座標を持たせたライン
  2. 横断測線：各測線の始点と終点を結び、連番を振ったライン 

1. 流下方向（縦断）の作成ツール

縦断データのCSVを読み込み、以下の処理を行います。

  • mm単位の座標をm単位に変換してラインを作成
  • 勾配（角度）を計算し、自動で色分け（0～2度、2～3度...30度以上）
  • 属性テーブルに「勾配角度」と「始点・終点座標(X,Y)」を出力

スクリプトコード

import csv
import math
from qgis.core import (
    QgsProject,
    QgsVectorLayer,
    QgsFeature,
    QgsGeometry,
    QgsField,
    QgsPointXY,
    QgsGraduatedSymbolRenderer,
    QgsRendererRange,
    QgsSymbol,
    QgsCoordinateReferenceSystem
)
from qgis.utils import iface
from PyQt5.QtWidgets import QFileDialog, QMessageBox
from PyQt5.QtCore import QVariant
from PyQt5.QtGui import QColor

def main_process():
    # 1. ファイル選択
    file_path, _ = QFileDialog.getOpenFileName(
        None, "CSVファイルを選択してください", "", "CSV Files (*.csv);;All Files (*)"
    )
    if not file_path: return

    # 2. データを読み込み
    points = []
    encodings = ['cp932', 'utf-8']
    data_rows = []
    
    for enc in encodings:
        try:
            with open(file_path, 'r', encoding=enc, newline='') as f:
                reader = csv.reader(f)
                data_rows = list(reader)
            break
        except UnicodeDecodeError: continue
            
    if not data_rows: return

    # 座標データ列探索
    target_col_index = -1
    start_row = 1 if len(data_rows) > 1 else 0
    for i in range(start_row, min(len(data_rows), 5)):
        row = data_rows[i]
        for idx, val in enumerate(row):
            if val.count(',') >= 3: 
                target_col_index = idx
                break
        if target_col_index != -1: break
    
    if target_col_index == -1: return

    # データ解析
    for i in range(start_row, len(data_rows)):
        row = data_rows[i]
        if len(row) <= target_col_index: continue
        val_str = row[target_col_index]
        try:
            parts = val_str.split(',')
            if len(parts) < 4: continue
            # [1]=Z, [2]=X, [3]=Y (mm->m)
            z = float(parts[1]) / 1000.0
            x = float(parts[2]) / 1000.0
            y = float(parts[3]) / 1000.0
            points.append({'x': x, 'y': y, 'z': z})
        except: continue

    if len(points) < 2: return

    # 3. レイヤ作成（流下方向）
    crs = QgsProject.instance().crs().authid()
    if not crs: crs = "EPSG:6668"
    line_layer = QgsVectorLayer(f"LineString?crs={crs}", "流下方向", "memory")
    pr = line_layer.dataProvider()
    
    # 属性追加
    pr.addAttributes([
        QgsField("slope_angle", QVariant.Double),
        QgsField("Start_X", QVariant.Double),
        QgsField("Start_Y", QVariant.Double),
        QgsField("End_X", QVariant.Double),
        QgsField("End_Y", QVariant.Double)
    ])
    line_layer.updateFields()

    new_features = []
    
    # 4. ライン生成と勾配計算
    for i in range(len(points) - 1):
        p1 = points[i]
        p2 = points[i+1]
        dx = p2['x'] - p1['x']
        dy = p2['y'] - p1['y']
        dz = p2['z'] - p1['z']
        dist = math.sqrt(dx*dx + dy*dy)
        angle = 0.0
        if dist != 0:
            angle = math.degrees(math.atan2(dz, dist))
        elif dz != 0:
            angle = 90.0 if dz > 0 else -90.0
            
        geom = QgsGeometry.fromPolylineXY([QgsPointXY(p1['x'], p1['y']), QgsPointXY(p2['x'], p2['y'])])
        f = QgsFeature()
        f.setGeometry(geom)
        f.setAttributes([angle, p1['x'], p1['y'], p2['x'], p2['y']])
        new_features.append(f)

    pr.addFeatures(new_features)
    line_layer.updateExtents()
    
    # 5. 色分け適用
    target_field = 'abs("slope_angle")'
    categories = [
        (0.0,  2.0,  QColor(0, 255, 255), "0度 ～ 2度未満"),
        (2.0,  3.0,  QColor(0, 0, 255),   "2度 ～ 3度未満"),
        (3.0,  7.0,  QColor(0, 255, 0),   "3度 ～ 7度未満"),
        (7.0,  10.0, QColor(139, 0, 0),   "7度 ～ 10度未満"),
        (10.0, 15.0, QColor(255, 0, 255), "10度 ～ 15度未満"),
        (15.0, 20.0, QColor(205, 133, 63), "15度 ～ 20度未満"),
        (20.0, 30.0, QColor(255, 165, 0),  "20度 ～ 30度未満"),
        (30.0, 90.0, QColor(255, 0, 0),    "30度以上")
    ]
    ranges = []
    for lower, upper, color, label in categories:
        symbol = QgsSymbol.defaultSymbol(line_layer.geometryType())
        symbol.setColor(color)
        symbol.setWidth(0.8)
        ranges.append(QgsRendererRange(lower, upper, symbol, label))

    renderer = QgsGraduatedSymbolRenderer(target_field, ranges)
    renderer.setMode(QgsGraduatedSymbolRenderer.Custom)
    line_layer.setRenderer(renderer)

    QgsProject.instance().addMapLayer(line_layer)
    QMessageBox.information(None, "完了", "CSVから色分けラインを作成しました。\n属性テーブルに座標(m)が入っています。")

main_process()

勾配ごとに色分けされた流下方向ライン

2. 横断測線の作成ツール

横断測線のCSVデータから、各測線（IDグループ：10000, 20000...）ごとの始点と終点を抽出し、ラインを作成します。

  • IDごとにグループ化し、最初と最後を自動抽出
  • ラインを作成し、属性に連番（Seq_No）を付与
  • 始点・終点のX,Y座標（m単位）を属性テーブルに格納

スクリプトコード

import csv
from qgis.core import (
    QgsProject,
    QgsVectorLayer,
    QgsFeature,
    QgsGeometry,
    QgsField,
    QgsPointXY
)
from PyQt5.QtWidgets import QFileDialog, QMessageBox
from PyQt5.QtCore import QVariant

def create_transverse_lines_final():
    # 1. ファイル選択
    file_path, _ = QFileDialog.getOpenFileName(
        None, "CSVファイルを選択してください", "", "CSV Files (*.csv);;All Files (*)"
    )
    if not file_path: return

    # 2. データ読み込み
    groups = {}
    encodings = ['cp932', 'utf-8']
    rows = []
    
    loaded = False
    for enc in encodings:
        try:
            with open(file_path, 'r', encoding=enc, newline='') as f:
                reader = csv.DictReader(f)
                if not reader.fieldnames: continue
                field_map = {name.strip(): name for name in reader.fieldnames}
                id_col = next((k for k in field_map if 'ID' in k), None)
                il_col = next((k for k in field_map if 'ILData' in k), None)
                if not id_col or not il_col: continue
                for row in reader:
                    rows.append({'id_col': id_col, 'il_col': il_col, 'data': row})
                loaded = True
            break
        except UnicodeDecodeError: continue
            
    if not loaded: return

    # 3. グループごとの始点・終点抽出
    for item in rows:
        row = item['data']
        id_key = item['id_col']
        try:
            current_id_val = int(row[id_key])
            group_key = current_id_val // 10000
            if group_key == 0: continue 

            if group_key not in groups:
                groups[group_key] = {
                    'min_id': current_id_val, 'min_row': row,
                    'max_id': current_id_val, 'max_row': row
                }
            else:
                if current_id_val < groups[group_key]['min_id']:
                    groups[group_key]['min_id'] = current_id_val
                    groups[group_key]['min_row'] = row
                if current_id_val > groups[group_key]['max_id']:
                    groups[group_key]['max_id'] = current_id_val
                    groups[group_key]['max_row'] = row
        except ValueError: continue

    # 4. レイヤ作成（横断測線）
    crs = QgsProject.instance().crs().authid()
    if not crs: crs = "EPSG:6668"
    vl = QgsVectorLayer(f"LineString?crs={crs}", "横断測線", "memory")
    pr = vl.dataProvider()
    
    # 属性定義
    pr.addAttributes([
        QgsField("Seq_No", QVariant.Int),
        QgsField("Group_ID", QVariant.Int),
        QgsField("Start_ID", QVariant.Int),
        QgsField("End_ID", QVariant.Int),
        QgsField("Start_X", QVariant.Double),
        QgsField("Start_Y", QVariant.Double),
        QgsField("End_X", QVariant.Double),
        QgsField("End_Y", QVariant.Double)
    ])
    vl.updateFields()
    
    new_features = []
    il_col_name = rows[0]['il_col']
    sorted_group_keys = sorted(groups.keys())

    # 5. ライン作成
    for seq_index, g_key in enumerate(sorted_group_keys):
        g_val = groups[g_key]
        
        def get_xy(row_data):
            try:
                val_str = row_data[il_col_name]
                parts = val_str.split(',')
                if len(parts) < 4: return None
                # [2]=X, [3]=Y, mm->m
                x = float(parts[2]) / 1000.0
                y = float(parts[3]) / 1000.0
                return QgsPointXY(x, y)
            except: return None

        pt_start = get_xy(g_val['min_row'])
        pt_end = get_xy(g_val['max_row'])
        
        if pt_start and pt_end:
            line_geom = QgsGeometry.fromPolylineXY([pt_start, pt_end])
            feat = QgsFeature()
            feat.setGeometry(line_geom)
            feat.setAttributes([
                seq_index, g_key * 10000, g_val['min_id'], g_val['max_id'],
                pt_start.x(), pt_start.y(), pt_end.x(), pt_end.y()
            ])
            new_features.append(feat)

    pr.addFeatures(new_features)
    vl.updateExtents()
    QgsProject.instance().addMapLayer(vl)
    QMessageBox.information(None, "完了", "レイヤ「横断測線」を作成しました。\n属性テーブルを確認してください。")

create_transverse_lines_final()

自動作成された横断測線ライン

おわりに

これらのスクリプトにより、複雑な手順を踏まずに「CSVを選択するだけ」で、土石流支援システムと同等な形状をQGIS上に再現できるようになります。

コメント

コメントを投稿

このブログの人気の投稿

日本でよく使う EPSG コード 一覧

1. 経緯度（地理）座標系 区分 名称 EPSG 主な用途・備考 世界標準（地理） WGS 84（緯度・経度） 4326 GNSS/GPS の測位やデータ交換で使用。単位は「度」。 Web地図（投影） WGS 84 / Pseudo-Mercator（Web Mercator） 3857 Google/OSM/国土地理院タイルなどの 表示 に使用。単位は「メートル」。別名：Popular Visualization / Google Mercator。 日本最新 JGD2011（緯度・経度） 6668 測地成果2011。公共測量の基準。平面直角座標系は下表参照。 過去の全国基準 JGD2000（緯度・経度） 4612 2002–2011 年に使用。JGD2011 と混在時は変換推奨。 旧日本測地系 Tokyo Datum 4301 1940年代〜2002年まで使用。WGS84とは数百mのズレあり。 2. 平面直角座標系（Gauss–Krüger 投影） 系 (ローマ数字) 系 (数字) 中央経線 JGD2011 EPSG JGD2000 EPSG 適用都道府県（概要） I 1 129°30′E 6669 2443 長崎、鹿児島西部 II 2 131°00′E 6670 2444 福岡・佐賀・熊本・大分・宮崎、鹿児島...
続きを読む

石川県：土砂災害（特別）警戒区域＋ＣＳ立体図

石川県の土砂災害（特別）警戒区域まとめ 2025 年 5 月 16 日 時点 石川県砂防課が公表する 土砂災害（特別）警戒区域指定状況 （PDF）を整理し、市町別に 土石流・地すべり・急傾斜地崩壊 の警戒区域（イエロー）と 特別警戒区域（レッド）を集計しました。 市町 警戒区域 (イエロー) 特別警戒区域 (レッド) 合計 土石流 地すべり 急傾斜 土石流 地すべり 急傾斜 輪島市 374 217 375 313 0 372 966 金沢市 230 82 247 170 0 235 559 能登町 240 68 218 210 0 209 526 七尾市 192 45 201 165 0 201 438 珠洲市 91 76 164 83 0 163 331 小松市 181 13 66 162 0 64 260 白山市 139 20 91 73 0 84 250 志賀町 130 15 100 106 0 99 245 穴水町 93 12 115 79 0 113 220 加賀市 78 21 101 65 0 96 200 津幡町 78 48 57 56 0 56 183 中能登町 130 4 36 99 0 35 170 能美市 64 1 41 59 0 40 106 宝達志水町 33 30 53 26 0 52 116 かほく市 12 2 30 5 0 30 44 川北町 0 0 0 0 0 0 0 野々市市 0 0 0 0 0 0 0 内灘町 0 0 5 0 0 5 5 合計 2065 654 1900 1671 0 1854 4619 ＣＳ立体図＋土砂災害警戒区域 MAP を開く 関連資料・データ ▸ 土砂災害（特別）警戒区域指定状況 PDF（R7.5.16） ▸ 石川県土砂災害（特別）警戒区域図 オープンデータ（SHP・GeoJSON 等） ▸ CS 立体図を使った 地形判読マニュアル（林野庁） ハザードマップの入手方法 お住まいの 市役所・町役場の防災担当窓口 や自...
続きを読む

設計定数を求めるための代表N値について

代表 N 値（設計用採用値）の決め方 ― 要点まとめ 1．実務でよく尋ねられる「出典」の悩み（改訂） 発注者から「この整理手法の出典はどこか」と問われることが多い一方で、道路土工・道路橋示方書・地盤調査法など主要文献には 「平均値±1/2σで採用する」と明記した箇所はほとんど見当たりません。 しかし、 『報告書作成マニュアル［土質編］』（社団法人 全国地質調査業協会連合会） には、 ばらつきデータから設計用値を定める 三手法 （ a： 平均、 b： 最小／最大、 c： 平均と標準偏差による補正）が明記されています。 特に c) として次式の例が示され、 設計値 ＝ 平均値 − (1/2) × 標準偏差 （平均値を標準偏差で補正して設計値とする） と記載されています（第2章 2.7.3、p.74）。 したがって、本稿で用いる 「平均 − 1/2σ」 という整理は単なる経験則ではなく、 公的マニュアルに裏づけのある実務手法 です。 なお同マニュアルは、標準偏差の重み係数（1/2 など）の選び方には議論があること、すなわち k の選定はデータ量・ばらつき・重要度等に応じて調整する余地があることも併記しています。 出典： 『報告書作成マニュアル［土質編］』第2章 2.7.3「設計用土質定数の設定方法」p.74（a〜c の三手法と「平均 − (1/2) 標準偏差」の例示） 2．建築基礎構造設計規準・同解説（1974）の記述 地盤の場合、土質試験などによって得られる諸数値にかなりのばらつきがあるのがむしろ普通であり、このようにばらつきのある測定値の取り扱いについては、ばらつきの原因とその程度、基礎底面と各試料採取位置との相対的な関係など、複雑な要素が絡むため一概には言えない。しかし通常、1地点で同一土質層について得られた測定値は、その平均値を用いて差し支えないとする。そのうえで、とくにばらつきが著しい場合などは、平均値よりさらに安全側を見込んで 設計用採用値 ＝ （平均値） ± （標準偏差）／2 とする。ここで複号 ± は、算定結果が安全側となる符号をとる...
続きを読む