トポロジで都道府県を簡略化してみた (これがベストなのかは知らない)

PostGIS

今回のねらい

http://strk.keybit.net/projects/postgis/Paris2011_TopologyWithPostGIS_2_0.pdf にある "Topology with PostGIS 2.0"の7,8枚目を見て、都道府県ポリゴンのsimplifyにトポロジを使えるのではないかと考えた次第。

Typmodさまさま

ということで、既に格納している国土数値情報(行政区域)をゴニョゴニョしようとして、別データベースに移設しようとしました。
こういう場合はPostGIS 2.0では、PostgreSQLのpg_dumpを使って、データをダンプしてリストアしてやればOK。

pg_dump -t (テーブル名) (複写元データベース名) > (ファイル)
psql (複写先データベース名) -f (ファイル)

こんなかんじ。

これまでは、私の場合、pg_dump -a でデータのみ出して、CREATE TABLE と AddGeometryColumn を先に実行したうえでデータを叩き込んでいました*1。
なので、ダンプをそのまんま食ってくれるのにびっくりした。いや、ありがたい。

前提

CREATE TABLE pref (
  gid SERIAL PRIMARY KEY,
  area DOUBLE PRECISION,
  geom GEOMETRY(POLYGON,4612)
);

…があるとします。ポリゴンはST_Dumpでシングルなポリゴンにして叩き込んでます。

面積を出すのにGEOGRAPHYを使う

今回は小さい島等は無視したいので、面積で切ることにしました。
まずは面積を出さなければならないので、次のようにしました。

UPDATE pref SET area=ST_Area(ST_SetSRID(geom,4326)::GEOGRAPHY);

GEOGRAPHY使うとあまり深く考えなくていいです。投影法にはUTM使ってるんだけどね。

トポロジ作成とポリゴン追加

トポロジを作成します。

SELECT topology.CreateTopology('topo1',4612);

50000000m^2=50km^2以上の面積を持つポリゴンをトポロジに追加するには次のようにします。

SELECT TopoGeo_AddPolygon('topo1',geom,0) FROM pref WHERE area >= 50000000;

ここでTopoGeo_AddPolygonの第3引数は、許容差です。ここでは0にしていますが、ポリゴンの境界がズレてそうで、そのズレを許容したい場合に利用します。

この追加は、けっこう時間がかかりますので注意が必要です。

簡略化したラインストリングをtopo2に入れる

SELECT topology.CreateTopology('topo2',4612);

さきほどと同じように、TopoGeo_AddLineStringで ST_SimplifyPreserveTopology で生成したラインストリングを入れていきます。
今回は第2引数を2(単位は度)としました。これは相当削ってしまうので注意です。

SELECT TopoGeo_AddLineString('topo2',ST_SimplifyPreserveTopology(geom, 2)) FROM topo1.edge;

つつがなく終了したら、エッジ(とノード)が入ったことになります。これをポリゴン化します。

SELECT topology.Polygonize('topo2');

TopoGeometryの作成

テーブルとTopoGeometryカラムを作ります。

CREATE TABLE test (
  gid SERIAL PRIMARY KEY,
  face_id INT
);
SELECT AddTopoGeometryColumn('topo2', 'public', 'test', 'topogeom', 'POLYGON');
CREATE INDEX ix_test_face_id ON test (face_id);

testテーブルに挿入します。

INSERT INTO test(topogeom)
SELECT topology.CreateTopoGeom('topo2',3,1,Q.faces)
      FROM (
        SELECT ARRAY[ARRAY[f.face_id,3 ] ] as faces FROM topo2.face AS f
        WHERE NOT mbr IS NULL
      ) AS Q;

ビューを作成したうえでQGISで見る

TopoGeometryはGEOMETRYにキャストできますが、そのままではGEOMETRY型ではないようです。
キャストするためのビューを作成してGEOMETRYとして見てみることにします。

CREATE VIEW v_test AS SELECT gid, topogeom::geometry FROM test;

ポリゴンに対するsimplify

ポリゴン(GEOMETRYのポリゴンであって、トポロジのフェイスではない)ごとにsimplifyしてみましょう。

CREATE TABLE test2 (
  gid SERIAL PRIMARY KEY,
  geom GEOMETRY(POLYGON,4612)
);
CREATE INDEX ix_test2_geom ON test2 USING GiST(geom);

INSERT INTO test2(geom) SELECT ST_SimplifyPreserveTopology(geom, 2) FROM pref WHERE area >= 50000000;

結果を見てみましょう

まず、トポロジを使った場合です。

…ああ、まあ、山口県がすごく小さくなったのは、トポロジのノードが下関とかのあたりに存在しないからです、ていうかsimplifyのやりすぎ。

ポリゴンごとのsimplifyは次のようになりました。

こっちは、ポリゴン同士が重なったり隙間があいたり、となっています。

属性は持てない

トポロジのフェイスとかは属性は持てないので、ポリゴンは出たけど何県かわからない。

まとめ

  • トポロジばんざい

*1:ただ1.4から populate_geometry_columns が出てるので、一回それを実行してやるとOKになるはずなんで、ありがたみが薄いかも知んない