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

# はじめに
 今回はSQLでの面積按分を取り上げます。GISエンジニアの方々にSQLの面白さに目覚めていただきたい一心で投稿しております。
## どうしてSQLでやるの？
・ メリット：速度
データベースに投入されているデータを、RやQGISなどのアプリケーションソフトで処理するには、必ず通信が発生します。データベースの内部で実行できれば、通信のオーバヘッドが無くなり高速化が期待できます。
## 面積按分
 面積按分についての詳しい説明は[他のサイト](http://www.pasco.co.jp/recommend/word/word050/ "PASCO")を御覧ください。
大雑把な説明としては、「切り抜い形状と、もとの形状の面積の比率を求め、統計量もその比率で切り出す計算」といった感じです。
## 実行環境
この記事は以下の環境で実行しています。基本的にはPostGISが動けば実行可能であると思います。

+ データベース：PostgreSQL 9.4 (PostGISが動作すればOK）
+ データ：男女別人口及び世帯総数（e-Stat: 平成22年国勢調査 世界測地系１ｋｍメッシュ）
+ PC：Mac OS X（PostgreSQLとQGISが動作する環境ならOKです）

## データの準備
PostgreSQLにポリゴンデータが投入されている必要があります。
説明が長くなりそうなので、[準備編](http://qiita.com/takahi/items/9d44966b10d6257d8b2d)でご説明しております。

+ データ：１ｋｍ四方（3次メッシュ）の人口統計
 <img width="764" alt="スクリーンショット 2015-12-06 22.21.55.png" src="https://qiita-image-store.s3.amazonaws.com/0/53918/ca03a23f-c3ec-5498-c8b1-556d940566cb.png">


+ 人口を計算したいエリア
 春日通り→不忍通り→白山通り→明治通りで囲まれた範囲
 今後はこのエリアをポリゴンAと呼びます。
![スクリーンショット 2015-12-06 22.00.16.png](https://qiita-image-store.s3.amazonaws.com/0/53918/f75eca8b-f78b-1941-a9b4-fcde8c7ddd5a.png)
© OpenStreetMap contributors
+ ポリゴンAのWKT

```sql:wkt
SELECT '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
))' AS wkt 
```

+ 人口データをポリゴンAを重ねた図 
 暗算すると5万人以上になりそうです。
 ポリゴンAは９個のメッシュに重なっていることも分かります。<img width="979" alt="スクリーンショット 2015-12-06 22.34.22.png" src="https://qiita-image-store.s3.amazonaws.com/0/53918/e8fe4587-8851-5144-b30e-25e76603fe60.png">

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

ポリゴンAをgeometry型に変換し、人口データと重なりあう空間を検索

```sql:step1
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 
   m."KEY_CODE" AS meshcode  -- メッシュコード
 , m."tblT000608" AS population -- メッシュ内人口
 , ST_AREA( m.the_geom::geography ) AS mesh_area -- メッシュの面積（メートル法）
FROM mesh5339 AS m , polygon_a AS p
WHERE 
 ST_Intersects( m.the_geom , p.geo )
```
+ SQLの説明
    + WITH句：下のSELECTを実行する前のWITH句でデータを作成します。SELECT句をネストしても構いませんが、WITH句で明示することで分かりやすくなります。今回はST_GeometryFromText関数を用いて、ポリゴンAの形状をWKTからgeometry型に変換しています。WITH句作成したテーブル名は「polygon_a」になります。
    + ST_AREA：面積を計算する関数です。the_geomという１ｋｍメッシュの形状をgeography型にキャストすることにより、メートル単位の面積計算結果を得ることが出来ます。
    + ST_Intersects:形状が重なりあうかを判定する関数です。ポリゴンAと重なりあうメッシュを検索しています。
    
+ 実行結果
 9行（９個）のメッシュがヒットし、QGISで表示したものと同様の結果が得られている事が分かります。

| meshcode| population |mesh_area|
|:-----------:|------------:|:------------:|
|53394568|	15110|	1045744.35578108|
|53394569|	18057|	1045744.35578513|
|53394577|	17936|	1045637.01545906|
|53394578|	27046|	1045637.01545453|
|53394579|	20954|	1045637.01545906|
|53394587|	27577|	1045529.42712307|
|53394588|	26061|	1045529.42711902|
|53394589|	21776|	1045529.42712307|
|53394598|	23714|	1045421.92929935|

##STEP2: 重なり部分の空間形状取得
PostGISのST_Intersection関数を用いることで、２つの形状が重なりあう部分だけを抽出可能です。


```sql:step2
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 
   m."KEY_CODE" AS meshcode  -- メッシュコード
 , m."tblT000608" AS population -- メッシュ内人口
 , ST_AREA( m.the_geom::geography ) AS mesh_area -- メッシュの面積（メートル単位）
 , ST_AREA( ST_Intersection(m.the_geom , p.geo ) :: geography ) AS and_area -- 重なり部分の面積計算
 , ST_AsText( ST_Intersection(m.the_geom , p.geo ) ) AS and_polygon -- 重なり部分の形状をWKTで出力
FROM mesh5339 AS m , polygon_a AS p
WHERE 
 ST_Intersects( m.the_geom , p.geo )
```
+ SQLの解説
 + ST_Intersection：重なりあうエリアを返却する
 + ST_AREA(ST_Intersection(a,b)::geography): 重なりあうエリアの面積をメートル単位で計算する
 + ST_AsText：geometry型をWKT形式に変換する関数。
  QGISのディミテッドテキストレイヤの作成で表示が可能です
+ 結果をQGISで表示
 9個のポリゴンを色を変えて表示しています。
<img width="654" alt="スクリーンショット 2015-12-06 23.44.50.png" src="https://qiita-image-store.s3.amazonaws.com/0/53918/f3842ebc-f12a-eb56-0b48-b88f77af51cf.png">

これで面積按分に必要な、重なりあう空間の面積を得る事が出来ました。



##STEP３: 面積按分
STEP２で面積按分に必要なパラメータを得ることが出来ました。
面積按分による人口計算を以下の式で行います。

人口 × （ 重なり部分の面積 / メッシュの面積 ）

```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."tblT000608" * (
	 ST_AREA( ST_Intersection(m.the_geom , p.geo ) :: geography ) / ST_AREA( m.the_geom::geography ) 
	)
 ) AS population --ポリゴンAの全体の人口
FROM mesh5339 AS m , polygon_a AS p
WHERE 
 ST_Intersects( m.the_geom , p.geo )
```
+ SQLの解説
    + COUNT:行数をカウントする関数、今回は「行数＝メッシュ数」
    + SUM:合計値を計算する関数
+ 実行結果

| mesh_num | population |
|:-----------:|:------------:|
| 9|	64265.7839089188|
少数を切り上げすると、エリアAの人口は約64,266人である事が分かりました。