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

はじめに

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

SpatiaLiteって？

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

どうしてSQLでやるの？

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

SQL

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

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

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

実行結果

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: 重なり部分の空間形状取得

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_"を加えても同様に動作します。

実行結果

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

STEP3: 面積按分

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

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でもまったく同様に使えます。

実行結果

mesh_num	population
9	64265.783909

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

まとめ

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

おまけ：処理した地物をQGISで見たいよ！

• 以下の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」が返ります。

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

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

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