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` を実行し、レイヤ名やジオメトリカラム名を確認する必要があります。）

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

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

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

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

```bash
$ 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` で問い合わせた結果をもとに出力ファイルを生成することができます。

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

```sql
SELECT * FROM ne_10m_airports WHERE type='major'
```

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

```bash
$  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` オプションで条件のみを記載しても同等の結果が得られます。

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

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

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

```bash
$  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
```

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

```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'
))
```

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

```bash
$ 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 文中でファイル名、レイヤ名（テーブル名）を指定することができます。

### 空間演算

```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)
```

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

```bash
$ 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
```

### 空間集計

```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
```

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

```bash
$ 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 に変換するのもちょっと面倒。だけど条件で抽出したい。とかそういったケースで便利かもしれません。