CUI で GIS

今回は GIS で主題図を作成するときに前処理を CUI で行うと何かと良いことがある、という話です。
そこで登場するのが ogr2ogr というコマンドユーティリティーです。QGIS のユーザーであれば ogr2ogr はインストール済です。というより、QGIS でベクターデータを処理するときの GUI は ogr2ogr のコマンドラインを生成するために設けられていることが殆どです。このコマンドラインをバッチファイルに記述することで冗長な作業が大幅に省けます。

事例

以降、「（国名、領土、領主が現在と異なる）明治時代の東南アジアの地図を作る」という課題を例にとって説明します。

作業手順

作業の手順は以下の通りです。今回は 2～5 の作業をバッチファイルで一気に処理します。

基図のデータのダウンロード
東南アジア域のデータの抽出
独立前の領地のポリゴンデータの作成
属性の変更
データの集約
QGIS による主題図の作成

環境変数の設定

ogr2ogr の具体的なインストール先を PATH 環境変数に加えます。また、ogr2ogr で用いる他の変数（今回未使用の分を含む）を設定します。これらはバッチファイルの先頭に記述する必要があります。
以下は ogr2ogr が C:\OSGeo4W\bin にインストールされている場合の例です。これに該当しない場合は実状にあわせて変数 OSGEO4W を書き換える必要があります。

@echo off
setlocal
set OSGEO4W_ROOT=C:\OSGEO4W
set PATH=%OSGEO4W_ROOT%\bin;%PATH%
call %OSGEO4W_ROOT%\etc\ini\gdal.bat

今回は「手順 4」で Python のプログラムを実行するので、さら以下を追加します。

call %OSGEO4W_ROOT%\etc\ini\python3.bat

手順 1　基図データのダウンロード

今回使用した基図はフリーの世界地図を提供している Natural Earth というサイトの国境データ「ne_110m_admin_0_countries.zip」です。
ogrinfo.exe でこのデータの仕様を確認します。このように ogr ユーティリティーは /vsizip/ を接頭すると zip ファイルの解凍が内部的に処理されます。

c:\OSGeo4W\bin\ogrinfo -so /vsizip/ne_110m_admin_0_countries.zip ne_110m_admin_0_countries

　これの出力から以下が確認できます。

データフォーマットは ESRI Shapefile
座標参照系は WGS 84
国名、国名を表示する軽度と緯度の属性名は NAME_JA、LABEL_X、LABEL_Y

手順 2　東南アジア域のデータの抽出

以下のコマンドで基図のデータから所望の東南アジア域のデータを southEastAsia.json というファイル名で抽出します。GeoJson はベクターデータのための JSON のフォーマットです。

If exist southEastAsia.json del southEastAsia.json
ogr2ogr ^
    -nln southEastAsia ^
    -lco COORDINATE_PRECISION=5 ^
    -sql "SELECT NAME_JA AS NAME, LABEL_X AS X, LABEL_Y AS Y FROM ne_110m_admin_0_countries" ^
    -clipsrc 45 -12 150 45 ^
    southEastAsia.json /vsizip/ne_110m_admin_0_countries.zip

ogr2ogr のオプションは以下で与えています。

レイヤー名は southEastAsia
国境線の座標は小数点以下 5 桁
属性は NAME_JA、LABEL_X、LABEL_Y だけを取出して、属性名を NAME、X、Y に変更
抽出する領域は南西端と北東端の緯度・経度を (-12, 45)、(45, 150) で指定

手順 3　独立前の領地のポリゴンデータの作成

以下のバッチファイルは次の操作を行います。

カンボジア、ラオス、ベトナムのデータからインドシナ連邦のポリゴンデータ indochina.json を作成
インド、バングラデシュのデータからインド帝国のポリゴンデータ india.json を作成

set UNION=SELECT ST_Union(ST_buffer(geometry, 0.0001)) AS geometry
set UNION=%UNION% FROM southEastAsia WHERE NAME IN

ogr2ogr ^
    -lco COORDINATE_PRECISION=5 ^
    -dialect sqlite ^
    -sql "%UNION% ('カンボジア', 'ラオス', 'ベトナム')" ^
    indochina.json southEastAsia.json

ogr2ogr ^
    -lco COORDINATE_PRECISION=5 ^
    -dialect sqlite ^
    -sql "%UNION% ('インド', 'バングラデシュ')" ^
    india.json southEastAsia.json

ポリゴンの合併（union）処理には SQLite の ST_Union 関数を用います。この関数は標準 SQLにはないので dialect（方言）オプションを与えています。

手順 4　属性の変更

Python のプログラムで属性を主題図用のものに変更します。

python edit_props.py

edit_props.py は以下の通りです。southEastAsia.json、indochina.json、india.json の国名を旧国名に変更し、新たな属性として統治国の頭文字（図中の凡例のカテゴリ）を追加します。後二者についてはラベルの位置も変更します。また元の southEastAsia.json からインドシナ連邦とインド帝国の構成国が除かれます。(Geo)Json は Python の辞書型あるいは辞書型のリストに変換可能なので、このような処理が比較的容易に行えます。

import json

def get_json(path):

    with open(path, encoding="utf-8") as f:
        return json.load(f)

def over_write(path, js):

    with open(path, "w", encoding="utf-8") as f:
        json.dump(js, f, ensure_ascii=False)

