SpatiaLiteを使ってSQLだけで面積按分

# はじめに
初のAdventCalendar参加、「PythonからGRASS6の処理を動かす際の下ごしらえ」的なのを書こうかなあと思っていたのですが、[takahi](http://qiita.com/takahi)さんが12/07に投稿された、「[PostGISを使ってSQLだけで面積按分](http://qiita.com/takahi/items/b5823a3fe1c48652e3d2)」の記事が大変参考になったので、同じことをSpatiaLiteでやったらどうなるか？　…という内容にしてみます。
（takahiさん、有益な記事をありがとうございます！）

# SpatiaLiteって？
- データベース管理システム（DBMS）のSQLiteを空間データベースとして扱える拡張機能です。
- 属性データや地物の条件に応じてデータを選択したり、抽出したり、結合したり、まとめたり、空間処理したり、という操作がSQLという言語を使って行えます。
- データは一つのファイル内に収まるので、サイズはともかく見た目はコンパクトです。


# どうしてSQLでやるの？
- メリット①：速度　⇒takahiさんの記事のとおりです。
- メリット②：トランザクション
    - 個人的には、SQLで空間処理を行う大きなメリットとして、トランザクションが利用できることを挙げたい！
    - トランザクションとは、PostgreSQLやSQLiteなどのデータベース管理システムの機能で、一連のSQLによる処理（クエリ）の結果を確定（コミット）するまでは、以前の状態に戻す（ロールバック）ことができる、というものです。
    - 例えば、<br>①たくさんの地物の中から条件に合う地物を抽出して、<br>②抽出した地物を特定の条件でディゾルブ（融合）して、<br>③新しい地物としてデータを出力しつつ、<br>④元のデータから抽出された地物だけを除去する<br>...という処理を行うとすると、処理途中に一時的なデータが出てきます。（例えば extract.shp とか、dissolved.shp とか？）<br>さらに元のデータから地物が一部削除されるので、処理が失敗した！　…という時のために、あらかじめバックアップを作っておく事になると思います。
    - しかーし、トランザクションを活用すれば、処理結果を確認して (ﾉ∀`)ｱﾁｬｰ となってしまっても、**一発ロールバックで全部なかったことにできる**のです。一時ファイルもできません。なんてべんり。（でもごめんなさい、この記事ではトランザクション関係ありません^^;）
- というわけで、あなたもSQLに挑戦してみませんか？

# どうしてSpatiaLiteでやるの？
- ~~ぽすぐれぽすじすへの対抗心~~
- PostgreSQL/PostGISの場合、サーバーへの接続が必要になるが、SQLite/SpatiaLiteはローカルで利用できるため、容易にデータを作成できます。
- Shapeファイルの場合、shxファイルやdbfファイルなどがセットになりファイルの数が多くなりがちですが、SpatiaLiteではひとつのデータベースファイルの中に複数のベクタやテーブルを格納できるので、ファイル管理が簡単です。
- また、Shapeの場合、属性テーブル（dbfファイル）の列名の文字数等に制限がありますが、SpatiaLiteの場合はあまり気にする必要がありません。
- デメリットとしては、SpatiaLiteで作成した地物データをQGISでレイヤをインポートするには、データ作成後に「ジオメトリを登録する」作業が必要で、少しだけ手間だということと、データベースの中身はファイルに接続しなければ確認できないので、データベースファイルからはファイルの内容（データがいくつ入っているのか）が直感的に把握できないことでしょうか。。



# 実行環境
この記事は以下の環境で実行しています。

- PC：Windows8.1（他のプラットフォームでも利用できると思いますが試したことはありません。）
- SpatiaLite：Version 4.3.0a（SpatiaLite-GUI 2.0.0-devel から使用しています。設定が要らないのでやってみるならこれが楽ちん。）

![SpatiaLite-GUI](https://qiita-image-store.s3.amazonaws.com/0/101923/284d6883-7748-73ca-008d-e61e7dd16ad6.png)


# データの準備
SpatiaLiteの場合も、データベースファイルにポリゴンデータが投入されている必要があります。
データベースファイルの作成方法は[準備編](http://qiita.com/yoh_chan/items/82fa113cd2806a384ae4)に載せてあります。

# SQL
PostGISの場合と表記が違う部分は、コメントアウトして下の行に記載してあります。ただし、元の記事から比較しやすいように列名を一部変更しています。（the_geom ⇒ Geometry）

## STEP1:ポリゴンA形状の取り込みと空間検索

```sql:step1
WITH polygon_a AS (
	SELECT
	GeomFromText( -- ST_GeomFromText,ST_GeometryFromTextでも可
		'POLYGON((
			139.72798347473145 35.7437954343205,
			139.71708297729492 35.733763275056724,
			139.72987174987793 35.72588993164637,
			139.73201751708984 35.7215697134585,
			139.73459243774414 35.7215697134585,
			139.74343299865723 35.72839833791958,
			139.73648071289062 35.735226377009106,
			139.72798347473145 35.7437954343205
		))' , 4612
	) AS geo -- ポリゴンAの形状
)
SELECT
	m.KEY_CODE AS meshcode  -- メッシュコード
	, m.T000608001 AS population -- メッシュ内人口
	, Area( m.Geometry, 1 ) AS mesh_area -- メッシュの面積（メートル法） ST_Areaでも可
	-- , ST_AREA( m.Geometry::geography ) AS mesh_area　-- PostGISの場合
