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

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

はじめに

初のAdventCalendar参加、「PythonからGRASS6の処理を動かす際の下ごしらえ」的なのを書こうかなあと思っていたのですが、takahiさんが12/07に投稿された、「PostGISを使ってSQLだけで面積按分」の記事が大変参考になったので、同じことをSpatiaLiteでやったらどうなるか? …という内容にしてみます。
(takahiさん、有益な記事をありがとうございます!)

SpatiaLiteって?

  • データベース管理システム(DBMS)のSQLiteを空間データベースとして扱える拡張機能です。
  • 属性データや地物の条件に応じてデータを選択したり、抽出したり、結合したり、まとめたり、空間処理したり、という操作がSQLという言語を使って行えます。
  • データは一つのファイル内に収まるので、サイズはともかく見た目はコンパクトです。

どうしてSQLでやるの?

  • メリット①:速度 ⇒takahiさんの記事のとおりです。
  • メリット②:トランザクション
    • 個人的には、SQLで空間処理を行う大きなメリットとして、トランザクションが利用できることを挙げたい!
    • トランザクションとは、PostgreSQLやSQLiteなどのデータベース管理システムの機能で、一連のSQLによる処理(クエリ)の結果を確定(コミット)するまでは、以前の状態に戻す(ロールバック)ことができる、というものです。
    • 例えば、
      ①たくさんの地物の中から条件に合う地物を抽出して、
      ②抽出した地物を特定の条件でディゾルブ(融合)して、
      ③新しい地物としてデータを出力しつつ、
      ④元のデータから抽出された地物だけを除去する
      ...という処理を行うとすると、処理途中に一時的なデータが出てきます。(例えば extract.shp とか、dissolved.shp とか?)
      さらに元のデータから地物が一部削除されるので、処理が失敗した! …という時のために、あらかじめバックアップを作っておく事になると思います。
    • しかーし、トランザクションを活用すれば、処理結果を確認して (ノ∀`)アチャー となってしまっても、一発ロールバックで全部なかったことにできるのです。一時ファイルもできません。なんてべんり。(でもごめんなさい、この記事ではトランザクション関係ありません^^;)
  • というわけで、あなたもSQLに挑戦してみませんか?

どうしてSpatiaLiteでやるの?

  • ぽすぐれぽすじすへの対抗心
  • PostgreSQL/PostGISの場合、サーバーへの接続が必要になるが、SQLite/SpatiaLiteはローカルで利用できるため、容易にデータを作成できます。
  • Shapeファイルの場合、shxファイルやdbfファイルなどがセットになりファイルの数が多くなりがちですが、SpatiaLiteではひとつのデータベースファイルの中に複数のベクタやテーブルを格納できるので、ファイル管理が簡単です。
  • また、Shapeの場合、属性テーブル(dbfファイル)の列名の文字数等に制限がありますが、SpatiaLiteの場合はあまり気にする必要がありません。
  • デメリットとしては、SpatiaLiteで作成した地物データをQGISでレイヤをインポートするには、データ作成後に「ジオメトリを登録する」作業が必要で、少しだけ手間だということと、データベースの中身はファイルに接続しなければ確認できないので、データベースファイルからはファイルの内容(データがいくつ入っているのか)が直感的に把握できないことでしょうか。。

実行環境

この記事は以下の環境で実行しています。

  • PC:Windows8.1(他のプラットフォームでも利用できると思いますが試したことはありません。)
  • SpatiaLite:Version 4.3.0a(SpatiaLite-GUI 2.0.0-devel から使用しています。設定が要らないのでやってみるならこれが楽ちん。)

SpatiaLite-GUI

データの準備

SpatiaLiteの場合も、データベースファイルにポリゴンデータが投入されている必要があります。
データベースファイルの作成方法は準備編に載せてあります。

SQL

PostGISの場合と表記が違う部分は、コメントアウトして下の行に記載してあります。ただし、元の記事から比較しやすいように列名を一部変更しています。(the_geom ⇒ Geometry)

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

step1
WITH polygon_a AS (
    SELECT
    GeomFromText( -- ST_GeomFromText,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.T000608001 AS population -- メッシュ内人口
    , Area( m.Geometry, 1 ) AS mesh_area -- メッシュの面積(メートル法) ST_Areaでも可
    -- , ST_AREA( m.Geometry::geography ) AS mesh_area -- PostGISの場合
FROM
    MESH05339 AS m
    , polygon_a AS p
WHERE
    Intersects( m.Geometry, p.geo ) -- ST_Intersectsでも可
;

SQLの説明:

基本的な説明はtakahiさんのポストで書いてくださっているので、ここではPostGISとSpatiaLiteで違う部分(+α)について記述したいと思います。

  • WITH句はSQLite(SpatiaLite)でも同様に使用することができます。
  • WKTからGeometryを生成する関数GeomFromText()は、PostGISと同名のST_GeometryFromText()表記でも利用できます。(SpaitaLiteの大半の関数は、PostGISの関数と同じ名称で利用できるらしい。)
  • メッシュ面積の算出
    • 地物の面積を算出する部分
      Area( m.Geometry, 1 )
      は、takahiさんのポストでは
      ST_AREA( m.Geometry::geography )
      でした。
    • 面積を算出するArea()は、PostGIS表記のST_Area()としても動作します。
    • takahiさんの記事ではさらっと書いてありますが、PostGISの方のカッコ内 "::geography" で(緯度経度に基づく)geometry型から(メートル法に基づく)geography型へキャストした上で面積を(メートル単位で)算出しています。(へええ、こんなことができるんですね!)
    • しかしSpatiaLiteには "geometry型" と "geography型" の区別はありません。(ので、同じように書いても処理できません。)その代わり、Area関数にはメートル法で面積を算出するオプションが用意されていますので、ここではこれを利用しました。(二番目の引数(=1)。ただし、SpatiaLiteにLWGEOMライブラリが実装されている必要があります。)
    • もう一つの方法としては、ジオメトリを適切なCRS(例えばJGD2000-平面直角座標系IX系:EPSG2451)に投影変換して面積を算出する方法もあります。その場合は
      Area( Transform( m.Geometry, 2451 ) )
      となります。
  • 地物の重なりを確認する関数はIntersects()ですが、やはりPostGIS表記のST_Intersects()でも動作します。

実行結果

image

meshcode population mesh_area
53394568 15110 1045744.355781
53394569 18057 1045744.355785
53394577 17936 1045637.015459
53394578 27046 1045637.015455
53394579 20954 1045637.015459
53394587 27577 1045529.427123
53394588 26061 1045529.427119
53394589 21776 1045529.427123
53394598 23714 1045421.929299

...SpatiaLiteでもちゃんと重なる地物を選択できています。

STEP2: 重なり部分の空間形状取得

step2
WITH polygon_a AS (
    SELECT
    GeomFromText(
        '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  -- メッシュコード
    , Area( m.Geometry, 1 ) AS mesh_area -- メッシュの面積(メートル単位)
    -- , ST_AREA( m.Geometry::geography ) AS mesh_area -- PostGISの場合
    , Area( Intersection( m.Geometry , p.geo ), 1 ) AS and_area -- 重なり部分の面積計算
    -- , ST_AREA( ST_Intersection( m.Geometry, p.geo ) :: geography ) AS and_area -- PostGISの場合
    , AsText( Intersection( m.Geometry , p.geo ) ) AS and_polygon -- 重なり部分の形状をWKTで出力
    -- , ST_AsText( ST_Intersection(m.Geometry , p.geo ) ) AS and_polygon -- PostGISの場合
FROM
    MESH05339 AS m
    , polygon_a AS p
WHERE
    Intersects( m.Geometry, p.geo )
    -- ST_Intersects( m.Geometry, p.geo )
;

SQLの説明:

  • ご覧のとおり、上で説明したArea()内のキャストの部分以外はほとんど同じになっていますね。
  • そして例によって、PostGISでの関数表記の"ST_"を加えても同様に動作します。

実行結果

image

ちゃんと地物の形状がWKT(Well Known Text)形式で得られています。

STEP3: 面積按分

いよいよ面積按分です。
といっても、、ほとんどtakahiさんの記事と同じなので書くことがない。汗

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.T000608001 *
        Area( Intersection( m.Geometry, p.geo ), 1 ) / Area( m.Geometry, 1 )
        -- ST_AREA( ST_Intersection( m.Geometry, p.geo ) :: geography ) / ST_AREA( m.Geometry::geography ) -- PostGISの場合
        ) AS population --ポリゴンAの全体の人口
FROM
    MESH05339 AS m
    , polygon_a AS p
WHERE
    Intersects( m.Geometry , p.geo )
;

SQLの説明:

  • ご覧のとおり、面積を計算する部分以外h(略
  • 行数をカウントするCOUNT()、合計を計算するSUM()は、SQLiteに標準で実装されている集計関数で、SpatiaLiteでもまったく同様に使えます。

実行結果

image

mesh_num population
9 64265.783909
  • PostGISでの実行結果と同じ結果が得られました!わーい!

まとめ

  • SpatiaLiteでは、PostGISとよく似た文法で各種の空間関数が利用できます。
  • でも少し作法が違うところもあります。
  • SQL初めて!の人も、ちょっと触ってもらえるとうれしい。(それでトラウマになったらごめんなさい)

おまけ:処理した地物をQGISで見たいよ!

  • 以下のSQLを実行します。
-- 重なり部分の地物と人口で新しくテーブルを作る
CREATE TABLE intersect_meshgeo AS
WITH polygon_a AS (
    SELECT
    GeomFromText(
        '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.T000608001 *
        Area( Intersection( m.Geometry, p.geo ), 1 ) / Area( m.Geometry, 1 ) AS population
    , Intersection( m.Geometry , p.geo ) AS Geometry -- 重なり部分の形状をジオメトリとして得る
FROM
    MESH05339 AS m
    , polygon_a AS p
WHERE
    Intersects( m.Geometry, p.geo )
;
-- 作成したテーブルのジオメトリをSpatiaLiteに登録する
SELECT RecoverGeometryColumn(
    'intersect_meshgeo'
    , 'Geometry'
    , 4612
    , 'POLYGON'
    , 'XY'
)
;
  • 無事に登録されれば、以下のように「1」が返ります。

image

  • QGISで
    レイヤ->レイヤの追加->SpatiaLiteレイヤの追加 と進みます。

image

  • 「新規」ボタンからデータベースファイルを指定します。
  • 指定したら、「接続」でレイヤが表示されます。
  • 先ほど作成した intersect_meshgeo を選択して「追加」とすれば、レイヤが追加できます。

image

(あえて結果は見せないでおこうそうしよう。)

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