はじめに

今回は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人である事が分かりました。

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