CUI で GIS (5)

　当ブログ、長い夏休みを経て久方ぶりの更新です。

　今回は「CUI で GIS」の第５弾で、小流域分割をテーマに取り上げます。
　仕上がりのイメージは下図のようです。小流域分割を表す紙の地図からこれの多角形の GIS データに作成しました。

　大まかな作業の流れは

スキャンした紙の地図に地上基準点（GCP：Ground Control Point）の座標を与えて幾何補正した地図（GeoTif 形式）に変換する
幾何補正後の地図上の小流域の多角形の頂点群のデータを作成する
上の頂点を結んで小流域毎のデータを作成する

となります。

地図のディジタライズ

　地図のディジタライズは QGIS のジオリファレンサを用いて行いました。ジオリファレンサの使い方は QGIS のオンラインマニュアルをご覧下さい。
　GCP の緯度・経度は地理院地図のセンターマーク（＋）をその点に移動させて表示される URL から取得しました。

頂点の座標の取得

　これも QGIS 上での作業です。
　属性として頂点番号 pid を持つ GeoJSON 形式の nodes.json というファイルにしました。

小流域毎のデータ

　隣接する小流域は頂点を共有しますが、手作業でポリゴンを拾うと、これが正確にできないので「キーーッ！」となります。そこで、今回は上の作業を挟んで Python でこの処理を行うことにしました。モジュール名は watersheds.py としました。モジュール内で ogr および osr というライブラリを用いているので QGIS 同梱の Python を使います。このための環境設定を watersheds.bat で行っています（8 行目の Python39 は適宜変更して下さい）。

: watersheds.bat

@echo off
setlocal

set OSGEO4W_ROOT=C:\OSGEO4W
set PATH=%OSGEO4W_ROOT%\bin;%PATH%
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python39
call %OSGEO4W_ROOT%\etc\ini\gdal.bat
call %OSGEO4W_ROOT%\etc\ini\python3.bat

if exist watershed.json (del watershed.json) 
python watersheds.py

# watersheds.py 

from osgeo import ogr
from osgeo import osr

# https://gdal.org/en/latest/api/python/osgeo.ogr.html
# https://gdal.org/en/latest/api/python/osgeo.osr.html

watersheds = (
    {"name": "A", "nodes": [4, 5, 6, 7, 15, 16, 17, 18,
                            8, 36, 35, 34, 4]},
    {"name": "B", "nodes": [1, 3, 2, 4, 34, 35, 36, 8, 18,
                            9, 10, 11, 12, 13, 14, 1]},
    {"name": "C", "nodes": [2, 19, 20, 21, 22, 23, 24, 16,
                            15, 7, 6, 5, 4, 2]},
    {"name": "D", "nodes": [21, 25, 26, 27, 28, 29, 24, 23, 22, 21]},
    {"name": "E", "nodes": [26, 30, 31, 32, 33, 27, 26]}
)

driver = ogr.GetDriverByName("GeoJSON")
dataSource = driver.Open("nodes.json")
layer = dataSource.GetLayer("nodes")
feature_count = layer.GetFeatureCount()

pts = []
for _ in range(feature_count):
    feature  = layer.GetNextFeature()
    geometry = feature.GetGeometryRef()
    x = geometry.GetX()
    y = geometry.GetY()
    pts.append([x, y])

layer = None  
dataSource = None

srs = osr.SpatialReference()
srs.ImportFromEPSG(4612)

driver = ogr.GetDriverByName("GeoJSON")
dataSource = driver.CreateDataSource("watersheds.json")
layer = dataSource.CreateLayer("watersheds", srs,
                        geom_type=ogr.wkbPolygon)

nameField = ogr.FieldDefn("name", ogr.OFTString)
layer.CreateField(nameField)

featureDefn = layer.GetLayerDefn()
for watershed in watersheds:
    ring = ogr.Geometry(ogr.wkbLinearRing)
    for node in watershed["nodes"]:
        pt = pts[node - 1]
        ring.AddPoint_2D(pt[0], pt[1])
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)
    feature = ogr.Feature(featureDefn)
    feature.SetGeometry(poly)
    feature.SetField("name", watershed["name"])
    layer.CreateFeature(feature)
    feature = None

layer = None
dataSource = None

　watersheds.py について簡単に説明すると、

9 ～ 18 行目は各小流域の多角形を構成する頂点を QGIS で表示した nodes の頂点番号を確認しながら、反時計回りで拾った配列と流域名の辞書の配列 watersheds に記憶します
20 ～ 34 行目は nodes.json を読み込んで全ての頂点の緯度・経度の配列の配列を pts に収めます
36 ～ 37 行目は EPSG コードが 4612 である JGD2000 の空間参照系のオブジェクトを srs に設定します
39～45 行目ではレイヤー名、地物のタイプ、空間参照系を指定し、name 属性を持つ waterhseds.json を新規に登録します
47～59 行目では各小流域について辞書の配列 watersheds に基づいてポリゴンを作成し、name 属性を与えます。つまり waterhseds.json に各小流域の情報を書き込んでいます

　以前の記事で「定型の処理が多いので実際に記述するコードはわずか」的なことを述べましたが、ここでも他所で書いたコードに手を加えて作成しています。なので汎用性を考えずに気軽に使い捨てのコードを書くことができます。
　そして毎回強調しますが、頂点の位置がづれていた、小流域を構成する頂点が間違えていたなどが発生しても簡単に直すことができます（実際にそうなった）。