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

  • 13
    いいね
  • 0
    コメント
この記事は最終更新日から1年以上が経過しています。

はじめに

今回はSQLでの面積按分を取り上げます。GISエンジニアの方々にSQLの面白さに目覚めていただきたい一心で投稿しております。

どうしてSQLでやるの?

・ メリット:速度
データベースに投入されているデータを、RやQGISなどのアプリケーションソフトで処理するには、必ず通信が発生します。データベースの内部で実行できれば、通信のオーバヘッドが無くなり高速化が期待できます。

面積按分

面積按分についての詳しい説明は他のサイトを御覧ください。
大雑把な説明としては、「切り抜い形状と、もとの形状の面積の比率を求め、統計量もその比率で切り出す計算」といった感じです。

実行環境

この記事は以下の環境で実行しています。基本的にはPostGISが動けば実行可能であると思います。

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

データの準備

PostgreSQLにポリゴンデータが投入されている必要があります。
説明が長くなりそうなので、準備編でご説明しております。

  • データ:1km四方(3次メッシュ)の人口統計
    スクリーンショット 2015-12-06 22.21.55.png

  • 人口を計算したいエリア
    春日通り→不忍通り→白山通り→明治通りで囲まれた範囲
    今後はこのエリアをポリゴンAと呼びます。
    スクリーンショット 2015-12-06 22.00.16.png
    © 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は9個のメッシュに重なっていることも分かります。スクリーンショット 2015-12-06 22.34.22.png

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という1kmメッシュの形状をgeography型にキャストすることにより、メートル単位の面積計算結果を得ることが出来ます。
    • ST_Intersects:形状が重なりあうかを判定する関数です。ポリゴンAと重なりあうメッシュを検索しています。
  • 実行結果
    9行(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関数を用いることで、2つの形状が重なりあう部分だけを抽出可能です。

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のディミテッドテキストレイヤの作成で表示が可能です
  • 結果をQJISで表示 9個のポリゴンを色を変えて表示しています。 スクリーンショット 2015-12-06 23.44.50.png

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

STEP3: 面積按分

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

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

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

この投稿は FOSS4G Advent Calendar 20157日目の記事です。