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

はじめに

土石流区域設定支援システムのデータベース（dem.mdb）から、縦断情報と横断情報を読み取り、QGIS上に測線を作成するPythonスクリプトを作成しました。

当初はMDBを直接開く部分でうまくいかなかったため、次の手順で処理していました。

これまでの方法
MDBをAccessで開く → Excelの各シートへコピー → CSV出力 → QGISでCSVを選択 → Pythonを実行

今回、Microsoft AccessのODBC接続方法を見直したことで、QGISからMDBを直接読み込み、1回の操作で2つのレイヤを一括作成できるようになりました。

新しい方法
QGISでdem.mdbを選択 → 「流下方向」と「横断測線」を一括作成

MDB内の独自形式である「カンマ区切り座標」と「mm単位」を自動変換し、以下の2レイヤを生成します。

  1. 流下方向（縦断）：区間勾配を計算し、勾配に応じて色分けしたライン
  2. 横断測線：各横断グループの始点と終点を結んだライン
QGISで作成した縦断・横断測線

1. 更新版：MDBを直接読み込んで一括作成

この更新版では、MDB内の次の2テーブルを直接読み込みます。

  • 縦断情報：GISKey、ID、ILData
  • 横断情報：GISKey、ID、ILData

ILDataには、標高・X座標・Y座標がカンマ区切りかつmm単位で記録されています。スクリプトはこれをm単位へ変換してQGISのラインを作成します。

一括処理の内容

  • MDBファイルを1回選択するだけで、縦断・横断をまとめて処理
  • GISKeyごとにデータを分離し、異なる箇所同士の誤接続を防止
  • 縦断はID順に接続し、距離・標高差・勾配角度を計算
  • 横断はIDの10000単位のグループごとに、最小IDと最大IDを接続
  • 作成した2レイヤをMDB名のレイヤグループに格納
  • 処理件数とスキップ件数を完了画面に表示

作成される属性

レイヤ 主な属性
流下方向 GISKey、始点ID、終点ID、勾配角、絶対勾配角、水平距離、標高差、始点・終点座標
横断測線 連番、GISKey、グループID、始点ID、終点ID、点数、始点・終点座標

実行環境

  • Windows版QGIS
  • QGISと同じビット数のMicrosoft Access ODBCドライバー
  • MDB内に「縦断情報」「横断情報」テーブルが存在すること

通常の64bit版QGISを使用している場合は、64bit版のMicrosoft Access Database Engineが必要です。

QGISでの実行方法

  1. 下記のコードをUTF-8形式のPythonファイルとして保存します。
  2. QGISで「プラグイン」→「Pythonコンソール」を開きます。
  3. 「エディタを表示」を押して、保存したPythonファイルを開きます。
  4. 実行ボタンを押し、対象のdem.mdbを選択します。
  5. レイヤパネルに「流下方向」と「横断測線」が追加されます。

Pythonコンソールから直接実行する場合は、次のように入力します。

exec(open(r"C:\保存先\QGIS_MDB_batch_v2.py", encoding="utf-8").read())

MDB接続でつまずいた原因

MDBを直接開く処理では、データ解析より先にODBC接続でつまずきやすい点に注意が必要です。今回も当初は「MDBを開けませんでした」と表示されましたが、Accessドライバーそのものは認識されていました。

原因は、接続文字列のMDBパスを波括弧で囲んでいたことでした。

うまくいかなかった例
DBQ={C:\フォルダ\dem.mdb}
正しい例
DBQ=C:\フォルダ\dem.mdb

Accessドライバーによっては、波括弧までファイル名として解釈され、「ファイル名が正しくありません」となる場合があります。また、日本語や空白を含むOneDrive上のパスに備え、更新版ではWindowsの短縮パスでも再試行するようにしています。

更新版スクリプトコード

