はじめに
今回はSQLでの面積按分を取り上げます。GISエンジニアの方々にSQLの面白さに目覚めていただきたい一心で投稿しております。
どうしてSQLでやるの？
・ メリット：速度
データベースに投入されているデータを、RやQGISなどのアプリケーションソフトで処理するには、必ず通信が発生します。データベースの内部で実行できれば、通信のオーバヘッドが無くなり高速化が期待できます。
面積按分
面積按分についての詳しい説明は他のサイトを御覧ください。
大雑把な説明としては、「切り抜い形状と、もとの形状の面積の比率を求め、統計量もその比率で切り出す計算」といった感じです。
実行環境
この記事は以下の環境で実行しています。基本的にはPostGISが動けば実行可能であると思います。
- データベース：PostgreSQL 9.4 (PostGISが動作すればOK）
- データ：男女別人口及び世帯総数（e-Stat: 平成22年国勢調査 世界測地系１ｋｍメッシュ）
- PC：Mac OS X（PostgreSQLとQGISが動作する環境ならOKです）
データの準備
PostgreSQLにポリゴンデータが投入されている必要があります。
説明が長くなりそうなので、準備編でご説明しております。
- データ：１ｋｍ四方（3次メッシュ）の人口統計
- 人口を計算したいエリア
春日通り→不忍通り→白山通り→明治通りで囲まれた範囲
今後はこのエリアをポリゴンAと呼びます。
© OpenStreetMap contributors - ポリゴンAのWKT
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は９個のメッシュに重なっていることも分かります。
#SQL
##STEP1: ポリゴンA形状の取り込みと空間検索
ポリゴンAをgeometry型に変換し、人口データと重なりあう空間を検索
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関数を用いることで、２つの形状が重なりあう部分だけを抽出可能です。
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個のポリゴンを色を変えて表示しています。
これで面積按分に必要な、重なりあう空間の面積を得る事が出来ました。
##STEP３: 面積按分
STEP２で面積按分に必要なパラメータを得ることが出来ました。
面積按分による人口計算を以下の式で行います。
人口 × （ 重なり部分の面積 / メッシュの面積 ）
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人である事が分かりました。 |