OGR de SQL

GIS でデータベースを扱うことができるソフトといえば、そう、 OGR ですね。

MySQL や PostgreSQL のように効率よく処理するために実データが複数ファイルに分かれていることもあれば、 SQLite のように単ファイルで保存されている形式もあります。広義の定義では構造化されていないテキストファイルもデータベースに含むと考えることもできます。

つまり ESRI Shapefile や GeoJSON であっても十分にデータベースだよね。そんでもって、複数のデータベースを結合させたり SQL などで問い合わせることができるシステムなら、それはもう RDBMS だよね。

ということで RDBMS である OGR ですが、どの程度 SQL が使えるのか確認してみます。

OGR における SQL

ベクタデータを扱う OGR ですが、 ESRI Shapefile をはじめ種々のデータを読み込んだ際に、 SQLite 形式の仮想データとして扱われます。そのため ESRI Shapefile であっても何であっても、データを読み込み後 SQL を実行することができます。

使用することができる SQL の方言は -dialect オプションで指定します。デフォルトで -dialect OGRSQL と -dialect SQLITE をサポートしています。

SQLite の方は SpatiaLite の空間関数 ST_* も利用することができます。

ogrinfo

ogrinfo は一般的にはベクタデータの情報を確認するためのコマンドですが、 SQL を利用してデータ挿入なども可能です。

たとえば ne_10m_airports.shp に対し、地物を追加することもできます。

（レイヤ名を個別に指定できない） ESRI Shapefile 形式なので、レイヤ名（テーブル名）はファイル名と同一である ne_10m_airports 、ジオメトリカラムは GEOMETRY として扱われます。（ GeoPackage など、レイヤ名やジオメトリカラム名を指定することができる形式の場合は ogrinfo を実行し、レイヤ名やジオメトリカラム名を確認する必要があります。）

INSERT INTO ne_10m_airports (GEOMETRY, name) VALUES (ST_Point(135, 35), 'insertion test')

という INSERT 文を実行するには下記のコマンドを実行します。

$ ogrinfo -dialect SQLITE -sql "INSERT INTO ne_10m_airports (GEOMETRY, name) VALUES (ST_Point(135, 35), 'insertion test')" ne_10m_airports.shp

ちなみに INSERT がうまくいったか、地物数の確認をするには

$ ogrinfo -ro -dialect SQLITE -sql "SELECT COUNT(*) FROM ne_10m_airports" ne_10m_airports.shp
INFO: Open of `ne_10m_airports.shp'
      using driver `ESRI Shapefile' successful.

Layer name: SELECT
Geometry: None
Feature Count: 1
Layer SRS WKT:
(unknown)
COUNT(*): Integer (0.0)
OGRFeature(SELECT):0
  COUNT(*) (Integer) = 894

（こんなことをしなくても、普通に ogrinfo をたたけば地物数は確認できますが）

このほか UPDATE や、 -dialect OGRSQL だと ALTER TABLE なども対応しているようです。

ogr2ogr

ベクタデータの形式変換や座標系変換などを行う ogr2ogr ですが、 -sql を使って SELECT で問い合わせた結果をもとに出力ファイルを生成することができます。

条件によるフィルタリング

SELECT * FROM ne_10m_airports WHERE type='major'

を実行し、 WHERE 句を使ったデータの絞り込みを行うには

$  ogr2ogr -f "ESRI Shapefile" -oo ENCODING=UTF-8 -lco ENCODING=UTF-8 -dialect SQLITE -sql "SELECT * FROM ne_10m_airports WHERE type='major'" extract_airports.shp ne_10m_airports.shp

糖衣構文（？）として、 -where オプションで条件のみを記載しても同等の結果が得られます。

空間関数によるジオメトリ生成

SELECT ST_Buffer(GEOMETRY, 0.1) AS GEOMETRY, name FROM ne_10m_airports

を実行し、空港データ（点データ）に対し、バッファを得るには

$  ogr2ogr -f "ESRI Shapefile" -oo ENCODING=UTF-8 -lco ENCODING=UTF-8 -dialect SQLITE -sql "SELECT ST_Buffer(GEOMETRY, 0.1) AS GEOMETRY, name FROM ne_10m_airports" airports_buffer.shp ne_10m_airports.shp