FROM
	MESH05339 AS m
	, polygon_a AS p
WHERE
	Intersects( m.Geometry, p.geo ) -- ST_Intersectsでも可
;
```

### SQLの説明：
基本的な説明はtakahiさんのポストで書いてくださっているので、ここではPostGISとSpatiaLiteで違う部分（＋α）について記述したいと思います。

- WITH句はSQLite（SpatiaLite）でも同様に使用することができます。
- WKTからGeometryを生成する関数**`GeomFromText()`**は、PostGISと同名の**`ST_GeometryFromText()`**表記でも利用できます。（SpaitaLiteの大半の関数は、PostGISの関数と同じ名称で利用できるらしい。）
- メッシュ面積の算出
    - 地物の面積を算出する部分<br>**`Area( m.Geometry, 1 )`**<br>は、takahiさんのポストでは<br>**`ST_AREA( m.Geometry::geography )`**<br>でした。
    - 面積を算出する**`Area()`**は、PostGIS表記の**`ST_Area()`**としても動作します。
    - takahiさんの記事ではさらっと書いてありますが、PostGISの方のカッコ内　"**::geography**"　で（緯度経度に基づく）**geometry型**から（メートル法に基づく）**geography型**へキャストした上で面積を（メートル単位で）算出しています。（へええ、こんなことができるんですね！）
    - しかし**SpatiaLiteには "geometry型" と "geography型" の区別はありません**。（ので、同じように書いても処理できません。）その代わり、Area関数にはメートル法で面積を算出するオプションが用意されていますので、ここではこれを利用しました。（二番目の引数（＝1）。ただし、SpatiaLiteにLWGEOMライブラリが実装されている必要があります。）
    - もう一つの方法としては、ジオメトリを適切なCRS（例えばJGD2000-平面直角座標系IX系：EPSG2451）に投影変換して面積を算出する方法もあります。その場合は<br> **`Area( Transform( m.Geometry, 2451 ) )`**<br>
となります。
- 地物の重なりを確認する関数は**`Intersects()`**ですが、やはりPostGIS表記の**`ST_Intersects()`**でも動作します。

### 実行結果
![image](https://qiita-image-store.s3.amazonaws.com/0/101923/36aeaed6-2c4f-46be-0e85-828fb9d972d4.png)


|meshcode|population|mesh_area|
|:--:|--:|:--:|
|53394568|15110|1045744.355781|
|53394569|18057|1045744.355785|
|53394577|17936|1045637.015459|
|53394578|27046|1045637.015455|
|53394579|20954|1045637.015459|
|53394587|27577|1045529.427123|
|53394588|26061|1045529.427119|
|53394589|21776|1045529.427123|
|53394598|23714|1045421.929299|

...SpatiaLiteでもちゃんと重なる地物を選択できています。


## STEP2: 重なり部分の空間形状取得

```sql:step2
WITH polygon_a AS (
    SELECT
    GeomFromText(
        'POLYGON((
            139.72798347473145 35.7437954343205,
            139.71708297729492 35.733763275056724,
            139.72987174987793 35.72588993164637,
            139.73201751708984 35.7215697134585,
            139.73459243774414 35.7215697134585,
            139.74343299865723 35.72839833791958,
            139.73648071289062 35.735226377009106,
            139.72798347473145 35.7437954343205
        ))' , 4612
    ) AS geo -- ポリゴンAの形状
)
SELECT
	m.KEY_CODE AS meshcode  -- メッシュコード
	, Area( m.Geometry, 1 ) AS mesh_area -- メッシュの面積（メートル単位）
	-- , ST_AREA( m.Geometry::geography ) AS mesh_area -- PostGISの場合
	, Area( Intersection( m.Geometry , p.geo ), 1 ) AS and_area -- 重なり部分の面積計算
	-- , ST_AREA( ST_Intersection( m.Geometry, p.geo ) :: geography ) AS and_area -- PostGISの場合
	, AsText( Intersection( m.Geometry , p.geo ) ) AS and_polygon -- 重なり部分の形状をWKTで出力
	-- , ST_AsText( ST_Intersection(m.Geometry , p.geo ) ) AS and_polygon -- PostGISの場合
FROM
	MESH05339 AS m
	, polygon_a AS p
WHERE
	Intersects( m.Geometry, p.geo )
	-- ST_Intersects( m.Geometry, p.geo )