# -*- coding: utf-8 -*-
"""
QGIS Pythonコンソール用
Microsoft Access MDBから「縦断情報」「横断情報」を直接読み込み、
  1) 流下方向（区間勾配で色分け）
  2) 横断測線
を一括作成します。

前提:
- Windows版QGIS
- QGISと同じビット数のMicrosoft Access ODBCドライバーがインストール済み
- MDB内に [縦断情報] / [横断情報] テーブルがあり、
  GISKey, ID, ILData の各列が存在すること
"""

import math
import os
import uuid
from collections import defaultdict

from qgis.core import (
    QgsFeature,
    QgsField,
    QgsGeometry,
    QgsGraduatedSymbolRenderer,
    QgsPointXY,
    QgsProject,
    QgsRendererRange,
    QgsSymbol,
    QgsVectorLayer,
)
from qgis.PyQt.QtCore import QVariant
from qgis.PyQt.QtGui import QColor
from qgis.PyQt.QtSql import QSqlDatabase, QSqlQuery
from qgis.PyQt.QtWidgets import QFileDialog, QMessageBox
from qgis.utils import iface


# 岐阜県（平面直角座標系第7系）の既定値。
# QGISプロジェクトが投影座標系なら、プロジェクトCRSを優先します。
DEFAULT_CRS = "EPSG:6675"
VERTICAL_TABLE = "縦断情報"
TRANSVERSE_TABLE = "横断情報"


def _output_crs_authid():
    """メートル座標として使える出力CRSを返す。"""
    project_crs = QgsProject.instance().crs()
    if project_crs.isValid() and not project_crs.isGeographic():
        authid = project_crs.authid()
        if authid:
            return authid
    return DEFAULT_CRS


def _installed_access_driver_names():
    """Windowsレジストリと一般的な名称からAccess ODBCドライバー候補を返す。"""
    names = []

    try:
        import winreg

        registry_paths = [
            r"SOFTWARE\ODBC\ODBCINST.INI\ODBC Drivers",
            r"SOFTWARE\WOW6432Node\ODBC\ODBCINST.INI\ODBC Drivers",
        ]
        for root in (winreg.HKEY_LOCAL_MACHINE, winreg.HKEY_CURRENT_USER):
            for reg_path in registry_paths:
                try:
                    with winreg.OpenKey(root, reg_path) as key:
                        index = 0
                        while True:
                            try:
                                driver_name, value, _ = winreg.EnumValue(key, index)
                                index += 1
                                text = driver_name.lower()
                                if (
                                    "access" in text
                                    and ".mdb" in text
                                    and str(value).lower() == "installed"
                                ):
                                    names.append(driver_name)
                            except OSError:
                                break
                except OSError:
                    pass
    except ImportError:
        pass

    # 通常の日本語Windowsでも、ドライバー名はこの英語名で登録されます。
    names.extend([
        "Microsoft Access Driver (*.mdb, *.accdb)",
        "Microsoft Access Driver (*.mdb)",
    ])

    # 重複除去（順序保持）
    return list(dict.fromkeys(names))


