Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
20
Help us understand the problem. What is going on with this article?
@takahi

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

はじめに

今回は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のディミテッドテキストレイヤの作成で表示が可能です
  • 結果をQGISで表示 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人である事が分かりました。

20
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
takahi
このサイトの掲載内容は私個人の見解であり、必ずしも私が所属する会社、組織、団体の立場、戦略、意見を代表するものではありません。

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
20
Help us understand the problem. What is going on with this article?