def edit_multi():

    skipes = ("カンボジア", "ラオス", "ベトナム", "インド", "バングラデシュ")
    countries = {
        "インドネシア": {"name_old": "オランダ領東インド", "ADMIN": "N"},
        "タイ王国": {"name_old": "シャム", "ADMIN": "O"},
        "ミャンマー": {"name_old": "ビルマ", "ADMIN": "E"},
        "パキスタン": {"name_old": "バルチスタン", "ADMIN": "E"},
        "アフガニスタン": {"name_old": "アフガニスタン", "ADMIN": "E"},
        "スリランカ": {"name_old": "セイロン", "ADMIN": "E"},
        "中華人民共和国": {"name_old": "清", "ADMIN": "O"},
        "中華民国": {"name_old": "台湾", "ADMIN": "J"},
        "フィリピン": {"name_old": "フィリピン", "ADMIN": "S"},
        "マレーシア": {"name_old": "マライ連邦", "ADMIN": "E"},
        "ネパール": {"name_old": "ネパール", "ADMIN": "O"}
    }

    js = get_json("southEastAsia.json")

    pop_indexes = []
    for idx, feat in enumerate(js["features"]):
        props = feat["properties"] 
        name  = props["NAME"]
        if name in skipes:
            pop_indexes.append(idx)
        elif name in countries:
            obj = countries[name]
            props["NAME"]  = obj["name_old"] # 旧国名
            props["ADMIN"] = obj["ADMIN"]    # 統治国の頭文字
        else: # 図中に国境のみ表示（国名非表示）
            props["NAME"]  = ""
            props["ADMIN"] = "O"

    # southEastAsia.json からカンボジア、ラオス、ベトナム
    # インド、バングラデシュのデータを消去
    for idx in reversed(pop_indexes):
        js["features"].pop(idx)

    over_write("southEastAsia.json", js)

def edit_single(path, NAME, X, Y, ADMIN):

    js = get_json(path)

    props = js["features"][0]["properties"] 
    props["NAME"]  = NAME
    props["X"]     = X # ラベルの表示位置（経度）
    props["Y"]     = Y # ラベルの表示位置（緯度）
    props["ADMIN"] = ADMIN

    over_write(path, js)

if __name__ == "__main__":

    edit_multi()
    edit_single("indochina.json", "インドシナ連邦", 107.6, 13.5, "F")
    edit_single("india.json",     "インド帝国",      76.1, 22.0, "E")

手順 5　データの集約

改めて indochina.json と india.json を southEastAsia.json に追加して主題図用のデータを完成します。最後に無用になった india.json は southEastAsia.json は削除します。

ogr2ogr -update -append southEastAsia.json indochina.json
ogr2ogr -update -append southEastAsia.json india.json
del indochina.json india.json

手順 2~5 を再掲したバッチファイルは以下のようになります。これを適当な名前、例えば make_map.bat にセーブします（全角文字を含むので文字コードは shift-jis）。

: ne_110m_admin_0_countries.zip
:    https://www.naturalearthdata.com/downloads/110m-cultural-vectors/110m-admin-0-countries/

@echo off
setlocal
set OSGEO4W=C:\OSGeo4W
set PATH=%OSGEO4W%\bin;%PATH%
call %OSGEO4W_ROOT%\etc\ini\gdal.bat
call %OSGEO4W_ROOT%\etc\ini\python3.bat

If exist southEastAsia.json del southEastAsia.json
ogr2ogr ^
    -nln southEastAsia ^
    -lco COORDINATE_PRECISION=5 ^
    -sql "SELECT NAME_JA AS NAME, LABEL_X AS X, LABEL_Y AS Y FROM ne_110m_admin_0_countries" ^
    -clipsrc 45 -12 150 45 ^
    clipped.json /vsizip/ne_110m_admin_0_countries.zip

set UNION=SELECT ST_Union(ST_buffer(geometry, 0.0001)) AS geometry
set UNION=%UNION% FROM southEastAsia WHERE NAME IN

ogr2ogr ^
    -lco COORDINATE_PRECISION=5 ^
    -dialect sqlite ^
    -sql "%UNION% ('カンボジア', 'ラオス', 'ベトナム')" ^
    indochina.json southEastAsia.json

ogr2ogr ^
    -lco COORDINATE_PRECISION=5 ^
    -dialect sqlite ^
    -sql "%UNION% ('インド', 'バングラデシュ')" ^
    india.json southEastAsia.json

python edit_layers.py

ogr2ogr -update -append southEastAsia.json indochina.json
ogr2ogr -update -append southEastAsia.json india.json
del indochina.json india.json

手順 6　主題図の作成

make_map

で作成される southEastAsia.json を QGIS に読み込んでお好みの主題図に仕上げます。

ラベル位置の微調整
凡例の追加（ハッチングパターン、カテゴリ名）
印刷レイアウト
その他細々

まとめ

make_map.bat は以前どこかで作ったバッチファイルをかき集めて今回の問題に合わせて微修正したものです。なので、ここで書き下ろしたのは edit_props.py だけです。また当然、一回で主題図のデータができた訳ではなくて何度かの試行錯誤を伴いました。
　つまり、バッチファイルを利用するメリットとしては以下が挙げられます。

何度でも失敗できる
バッチファイルを書き換えることで類似の問題で再利用できる
データ加工の履歴が保存できる

GUI の場合は失敗したら最初からやり直しになるので、一番目の「何度でも失敗できる」は特に大きなメリットです。

CUI で GIS

事 例

作業手順

環境変数の設定

手順 1 基図データのダウンロード

手順 2 東南アジア域のデータの抽出

手順 3 独立前の領地のポリゴンデータの作成

手順 4 属性の変更

手順 5 データの集約

手順 6 主題図の作成

まとめ

事例