def _open_mdb(mdb_path):
    """QODBCでMDBを開き、(database, connection_name) を返す。"""
    if "QODBC" not in QSqlDatabase.drivers():
        raise RuntimeError(
            "QGISのQt SQLにQODBCドライバーがありません。\n"
            "Windows版QGISのインストール構成を確認してください。"
        )

    errors = []

    # Access ODBC の DBQ はファイルパスを波括弧で囲まない。
    # DBQ={C:\...\file.mdb} とすると、Access ドライバーが波括弧まで
    # ファイル名として解釈し「ファイル名が正しくありません」になることがある。
    native_path = os.path.normpath(os.path.abspath(mdb_path))
    if not os.path.isfile(native_path):
        raise RuntimeError("MDBファイルが見つかりません。\n{}".format(native_path))

    # 通常パスに加え、取得できる場合はWindowsの短いパスも試す。
    # 日本語・空白を含むOneDrive上のパスで古いAccessドライバーが失敗する場合の保険。
    path_candidates = [native_path]
    try:
        import ctypes

        buffer_size = 32768
        buffer = ctypes.create_unicode_buffer(buffer_size)
        result = ctypes.windll.kernel32.GetShortPathNameW(
            native_path, buffer, buffer_size
        )
        if result and buffer.value and buffer.value not in path_candidates:
            path_candidates.append(buffer.value)
    except Exception:
        pass

    for driver_name in _installed_access_driver_names():
        for dbq_path in path_candidates:
            connection_name = "mdb_" + uuid.uuid4().hex
            db = QSqlDatabase.addDatabase("QODBC", connection_name)
            connection_string = "DRIVER={{{}}};DBQ={};READONLY=TRUE;".format(
                driver_name, dbq_path
            )
            db.setDatabaseName(connection_string)

            if db.open():
                return db, connection_name

            errors.append(
                "{} / {}: {}".format(
                    driver_name, dbq_path, db.lastError().text()
                )
            )
            db.close()
            del db
            QSqlDatabase.removeDatabase(connection_name)

    detail = "\n".join(errors)
    raise RuntimeError(
        "MDBを開けませんでした。\n"
        "QGISと同じビット数のMicrosoft Access ODBCドライバーが必要です。\n\n"
        + detail
    )


def _read_access_table(db, table_name):
    """GISKey, ID, ILDataをMDBテーブルから読み込む。"""
    available = db.tables()
    if table_name not in available:
        raise RuntimeError(
            "MDB内にテーブル「{}」がありません。\n\n存在するテーブル:\n{}".format(
                table_name, "\n".join(available)
            )
        )

    sql = (
        "SELECT [GISKey], [ID], [ILData] "
        "FROM [{}] ORDER BY [GISKey], [ID]".format(table_name)
    )
    query = QSqlQuery(db)
    exec_method = getattr(query, "exec", None) or query.exec_
    if not exec_method(sql):
        raise RuntimeError(
            "テーブル「{}」の読み込みに失敗しました。\n{}".format(
                table_name, query.lastError().text()
            )
        )

    rows = []
    skipped = 0
    while query.next():
        try:
            gis_key = int(query.value(0))
            point_id = int(query.value(1))
            il_data = str(query.value(2)).strip()
            if not il_data:
                skipped += 1
                continue
            rows.append((gis_key, point_id, il_data))
        except (TypeError, ValueError):
            skipped += 1

    return rows, skipped


def _parse_ildata(il_data):
    """ILDataから Z, X, Y（mm→m）を取り出す。"""
    try:
        parts = [part.strip() for part in str(il_data).split(",")]
        if len(parts) < 4:
            return None
        z = float(parts[1]) / 1000.0
        x = float(parts[2]) / 1000.0
        y = float(parts[3]) / 1000.0
        return z, x, y
    except (TypeError, ValueError):
        return None


