水位流量観測所を地図上で探したい ～QGIS編～

はじめに

iRICで河床変動計算などの解析を実施する際、境界条件として対象河川の水位や流量のデータが必要になりますよね。その際、取得先の候補と言えば水文水質データベースが挙げられると思います。

任意の観測所の水位や流量を取得できる便利なサイトなのですが、個人的には使い勝手を改善したいことがあります。それは、

　
$\huge{観測所の位置がわかりにくい...}$
　

ということです。実際に北海道の石狩川の水位流量観測所を検索すると、以下の結果が得られます。

土地勘のある方なら分かるかもしれませんが、普通の人は所在地だけだと、どこにある観測所なのか分かりにくいです。

そのため、「これらの観測所を地図上で可視化しよう」というのが本記事の主旨です。

必要となるツール等

こんな記事に興味を持つ方は、十中八九プログラマーではなく土木系の学生 or 建設コンサルの方だと思うので、「Python環境を構築しましょう」みたいなことを言うつもりはないです。安心してください。

本記事では、以下をご準備頂ければ大丈夫です。

• QGIS 3.40
• エクセル
• 複数回のコピペ作業をやり抜く強靭な精神力

手順

おおまかな手順は以下の通りです。

• 可視化したいデータを準備する
• データをポイントに変換する

では、早速やってみましょう。

可視化したいデータを準備する

ここが最も簡単で、最も大変な作業です。まず前提として、水文水質データベースの「利用上の注意」には以下の記載があります。

当ホームページは一般を対象としており、通常のブラウザで閲覧することを前提に情報を掲載しております。ツール等による、自動的なデータ収集等はサーバに負荷がかかり、情報提供できなくなる恐れがありますので原則としてご遠慮ください。ご理解・ご協力お願いします。

つまり、Seleniumを使った自動化やスクレイピングによるデータ取得は御法度です。データが欲しくば、己の手で切り開くしか道はありません。

ということで、可視化したい水位流量観測所の情報をコピペでエクセルに移す作業が必要になります。そのため、可視化対象は最小限にしましょう。間違っても「全国の水位流量観測所を可視化してやる！」とか考えないでください。普通に心が折れます。

では、頑張ってやっていきましょう...

水文水質データベースにアクセスし、「観測所諸元からの検索」をクリックします。

今回は静岡県内の水位流量観測所を可視化してみようと思うので、「観測項目」を水位流量、「都道府県」を静岡県にして検索します。すると、44件のデータが見つかりました。（つまり44回コピペ作業をする気力が求められています）

各観測所名はリンクになっているので、リンクをクリックして観測所情報の緯度経度よりも上の情報をコピーして、エクセルに貼り付けましょう。

貼り付ける際、項目名の列は一番最初以外要らないので、消しておきましょう。

作業していると時折、座標情報が書いてない観測所に遭遇しますが、こういう場合どうしようもないので無視します。

...

さて、数多のコピペ作業を終えたら、エクセルは以下の様になっていると思います。

1～7行目を全て選択した状態で、エクセルの別シートに行/列を入れ替えて貼り付けましょう。

貼り付けた直後は、貼り付けたセルが全て選択状態になっていると思うので、「セルを結合して中央揃え」を押してセルの結合を解除しておきます。

あとは並べ替え機能を実行して余計な空白行をなくしつつ、観測所記号のセルの書式設定を数値に変えます。（この辺はただのエクセル作業なので細かくは説明しないです）

ここまで出来たら、「CSV（コンマ区切り）」形式で保存しましょう。保存をミスるとマジで心が折れると思うので心配な方は先に「Excelブック(*.xlsx)」形式で保存しておくと安心です。

データをポイントに変換する

ここからはQGISを使って作業します。

QGISのダウンロード方法などは、本記事では紹介しませんが、以下を参考にすれば良いと思います。

では、早速QGISを開いて、csvファイルを読み込みます。データソースマネージャを開きます。

以下の手順で作成したcsvファイルを読み込みます。

csvファイルを読み込んだら、ツールボタンの中からPythonマークを探して押します。

Pythonコンソールが出てくるので、「エディタの表示」アイコンをクリックします。

Pythonコードの入力欄が表示されるので、以下のコードをコピペで貼り付けます。

# -*- coding: utf-8 -*-

import re
from qgis.core import (
    QgsVectorLayer,
    QgsField,
    QgsFeature,
    QgsPointXY,
    QgsProject,
    QgsFields,
    QgsWkbTypes,
    QgsGeometry,  # QgsGeometryをインポート
)
from PyQt5.QtCore import QVariant

