前回用いた e-Stat 統計地理情報システムのデータは 1 次メッシュ毎にアーカイブされています。今回は、前回の事例の対象流域（clipsrc）をカヴァーする 1 次メッシュをどうやって調べたか、について記します。
Python API
ogr2ogr は地理空間データを抽象化するクラスライブラリ GDAL を用いて C++ で作成されたプログラムです。同様に Python でも GDAL を用いてプログラムを自作することができます。このための環境も QGIS に含まれています。さらに言えば、ogrmerge.py などいくつかのプログラムが Python で実装されています。
GDAL は多機能なので全てを把握するのは困難ですが、主な使用法についてクックブック（英文 org）が提供されています。ニーズが低いのか 2.7 版の Python での記事になっていますが、それでも十分に参考になります。
meshCode.py
冒頭の用途のために Python API を用いて記述した meshCode.py は以下の通りです。
入力ファイル（file_input）内の特定のレイヤーの地物（feature）あるいはレイヤー内の地物群を包含する 1 次メッシュのコードを列挙します。file_input のフォーマットは Shapefile、SQLite、GeoJSON あるいは GeoPackge の何れかを想定しています。
33 行目以降の行で地物毎の北西端・南東端の座標（envelope）から 1 次メッシュのコードを割り出しています。
import os
import sys
from osgeo import ogr
if __name__ == '__main__':
file_input = sys.argv[1]
if not os.path.exists(file_input):
sys.exit(f"\n {file_input} not found.")
for driver_name in ("ESRI Shapefile", "SQLite", "GeoJSON", "GPKG"):
driver = ogr.GetDriverByName(driver_name)
layers = driver.Open(file_input, 0)
if not layers is None: break
if layers is None:
sys.exit(f"\n can not open {file_input}.")
if layers.GetLayerCount() == 1:
layer = layers.GetLayer(0)
layer_name = layer.GetName()
else:
if len(sys.argv) == 1:
sys.exit("\n please input layer name.")
layer_names = []
for layer in layers:
layer_names.append(layer.GetName())
layer_name = sys.argv[1]
if layer_name in layer_names:
layer = layers.GetLayer(layer_name)
else:
sys.exit(f"\n {layer_name} no such layer.")
feature_count = layer.GetFeatureCount()
code_set = []
for idx in range(feature_count):
feature = layer.GetNextFeature()
geometry = feature.GetGeometryRef()
envelope = geometry.GetEnvelope()
codes = []
for lon in [envelope[0], envelope[1]]:
for lat in [envelope[2], envelope[3]]:
code = f"{int(lat * 1.5)}{int(lon - 100)}" # 1st order
if not code in codes:
codes.append(code)
if not code in code_set:
code_set.append(code)
print(f"fid:{idx + 1:3d} {codes}")
if feature_count > 1:
print(f"\n union {code_set}")meshCode.bat
3 行目で宣言している ogr モジュールをインポートするため、QGIS 同梱の Python を使用するので、バッチファイル meshCode.bat は次のようになります。meshCode.py はレイヤーが 1 つしかない場合はレイヤー名を省略できる仕様にしています。
@echo off
setlocal
set OSGEO4W_ROOT=C:\OSGEO4W
set PATH=%OSGEO4W_ROOT%\bin;%PATH%
call %OSGEO4W_ROOT%\etc\ini\gdal.bat
call %OSGEO4W_ROOT%\etc\ini\python3.bat
if not "%1" == "" (
python meshCode.py %*
) else (
echo.
echo Usage : %~n0 src_dataset_name [layer_name]
)前回の事例では
meshCode someBasin.gpkgとしてレイヤー内の全ての小流域を包含する 1 次メッシュのコード 5239、5339 を得ました。
まとめ
先に「33 行目以降の行で～」と記しましたが、逆にこれより前のスクリプトは問題によらない定型の処理です。したがって、実質は 20 行足らずで記述されていることになります。それほど敷居は高くないのではないでしょうか。
なお、Python の心得のある人はとかく自作に走り勝ちですが、これはお勧めできません。出来合いのプログラムを最大限活用して、足らざる機能を自作で補足するというスタンスがタイパ的に優れていると思います。ま、飽くことなきマウス操作のタイパは最低レベルですが・・・