def _create_transverse_layer(rows, crs_authid):
    """横断情報から、各横断グループの最小ID～最大IDを結ぶ測線を作成。"""
    groups = {}
    invalid = 0

    # GISKeyもキーに含めるため、複数箇所入りMDBでも混線しない。
    for gis_key, point_id, il_data in rows:
        group_no = point_id // 10000
        if group_no <= 0:
            invalid += 1
            continue

        parsed = _parse_ildata(il_data)
        if parsed is None:
            invalid += 1
            continue

        _, x, y = parsed
        key = (gis_key, group_no)
        item = (point_id, x, y)

        if key not in groups:
            groups[key] = {"min": item, "max": item, "count": 1}
        else:
            groups[key]["count"] += 1
            if point_id < groups[key]["min"][0]:
                groups[key]["min"] = item
            if point_id > groups[key]["max"][0]:
                groups[key]["max"] = item

    layer = QgsVectorLayer(
        "LineString?crs={}".format(crs_authid), "横断測線", "memory"
    )
    provider = layer.dataProvider()
    provider.addAttributes([
        QgsField("Seq_No", QVariant.Int),
        QgsField("GISKey", QVariant.LongLong),
        QgsField("Group_ID", QVariant.Int),
        QgsField("Start_ID", QVariant.Int),
        QgsField("End_ID", QVariant.Int),
        QgsField("Pt_Count", QVariant.Int),
        QgsField("Start_X", QVariant.Double),
        QgsField("Start_Y", QVariant.Double),
        QgsField("End_X", QVariant.Double),
        QgsField("End_Y", QVariant.Double),
    ])
    layer.updateFields()

    features = []
    for seq_no, key in enumerate(sorted(groups.keys())):
        gis_key, group_no = key
        group = groups[key]
        start_id, start_x, start_y = group["min"]
        end_id, end_x, end_y = group["max"]

        feature = QgsFeature(layer.fields())
        feature.setGeometry(
            QgsGeometry.fromPolylineXY([
                QgsPointXY(start_x, start_y),
                QgsPointXY(end_x, end_y),
            ])
        )
        feature.setAttributes([
            seq_no,
            gis_key,
            group_no * 10000,
            start_id,
            end_id,
            group["count"],
            start_x,
            start_y,
            end_x,
            end_y,
        ])
        features.append(feature)

    provider.addFeatures(features)
    layer.updateExtents()
    layer.setCustomProperty("source_table", TRANSVERSE_TABLE)
    return layer, len(features), invalid


def _apply_slope_renderer(line_layer):
    """絶対勾配角による既存の色区分を適用。"""
    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, 999.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("slope_abs", ranges)
    renderer.setMode(QgsGraduatedSymbolRenderer.Custom)
    line_layer.setRenderer(renderer)


def _create_vertical_layer(rows, crs_authid):
    """縦断情報をGISKey別にID順で結び、区間勾配を計算。"""
    point_groups = defaultdict(list)
    invalid = 0

    for gis_key, point_id, il_data in rows:
        parsed = _parse_ildata(il_data)
        if parsed is None:
            invalid += 1
            continue
        z, x, y = parsed
        point_groups[gis_key].append((point_id, x, y, z))

    layer = QgsVectorLayer(
        "LineString?crs={}".format(crs_authid), "流下方向", "memory"
    )
    provider = layer.dataProvider()
    provider.addAttributes([
        QgsField("GISKey", QVariant.LongLong),
        QgsField("Start_ID", QVariant.Int),
        QgsField("End_ID", QVariant.Int),
        QgsField("slope_angle", QVariant.Double),
        QgsField("slope_abs", QVariant.Double),
        QgsField("distance_m", QVariant.Double),
        QgsField("dz_m", QVariant.Double),
        QgsField("Start_X", QVariant.Double),
        QgsField("Start_Y", QVariant.Double),
        QgsField("End_X", QVariant.Double),
        QgsField("End_Y", QVariant.Double),
    ])
    layer.updateFields()

    features = []
    for gis_key in sorted(point_groups.keys()):
        points = sorted(point_groups[gis_key], key=lambda value: value[0])
        for index in range(len(points) - 1):
            start_id, x1, y1, z1 = points[index]
            end_id, x2, y2, z2 = points[index + 1]

            dx = x2 - x1
            dy = y2 - y1
            dz = z2 - z1
            distance = math.hypot(dx, dy)

            if distance > 0:
                angle = math.degrees(math.atan2(dz, distance))
            elif dz > 0:
                angle = 90.0
            elif dz < 0:
                angle = -90.0
            else:
                angle = 0.0

            feature = QgsFeature(layer.fields())
            feature.setGeometry(
                QgsGeometry.fromPolylineXY([
                    QgsPointXY(x1, y1),
                    QgsPointXY(x2, y2),
                ])
            )
            feature.setAttributes([
                gis_key,
                start_id,
                end_id,
                angle,
                abs(angle),
                distance,
                dz,
                x1,
                y1,
                x2,
                y2,
            ])
            features.append(feature)

    provider.addFeatures(features)
    layer.updateExtents()
    layer.setCustomProperty("source_table", VERTICAL_TABLE)
    _apply_slope_renderer(layer)
    return layer, len(features), invalid