# --- 設定 ---
# QGISに読み込まれているソースレイヤ（テーブル）の名前
source_layer_name = "shizuoka_kansokujo"

# 作成する新しいポイントレイヤの名前
output_layer_name = "静岡県観測所_ポイント"

# 緯度経度情報が含まれているフィールド名
dms_field_name = "緯度経度"


# --- 関数定義 ---
def dms_to_dd(dms_str):
    """
    度分秒形式の緯度経度文字列を10進数に変換する
    例: "北緯 35度04分16秒 東経 138度22分12秒" -> (35.07111, 138.37000)
    """
    if not isinstance(dms_str, str):
        return None, None
    try:
        # 正規表現で緯度と経度の度分秒を抽出
        pattern = re.compile(
            r"北緯\s*(\d+)度\s*(\d+)分\s*(\d+)秒\s*東経\s*(\d+)度\s*(\d+)分\s*(\d+)秒"
        )
        match = pattern.match(dms_str.strip())
        if not match:
            return None, None

        lat_d, lat_m, lat_s, lon_d, lon_m, lon_s = map(float, match.groups())

        # 10進数に変換
        lat_dd = lat_d + (lat_m / 60) + (lat_s / 3600)
        lon_dd = lon_d + (lon_m / 60) + (lon_s / 3600)

        return lat_dd, lon_dd
    except (ValueError, TypeError):
        return None, None


# --- メイン処理 ---

# 1. QGISプロジェクトからソースレイヤを取得
layers = QgsProject.instance().mapLayersByName(source_layer_name)
if not layers:
    print(f'エラー: レイヤ "{source_layer_name}" がプロジェクト内に見つかりません。')
else:
    source_layer = layers[0]

    # 2. メモリ上に新しいポイントレイヤを作成 (CRSはWGS 84)
    vl = QgsVectorLayer("Point?crs=EPSG:4326", output_layer_name, "memory")
    pr = vl.dataProvider()

    # 3. ソースレイヤのフィールド定義をコピーし、新しいフィールドを追加
    source_fields = source_layer.fields()
    new_fields = QgsFields()
    for field in source_fields:
        new_fields.append(field)
    pr.addAttributes(new_fields)
    vl.updateFields()

    # 4. ソースレイヤの各フィーチャを処理
    dms_field_index = source_layer.fields().indexOf(dms_field_name)
    if dms_field_index == -1:
        print(
            f'エラー: フィールド "{dms_field_name}" がレイヤ "{source_layer_name}" に見つかりません。'
        )
    else:
        for source_feature in source_layer.getFeatures():
            # 緯度経度文字列を取得
            dms_string = source_feature[dms_field_index]

            # 10進数に変換
            lat, lon = dms_to_dd(dms_string)

            if lat is not None and lon is not None:
                # 新しいフィーチャ（ポイント）を作成
                new_feature = QgsFeature()

                # ジオメトリ（座標）を設定
                point = QgsPointXY(lon, lat)
                # QgsPointXYからQgsGeometryを作成して設定する
                new_feature.setGeometry(QgsGeometry.fromPointXY(point))

                # 元レイヤの属性をコピーし、新しい属性を追加
                attributes = source_feature.attributes()
                new_feature.setAttributes(attributes)

                # レイヤにフィーチャを追加
                pr.addFeature(new_feature)

        # 5. レイヤの更新とプロジェクトへの追加
        vl.updateExtents()
        QgsProject.instance().addMapLayer(vl)

        print(f'レイヤ "{output_layer_name}" の作成が完了しました。')

貼り付けたプログラムは一部修正が必要です。とはいえ、source_layer_nameに読み込んだレイヤの名前、output_layer_name にこれから作るポイントレイヤの名前を記入するだけです。後は何も変更しなくていいです。

レイヤ名を変更したら実行ボタンを押しましょう。csvファイル内の観測所にポイントが配置されるはずです。

ベースマップと重ねたり、ポイントのスタイルを変更すると、良い感じに可視化することができます。（QGISの基本的な操作なので、本記事では具体の手順は割愛します）

地理院地図の上に赤丸で観測所を重ねた図

これで可視化した範囲の水位流量観測所はすぐに確認できるようになりましたね。

おわりに

今回はQGIS上で水位流量観測所を可視化する方法を紹介しました。一度作ってしまえば（観測所の位置が変わらない限り）ずっと使えるので、観測所の位置を確認するのがだいぶ楽になるのではないでしょうか？

わざわざタイトルに「QGIS編」と入れたということは...