空間関数によるフィルタリング

SELECT * FROM ne_10m_airports
WHERE ST_intersects(GEOMETRY, (
    SELECT nation.GEOMETRY
    FROM 'ne_10m_admin_0_countries.shp'.'ne_10m_admin_0_countries' AS nation
    WHERE nation.NAME='Japan'
))

を実行し、複数データ（空港データと国データ）から日本国内にある空港を抽出してみます。

$ ogr2ogr -f "ESRI Shapefile" -oo ENCODING=UTF-8 -lco ENCODING=UTF-8 -dialect SQLITE -sql "SELECT * FROM ne_10m_airports WHERE ST_intersects(GEOMETRY, (SELECT nation.GEOMETRY FROM 'ne_10m_admin_0_countries.shp'.'ne_10m_admin_0_countries' AS nation WHERE nation.NAME='Japan'))" japan_airports.shp ne_10m_airports.shp

ogr2ogr の構文上は指定できる入力データソースはひとつですが、 SQL 文中でファイル名、レイヤ名（テーブル名）を指定することができます。

空間演算

SELECT ST_Intersection(nation.GEOMETRY, tz.GEOMETRY) AS GEOMETRY,
    nation.NAME, nation.ISO_A3, tz.time_zone
FROM ne_10m_admin_0_countries AS nation,
   'ne_10m_time_zones.shp'.'ne_10m_time_zones' AS tz
WHERE nation.ISO_A3='USA' AND ST_Intersects(nation.GEOMETRY, tz.GEOMETRY)

を実行し、複数データ（タイムゾーンデータと国データ）からアメリカのタイムゾーン分割を行ってみます。

$ ogr2ogr -f "ESRI Shapefile" -oo ENCODING=UTF-8 -lco ENCODING=UTF-8 -dialect SQLITE -sql "SELECT ST_Intersection(nation.GEOMETRY, tz.GEOMETRY) AS GEOMETRY, nation.NAME, nation.ISO_A3, tz.time_zone FROM ne_10m_admin_0_countries AS nation, 'ne_10m_time_zones.shp'.'ne_10m_time_zones' AS tz WHERE nation.ISO_A3='USA' AND ST_Intersects(nation.GEOMETRY, tz.GEOMETRY)" USA_tz_split.shp ne_10m_admin_0_countries.shp

空間集計

SELECT nation.GEOMETRY, nation.NAME,
    COUNT(ap.NAME) AS n_airports,
    MIN(ap.scalerank) AS min_scale
FROM ne_10m_admin_0_countries AS nation,
    'ne_10m_airports.shp'.'ne_10m_airports' AS ap
WHERE ST_Intersects(nation.GEOMETRY, ap.GEOMETRY)
GROUP BY nation.GEOMETRY

を実行し、複数データ（空港データと国データ）から各国ごとに空港がいくつあるか集計し、もっとも小さな scalerank の値を得るには。

$ ogr2ogr -f "ESRI Shapefile" -oo ENCODING=UTF-8 -lco ENCODING=UTF-8 -dialect SQLITE -sql "SELECT nation.GEOMETRY, nation.NAME, COUNT(ap.NAME) AS n_airports, MIN(ap.scalerank) AS min_scale FROM ne_10m_admin_0_countries AS nation, 'ne_10m_airports.shp'.'ne_10m_airports' AS ap WHERE ST_Intersects(nation.GEOMETRY, ap.GEOMETRY) GROUP BY nation.GEOMETRY" aggregate_airports.shp ne_10m_admin_0_countries.shp

空間インデックスが張られていないこともあるのかもしれませんが、処理に結構時間がかかりました。

まとめ

普段使いませんが今回は敢えて ESRI Shapefile のデータを利用し、任意の GIS データに対し SQL 文で処理を行ってみました。

パフォーマンス的にはあまりよくないとは思いますが、手元に Shapefile 形式などのデータがあり、わざわざ PostGIS や SpatiaLite に変換するのもちょっと面倒。だけど条件で抽出したい。とかそういったケースで便利かもしれません。