def create_layers_from_mdb():
    mdb_path, _ = QFileDialog.getOpenFileName(
        None,
        "MDBファイルを選択してください",
        "",
        "Microsoft Access Database (*.mdb *.accdb);;All Files (*)",
    )
    if not mdb_path:
        return

    db = None
    connection_name = None

    try:
        db, connection_name = _open_mdb(mdb_path)

        vertical_rows, vertical_read_skipped = _read_access_table(
            db, VERTICAL_TABLE
        )
        transverse_rows, transverse_read_skipped = _read_access_table(
            db, TRANSVERSE_TABLE
        )

        if not vertical_rows:
            raise RuntimeError("テーブル「縦断情報」に有効なデータがありません。")
        if not transverse_rows:
            raise RuntimeError("テーブル「横断情報」に有効なデータがありません。")

        crs_authid = _output_crs_authid()
        vertical_layer, vertical_count, vertical_invalid = _create_vertical_layer(
            vertical_rows, crs_authid
        )
        transverse_layer, transverse_count, transverse_invalid = (
            _create_transverse_layer(transverse_rows, crs_authid)
        )

        # レイヤパネル内でMDBごとにまとめる。
        project = QgsProject.instance()
        group_name = "MDB解析_{}".format(os.path.splitext(os.path.basename(mdb_path))[0])
        root = project.layerTreeRoot()
        group = root.addGroup(group_name)

        project.addMapLayer(vertical_layer, False)
        group.addLayer(vertical_layer)
        project.addMapLayer(transverse_layer, False)
        group.addLayer(transverse_layer)

        iface.mapCanvas().refresh()

        QMessageBox.information(
            None,
            "完了",
            "MDBを直接読み込み、2レイヤを作成しました。\n\n"
            "・流下方向: {:,} 区間\n"
            "・横断測線: {:,} 本\n"
            "・出力CRS: {}\n\n"
            "読み込み時スキップ: 縦断 {:,} / 横断 {:,}\n"
            "ILData解析スキップ: 縦断 {:,} / 横断 {:,}".format(
                vertical_count,
                transverse_count,
                crs_authid,
                vertical_read_skipped,
                transverse_read_skipped,
                vertical_invalid,
                transverse_invalid,
            ),
        )

    except Exception as error:
        QMessageBox.critical(None, "エラー", str(error))

    finally:
        if db is not None:
            db.close()
            del db
        if connection_name:
            QSqlDatabase.removeDatabase(connection_name)


create_layers_from_mdb()
注意：作成されるレイヤはメモリレイヤです。QGISを閉じる前に、GeoPackageやShapefileなどへエクスポートして保存してください。

2. 旧方式：CSVから流下方向（縦断）を作成

以下は、MDBを直接開けない環境で使用していたCSV経由の旧方式です。ODBCドライバーを導入できないパソコンでは、現在も代替手段として利用できます。

縦断データの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()

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

3. 旧方式：CSVから横断測線を作成

横断測線の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()

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

おわりに

当初はMDBを開くところでつまずき、ExcelとCSVを経由する方法で運用していました。この方法でも処理は可能でしたが、コピー、シート選択、CSV保存、ファイル選択といった作業が必要で、件数が多いほど手間や選択ミスが生じやすくなります。

今回の更新により、MDBを1回選択するだけで、縦断情報と横断情報を直接読み込み、2レイヤを一括作成できるようになりました。CSV作成が不要になったことで、作業時間の短縮だけでなく、古いCSVを使用するミスやコピー漏れの防止にもつながります。

なお、ODBCドライバーが利用できない環境に備え、従来のCSV版スクリプトも代替手段として残しています。

コメント

コメントを投稿

このブログの人気の投稿

IJCAD（AutoCAD互換CAD） コマンド早見表