;
```

### SQLの説明：
- ご覧のとおり、上で説明した**`Area()`**内のキャストの部分以外はほとんど同じになっていますね。
- そして例によって、PostGISでの関数表記の"ST_"を加えても同様に動作します。

### 実行結果
![image](https://qiita-image-store.s3.amazonaws.com/0/101923/602d3ac3-1c0b-599a-ee68-b58742b81ca0.png)

ちゃんと地物の形状がWKT（Well Known Text）形式で得られています。


## STEP3: 面積按分

いよいよ面積按分です。
といっても、、ほとんどtakahiさんの記事と同じなので書くことがない。汗

```sql:step3
WITH polygon_a AS (
    SELECT
    ST_GeometryFromText(
        'POLYGON((
            139.72798347473145 35.7437954343205,
            139.71708297729492 35.733763275056724,
            139.72987174987793 35.72588993164637,
            139.73201751708984 35.7215697134585,
            139.73459243774414 35.7215697134585,
            139.74343299865723 35.72839833791958,
            139.73648071289062 35.735226377009106,
            139.72798347473145 35.7437954343205
        ))' , 4612
    ) AS geo --ポリゴンAの形状
)
SELECT
    COUNT( m.'KEY_CODE' ) AS mesh_num , --ポリゴンAが重なったメッシュ数
    SUM(
        m.T000608001 *
        Area( Intersection( m.Geometry, p.geo ), 1 ) / Area( m.Geometry, 1 )
        -- ST_AREA( ST_Intersection( m.Geometry, p.geo ) :: geography ) / ST_AREA( m.Geometry::geography ) -- PostGISの場合
        ) AS population --ポリゴンAの全体の人口
FROM
    MESH05339 AS m
    , polygon_a AS p
WHERE
    Intersects( m.Geometry , p.geo )
;
```

### SQLの説明：
- ご覧のとおり、面積を計算する部分以外h（略
- 行数をカウントする**`COUNT()`**、合計を計算する**`SUM()`**は、SQLiteに標準で実装されている集計関数で、SpatiaLiteでもまったく同様に使えます。

### 実行結果
![image](https://qiita-image-store.s3.amazonaws.com/0/101923/541a3356-fb95-42f8-9741-79ecbc50ff08.png)


|mesh_num|population|
|:--:|:--:|
|9|64265.783909|


- PostGISでの実行結果と同じ結果が得られました！わーい！

# まとめ
- SpatiaLiteでは、PostGISとよく似た文法で各種の空間関数が利用できます。
- でも少し作法が違うところもあります。
- SQL初めて！の人も、ちょっと触ってもらえるとうれしい。（それでトラウマになったらごめんなさい）



# おまけ：処理した地物をQGISで見たいよ！
- 以下のSQLを実行します。

```sql
-- 重なり部分の地物と人口で新しくテーブルを作る
CREATE TABLE intersect_meshgeo AS
WITH polygon_a AS (
    SELECT
    GeomFromText(
        'POLYGON((
            139.72798347473145 35.7437954343205,
            139.71708297729492 35.733763275056724,
            139.72987174987793 35.72588993164637,
            139.73201751708984 35.7215697134585,
            139.73459243774414 35.7215697134585,
            139.74343299865723 35.72839833791958,
            139.73648071289062 35.735226377009106,
            139.72798347473145 35.7437954343205
        ))' , 4612
    ) AS geo -- ポリゴンAの形状
)
SELECT
	m.KEY_CODE AS meshcode  -- メッシュコード
	, m.T000608001 *
        Area( Intersection( m.Geometry, p.geo ), 1 ) / Area( m.Geometry, 1 ) AS population
	, Intersection( m.Geometry , p.geo ) AS Geometry -- 重なり部分の形状をジオメトリとして得る
FROM
	MESH05339 AS m
	, polygon_a AS p
WHERE
	Intersects( m.Geometry, p.geo )
;
-- 作成したテーブルのジオメトリをSpatiaLiteに登録する
SELECT RecoverGeometryColumn(
	'intersect_meshgeo'
	, 'Geometry'
	, 4612
	, 'POLYGON'
	, 'XY'
)
;
```

- 無事に登録されれば、以下のように「1」が返ります。

![image](https://qiita-image-store.s3.amazonaws.com/0/101923/04208fae-8f40-563b-ed25-e1b86357f169.png)

- QGISで<br>レイヤ->レイヤの追加->SpatiaLiteレイヤの追加　と進みます。

![image](https://qiita-image-store.s3.amazonaws.com/0/101923/8430eea1-98c5-c8d5-942e-321402d1fd82.png)

- 「新規」ボタンからデータベースファイルを指定します。
- 指定したら、「接続」でレイヤが表示されます。
- 先ほど作成した intersect_meshgeo を選択して「追加」とすれば、レイヤが追加できます。

![image](https://qiita-image-store.s3.amazonaws.com/0/101923/a2b2c018-d3f2-a31c-05dd-cab5bb19dc51.png)

（あえて結果は見せないでおこうそうしよう。）