AutoCAD LT 2025 コマンド早見表（完全版） AutoCAD LT 2025 主要コマンド早見表【日本語備考付き・完全版】 AutoCAD LT で頻出するコマンドを アルファベット順 にまとめました。 備考欄には「実務で役立つワンポイント」を一行で追加。★印はブロック・外部参照まわりで特に押さえておきたいコマンドです。 略号 正式コマンド 機能概要 備考 A ARC 円弧を作成 半径・３点・開始終了角で指定可 AA AREA 面積と周長を計算 ポリライン選択で複数領域合算 ADC ADCENTER コンテンツセンター ブロックや画像をドラッグ挿入 AL ALIGN 2D/3D で位置合わせ 3 点指定でスケール同時調整 AR ARRAY 配列複写 楕円・パス配列も可能 AS ASSISTANTOPEN Autodesk Assistant Web ヘルプ検索に便利 ATI ATTIPEDIT★ ブロック属性文字を直接編集 図面上でダブルクリックと同等 ATT ATTDEF★ ブロック属性を定義 タグ・既定値・プロンプトを設定 ATE ATTEDIT★ 挿入済み属性を編集 EATTEDIT GUI と同等 B BLOCK★ ブロック定義を作成 書き出し前に 基点 を忘れず指定 BC BCLOSE★ ブロックエディタを閉じる SAVEAS 不可状態から復帰 BE BEDIT★ ブロックエディタを開く 動的ブロック対応の本格編集 BH HATCH ハッチング 「固有原点」でパターンズレ防止 BO BOUNDARY 領域/ポリライン生成 内部島検出で高速ポリライン化 BR BREAK オブジェクトを分割 F オプションで１点除去 BS BSAVE★ ブロック定義を保存 BEDIT 中の上書き保存 BVS BVSTATE★ 動的ブロック可視状態 仕上違いを１定義で管理 C CIRCLE 円を作成 2P/3P 補助で簡易作図 CBAR CONSTRAIN...
続きを読む

設計定数を求めるための代表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 とする。ここで複号 ± は、算定結果が安全側となる符号をとる...
続きを読む

WEB地図（国土地理院など）のズームレベルについて

ズームレベルについて 0. ズームレベルとは？ Google Maps、国土地理院地図、OpenStreetMap (OSM)、Bing Maps などのWeb地図は、 ズームレベル (Zoom Level, ZL) で表示倍率を管理しています。 ZLが小さいと広域（世界全体）、ZLが大きいと詳細（街区レベル）になります。 これは世界中で共通の仕組み（Web Mercator投影 / EPSG:3857）に基づいています。 1. 計算式と一覧 ズームレベルごとの解像度（1ピクセルあたりの地上距離）は次の式で表されます： 解像度 [m/px] = (cos φ × 2πR) / (256 × 2^Z) φ = 緯度（ラジアン） R = 地球半径 6,378,137 m（WGS84） Z = ズームレベル 256 px = タイル1枚の大きさ 日本（北緯36度付近）の目安表： ZL 縮尺の目安 ZL 12 1:250,000 ZL 13 1:120,000 ZL 14 1:60,000 ZL 15 1:30,000 ZL 16 1:15,000 ZL 17 1:8,000 ZL 18 1:4,000 ZL 19 1:2,000 ZL 20 1:1,000 ZL 21 1:500 ※緯度によって多少変動しますが、日本の本州（北緯36°付近）では上記でほぼ問題ありません。 2. URLのどこにズームレベルが書かれているか Google Maps ：URL末尾の「◯◯z」がズームレベル 例： https://www.google.com/maps/@36.5844,136.6651, 11z 国土地理院地図 ：URLの #◯◯/緯度/経度/ の「◯◯」がズームレベル 例： https://maps.gsi.go.jp/vector/ #12 /36.6710/136.8362/ 3. WGS84について ズームレベルを扱うWeb地図は「WGS84」という地球基準座標系をもとにしています。 ただし表現方法に2種類あります： ✅ EPSG:4326 (WGS84) ：緯度・経度（度）。GNSSや測量で...
続きを読む