Database JUNKY

MySQL,MariaDBを中心としたブログです

PostgreSQL 位置参照データベースを作ってみる(GISデータ)

そもそも地図はよくわからないので、もう地図のことは書くまいと心に決めていたのですが、意外と見ていただいている方は多いので、やっぱり地図データベースのことを書きます。前回は、MySQLをベースに書きましたが、今回は、PostgreSQLで、PostGISを利用した位置参照データベースの構築、利用方法について書いてみたいと思います。PostGISとはPostgreSQL Geographic Information Systemの略で、まあなんでしょ?地図情報、位置情報などのライブラリを含めたツールの総称です。このページを見ているかたは、おおよそ理解しているかと思いますので細かい話は割愛しますね。逆におおよそ理解している方はこの内容を読むかどうかは別の話ですが・・・(汗

このページの目的としては、

・なんか地図のデータを使って ・なんかそれっぽい結果を出して ・それをどう使うかはその人次第

という方向で進めます、PostgreSQLPostGISの細かい機能、細かい意味についてはあまり触れません(というか全容は自分もよくわかっていません)、とにかくそれっぽく使えるまでの手順を書いていきたいと思います。でも初めて使ってみる方、必見!!!

[pgis] PostgreSQLといいましたが・・・ 実際この検証のページは、厳密にいいますと、PostgreSQLを使っているのではなく、Postgres Plusを利用しています。Postgres Plusを採用した大きな理由は、「PostGISや、PostgreSQLのインストールに時間をかけたくなかった」からです。Postgres Plusは、必要なサードパーティーモジュールがインプリメントされており、インストールも楽チンです。ダウンロードはSIOSテクノロジー株式会社様のPostgres Plusユーザーサイトからダウンロードできます(ユーザー登録は必要)を利用しています。PostGISに関する資料はありませんが、その他インストール方法等、わかりやすい資料が充実していますので、これから導入を検討する方はこれを機会にユーザー登録をしてみてはいかがでしょうか? 今回は以下のような環境で検証しました。

[shell] OS CentOS 5.5 x86-64 DISK 40GB Memory 4GB PostgreSQL バージョン PostgreSQL 8.4.1 PostGIS バージョン POSTGIS="1.4.0" GEOS="3.1.1-CAPI-1.6.0" [/shell]

ジオメトリデータベースの作成 PostGISを活用するためには、まず、ジオメトリデータベースを作成する必要があります。また、前述の通り、すでにPostgres PLUS がインストールされていることが前提で構築手順を説明します。

・pgisデータベースの構築 ・pgisという名前でデータベースを構築します(データベースの名前はなんでもかまいません) 。

[shell]

createdb pgis -E UTF-8 -U postgres

[/shell]

・ジオメトリ関連のテーブル設定(テンプレート情報の取り込み) PostGISを利用するために、当該データベースに、テンプレートテーブルを取り込む必要があります。下記手順は、Postgre PLUSのデフォルトインストールを想定して記載されておりますので、環境によりパスは適宜置き換えてください。 [shell]

psql pgis -U postgres -f /opt/PostgresPlus/8.4SS/share/postgresql/contrib/postgis.sql

psql pgis -U postgres -f /opt/PostgresPlus/8.4SS/share/postgresql/contrib/spatial_ref_sys.sql

[/shell]

・データベースが作成されたか確認する ちょっと手順が前後しましたが、pgisデータベースが生成されたかは?以下のコマンドで確認できます。

[sql]

psql -U postgres -l

List of databases Name       |  Owner   | Encoding | Collation  |   Ctype    |   Access privileges ------------------+----------+----------+------------+------------+----------------------- pgis             | postgres | UTF8     | ja_JP.utf8 | ja_JP.utf8 | postgres         | postgres | UTF8     | ja_JP.utf8 | ja_JP.utf8 | template0        | postgres | UTF8     | ja_JP.utf8 | ja_JP.utf8 | =c/postgres : postgres=CTc/postgres template1        | postgres | UTF8     | ja_JP.utf8 | ja_JP.utf8 | =c/postgres : postgres=CTc/postgres template_postgis | postgres | UTF8     | ja_JP.utf8 | ja_JP.utf8 | [/sql]

データベースの構築は以上で完了です。

・街区情報の設定 前回MySQLで説明したものとほぼ同様ですが、今回の手順も国土交通省の街区情報をベースに検討します。

・ 街区情報テーブルの作成 以下のような、DDLを作りました。方言がありますので、若干MySQLと定義が違いますが、大きな違いは無いでしょう。

[sql] CREATE TABLE geos ( id SERIAL PRIMARY KEY, prefecture varchar(64) DEFAULT NULL, city varchar(128) DEFAULT NULL, area varchar(128) DEFAULT NULL, jiban char(64) DEFAULT NULL, zaban int DEFAULT NULL, x decimal(7,1) DEFAULT NULL, y decimal(7,1) DEFAULT NULL, lat decimal(8,5) DEFAULT NULL, lng decimal(8,5) DEFAULT NULL, jukyo int DEFAULT NULL, daihyo int DEFAULT NULL, history16 int DEFAULT NULL, history17 int DEFAULT NULL, updated timestamp NOT NULL DEFAULT current_timestamp, created timestamp NOT NULL DEFAULT current_timestamp ) ; [/sql]

ジオメトリ型カラムを追加 ここだけでは、MySQLと違います。下記SQLは、geomという名前でgeosテーブルに二次元のジオメトリカラムを生成しています。この、「後で追加」する流儀が理解できないのですが、他のサイトを色々と参考にすると同じような手順でやっているので、その流儀に従ってみます。

[sql]

SELECT AddGeometryColumn('public', 'geos', 'geom', 4326, 'POINT', 2); [/sql]

POINTは、POINT型にキャストしています。中には緯度と経度が入っていますが、前述した部分でこのような記載をしました。下記部分の、4326と、2という部分に着目してください。

まず、"2"につきましては、このポインタは「X座標とY座標をもつ2次元」の型であることを示しています。"4326"という数字に関しては、WGS84のsridになります。WGS84世界測地系で番号で表すと4326であると覚えておいてください。また国土交通省の街区情報の緯度経度は世界測地系のため、そのような設定にしています。 @作成したテーブル定義を確認する @作成したテーブル定義が無事反映されているか確認します。 [sql] pgis=# \d geos; Table "public.geos" Column   |            Type             |                     Modifiers ------------+-----------------------------+--------------------------------------------------- id         | integer                     | not null default nextval('geos_id_seq'::regclass) prefecture | character varying(64)       | default NULL::character varying city       | character varying(128)      | default NULL::character varying area       | character varying(128)      | default NULL::character varying jiban      | character(64)               | default NULL::bpchar zaban      | integer                     | x          | numeric(7,1)                | default NULL::numeric y          | numeric(7,1)                | default NULL::numeric lat        | numeric(8,5)                | default NULL::numeric lng        | numeric(8,5)                | default NULL::numeric jukyo      | integer                     | daihyo     | integer                     | history16  | integer                     | history17  | integer                     | updated    | timestamp without time zone | not null default now() created    | timestamp without time zone | not null default now() geom       | geometry                    | Indexes: "geos_pkey" PRIMARY KEY, btree (id) Check constraints: "enforce_dims_geom" CHECK (st_ndims(geom) = 2) "enforce_geotype_geom" CHECK (geometrytype(geom) = 'POINT'::text OR geom IS NULL) "enforce_srid_geom" CHECK (st_srid(geom) = 4326) [/sql] ちゃんとできているようですね。^^

・街区情報のインポート 前回と変わらず街区情報を国土交通省からダウンロードし、そのデータをインポートします。 [sql] COPY geos (prefecture,city,area,jiban,zaban,x,y,lat,lng,jukyo,daihyo,history16,history17) from '/tmp/tokyo_utf8.csv' WITH CSV HEADER; [/sql]

・geomフィールドにデータを投入 先ほど作成したgeomカラムに対して値を放り込みます。また、例によってよくわからないSQLを書いてgeomカラムを更新するSQLを生成するという技を使います(もっと、きれいな方法をご存知の方、教えてください!!) [sql]

psql pgis -U postgres -c "select 'update geos set geom = GeometryFromText(''POINT(' || lng || ' ' || lat || ')'',4326) ' || 'where id= ' || id || ';' from geos;" -t > /tmp/geom_update_postgres.sql

[/sql]

・生成されたSQLを確認する 生成されたSQLをちょっと確認してみます。

[sql]

head /tmp/geom_update_postgres.sql

update geos set geom = GeometryFromText('POINT(139.24953 35.73422)',4326) where id= 2; update geos set geom = GeometryFromText('POINT(139.24953 35.73422)',4326) where id= 3; update geos set geom = GeometryFromText('POINT(139.24953 35.73422)',4326) where id= 4; [/sql]

GeometryFromTextは、MySQL同様で文字列表現をジオメトリ型に変換する関数です。

・生成されたSQLを実行 上記で生成した、SQLを、実行します(シェルスクリプト化しろよww)

[sql]

psql pgis -U postgres -f '/tmp/geom_update_postgres.sql'

[/sql]

以上でgeosテーブルはデータを含め、すべて設定されました。。。と思います。

・ためしてみよう!!!SQL では、giosテーブルを利用した、SQLを試してみます。うまくいくかな?ドキドキワクワク・・

検出された緯度と経度で半径100m以内の街区情報を近い順から出力する。 んと・・・ここ重要です。MySQLだと標準では実装されていません。MySQLですと、距離的なジオメトリ関数は存在しませんので独自に関数を作る必要があると思います。逆にPostGISは、距離的なジオメトリ関数が含まれておりますので当然のことながらこれを活用してみます。 うりゃ!!どうだ! [sql] SELECT * FROM ( SELECT prefecture,city,area,rtrim(jiban) as banchi,X(geom) as lng,Y(geom) as lat, distance_spheroid( geom, GeometryFromText('POINT(139.724501 35.687577)',4326), 'SPHEROID["GRS_1980",6378137,298.257222101]' ) AS KYORI FROM geos where distance_spheroid(geom,GeometryFromText('POINT(139.724501 35.687577)',4326),'SPHEROID["GRS_1980",6378137,298.257222101]') < 100 ) AS GISX ORDER BY GISX.KYORI; [/sql]

・・とこんな感じでSQLを書いてみました。緯度経度(139.724501 35.687577)は回答から先にいいますと、四ツ谷の二丁目です。 各SQL文で着目すべき点のみ、記載します。(あまり詳しく理解していませんので、こんなものかといったニュアンスで認識していただければと思います) ○distance_spheroid distance_spheroid(geom,GeometryFromText('POINT(139.724501 35.687577)',4326) 二つの緯度/経度の地点間の距離を,回転楕円体を指定して計算します、上記の例ですと、テーブル上の緯度/経度(geomカラム) から139.724501 35.687577の点間距離をメートル単位での結果を返します。 ○SPHEROID 'SPHEROID["GRS_1980",6378137,298.257222101]' 回転楕円形の情報を与える関数です。いろいろよくわからない文字と数字がありますが、これは地球楕円体の設定で、楕円体記号、赤道の半径(m)、扁平率の逆数(1/f)等の情報です。現在日本の測地系は、GRS_1980が主流のため、そのような設定にしました。もっと的確な説明があればよかったのですが、とにかく「このようなもの」と認識でスルーさせてください。 結果どうなの? 上記SQLを実行した結果は以下のようになりました。 [sql] ------------+--------+------------+--------+-----------+----------+------------------ prefecture |  city  |    area    | banchi |    lng    |   lat    |      kyori ------------+--------+------------+--------+-----------+----------+------------------ 東京都     | 新宿区 | 四谷二丁目 | 12     | 139.72438 | 35.68781 | 28.0764049790079 東京都     | 新宿区 | 四谷二丁目 | 11     | 139.72457 | 35.68731 | 30.2759248844976 東京都     | 新宿区 | 四谷二丁目 | 10     |   139.725 | 35.68763 | 45.5497254663094 東京都     | 新宿区 | 三栄町     | 18     | 139.72524 |  35.6879 | 75.8879717852212 東京都     | 新宿区 | 三栄町     | 25     | 139.72415 | 35.68827 | 83.1962062966614 東京都     | 新宿区 | 四谷二丁目 | 7      | 139.72487 | 35.68687 | 85.2589884379294 東京都     | 新宿区 | 三栄町     | 18     | 139.72509 | 35.68819 | 86.4200439986426 東京都     | 新宿区 | 三栄町     | 24     | 139.72472 | 35.68834 | 86.9472210027482 東京都     | 新宿区 | 四谷二丁目 | 13     | 139.72364 | 35.68802 | 92.1407737733267 東京都     | 新宿区 | 四谷二丁目 | 8      | 139.72534 | 35.68706 | 95.1741204078278 東京都     | 新宿区 | 三栄町     | 17     | 139.72552 | 35.68781 | 95.7920604004985 [/sql]

最も近いところで、(kyori)約28mのところで、 四谷二丁目12番地が見つかりました。 さらなる、パフォーマンスアップに! 上記で期待する結果はでてこれで話を終わりにしようと思ったのですが、パフォーマンスが若干悪い・・・ので、パフォーマンスをもっとよくするSQLを記載方法も追記します。 [sql] SELECT * FROM ( SELECT prefecture,city,area,rtrim(jiban) as banchi,X(geom) as lng,Y(geom) as lat, distance_spheroid( geom, GeometryFromText('POINT(139.724501 35.687577)',4326), 'SPHEROID["GRS_1980",6378137,298.257222101]' ) AS KYORI FROM geos where geom && Box2D(ST_GeomFromText('LINESTRING(139.650 35.600, 139.800 35.750)')) AND distance_spheroid(geom,GeometryFromText('POINT(139.724501 35.687577)',4326),'SPHEROID["GRS_1980",6378137,298.257222101]') < 100 ) AS GISX ORDER BY GISX.KYORI;

pgis-#   GISX.KYORI; prefecture |  city  |    area    | banchi |    lng    |   lat    |      kyori ------------+--------+------------+--------+-----------+----------+------------------ 東京都     | 新宿区 | 四谷二丁目 | 12     | 139.72438 | 35.68781 | 28.0764049790079 東京都     | 新宿区 | 四谷二丁目 | 11     | 139.72457 | 35.68731 | 30.2759248844976 東京都     | 新宿区 | 四谷二丁目 | 10     |   139.725 | 35.68763 | 45.5497254663094 東京都     | 新宿区 | 三栄町     | 18     | 139.72524 |  35.6879 | 75.8879717852212 東京都     | 新宿区 | 三栄町     | 25     | 139.72415 | 35.68827 | 83.1962062966614 東京都     | 新宿区 | 四谷二丁目 | 7      | 139.72487 | 35.68687 | 85.2589884379294 東京都     | 新宿区 | 三栄町     | 18     | 139.72509 | 35.68819 | 86.4200439986426 東京都     | 新宿区 | 三栄町     | 24     | 139.72472 | 35.68834 | 86.9472210027482 東京都     | 新宿区 | 四谷二丁目 | 13     | 139.72364 | 35.68802 | 92.1407737733267 東京都     | 新宿区 | 四谷二丁目 | 8      | 139.72534 | 35.68706 | 95.1741204078278 東京都     | 新宿区 | 三栄町     | 17     | 139.72552 | 35.68781 | 95.7920604004985 [/sql]

結果は同じですが、レスポンスが、5倍くらいよくなりました!私の書いた例ですと、前回MySQLSQLを書いたとおり、全件総当りでチェックしていたのですが、上記のSQLの例ですと、検出拠点は、決まっているので、ある程度レンジを抽出してから、正確な値を絞り込んでいくといったような手法になりますね。 ここで、新しくbox2Dというジオメトリ関数がでてきました。box2Dは2次元の矩形を生成します、前述で、「ある程度レンジを抽出して」の部分がここで利用されているわけです。

おまけ 最寄りの駅を探してみよう!!! もうすこし利用用途のありそうな例をおまけでつけます。今回のケースでは、上記地点(139.724501 35.687577)からもっとも近い駅を検索するSQLを例にしています。なお、設定方法に関する説明は、前述と変わらないため、割愛させていただきます。 ちょっと話はそれますが、無償で利用できるデータがあったりします。 駅データ.jp : http://www.ekidata.jp/st01.html 国土交通省同様、商用/非商用問わず、利用は自由です。ここからデータをダウンロードします(http://www.ekidata.jp/download/index.html) ダウンロードのフォーマットにつきましては、http://www.ekidata.jp/design/index.html をご覧ください。上記リンク先のフォーマットに従い、以下のようなテーブルを作成しました。

<<駅情報テーブル>> [sql] CREATE TABLE railway ( id SERIAL PRIMARY KEY, rr_cd decimal(2,0) NOT NULL, line_cd decimal(5,0) NOT NULL, station_cd decimal(7,0) NOT NULL, line_sort decimal(5,0) DEFAULT NULL, station_sort decimal(7,0) DEFAULT NULL, station_g_cd decimal(7,0) NOT NULL, r_type decimal(1,0) NOT NULL, rr_name varchar(95) NOT NULL, line_name varchar(192) NOT NULL, station_name varchar(192) NOT NULL, pref_cd decimal(2,0) NOT NULL, lng decimal(8,5) DEFAULT NULL, lat decimal(8,5) DEFAULT NULL, f_flag decimal(1,0) DEFAULT NULL, updated timestamp NOT NULL DEFAULT current_timestamp, created timestamp NOT NULL DEFAULT current_timestamp ) ; [/sql] ・geomカラム追加 [sql] SELECT AddGeometryColumn('public', 'railway', 'geom', 4326, 'POINT', 2);

addgeometrycolumn

public.railway.geom SRID:4326 TYPE:POINT DIMS:2 (1 row) [/sql]

・geomカラムの索引生成 [sql] create index idx_railway01 on railway using gist (geom); [/sql]

・テーブルの状態確認 [sql] pgis=# \d railway; Table "public.railway" Column    |            Type             |                      Modifiers --------------+-----------------------------+------------------------------------------------------ id           | integer                     | not null default nextval('railway_id_seq'::regclass) rr_cd        | numeric(2,0)                | not null line_cd      | numeric(5,0)                | not null station_cd   | numeric(7,0)                | not null line_sort    | numeric(5,0)                | default NULL::numeric station_sort | numeric(7,0)                | default NULL::numeric station_g_cd | numeric(7,0)                | not null r_type       | numeric(1,0)                | not null rr_name      | character varying(95)       | not null line_name    | character varying(192)      | not null station_name | character varying(192)      | not null pref_cd      | numeric(2,0)                | not null lng          | numeric(8,5)                | default NULL::numeric lat          | numeric(8,5)                | default NULL::numeric f_flag       | numeric(1,0)                | default NULL::numeric updated      | timestamp without time zone | not null default now() created      | timestamp without time zone | not null default now() geom         | geometry                    | Indexes: "railway_pkey" PRIMARY KEY, btree (id) "idx_railway01" gist (geom) Check constraints: "enforce_dims_geom" CHECK (st_ndims(geom) = 2) "enforce_geotype_geom" CHECK (geometrytype(geom) = 'POINT'::text OR geom IS NULL) "enforce_srid_geom" CHECK (st_srid(geom) = 4326) [/sql]

・ 駅名ダウンロードデータのutf-8化 [shell]

nkf -w /tmp/m_station.csv > /tmp/m_station_utf8.csv

[/shell]

・駅情報のインポート [shell] COPY railway (rr_cd,line_cd,station_cd,line_sort,station_sort,station_g_cd,r_type,rr_name,line_name,station_name,pref_cd,lng,lat,f_flag) from '/tmp/m_station_utf8.csv' WITH CSV HEADER; [/shell]

・geomフィールドの更新SQL作成 [sql]

psql pgis -U postgres -c "select 'update railway set geom = GeometryFromText(''POINT(' || lng || ' ' || lat || ')'',4326) ' || 'where id= ' || id || ';' from railway;" -t > /tmp/railway_update_postgres.sql

[/sql]

・geomフィールドの更新処理 [sql]

psql pgis -U postgres -f '/tmp/railway_update_postgres.sql'

[/sql]

・テーブルのバキューム [sql] vacuum analyze; [/sql]

SQLを発行する 四ツ谷二丁目から一番近い、駅を検索してみます。 [sql] SELECT * FROM ( SELECT rr_name,line_name,station_name,X(geom) as lng,Y(geom) as lat, distance_spheroid( geom, GeometryFromText('POINT(139.724501 35.687577)',4326), 'SPHEROID["GRS_1980",6378137,298.257222101]' ) AS KYORI FROM railway where distance_spheroid(geom,GeometryFromText('POINT(139.724501 35.687577)',4326),'SPHEROID["GRS_1980",6378137,298.257222101]') < 1000 ) AS GISX ORDER BY GISX.KYORI;

--------------+------------------------+--------------+-----------+----------+------------------ rr_name    |       line_name        | station_name |    lng    |   lat    |      kyori --------------+------------------------+--------------+-----------+----------+------------------ 東京メトロ   | 東京メトロ丸ノ内線     | 四谷三丁目   |  139.7201 | 35.68796 | 400.628780565716 東京メトロ   | 東京メトロ南北線       | 四ツ谷       | 139.72955 | 35.68603 | 488.198722616485 東京都交通局 | 都営新宿線             | 曙橋         | 139.72288 |  35.6924 | 554.878205966273 JR           | JR中央・総武線         | 四ッ谷       | 139.73064 | 35.68604 | 581.273897956608 JR           | JR中央線(快速)         | 四ッ谷       | 139.73064 | 35.68604 | 581.273897956608 JR           | JR中央本線(東京0塩尻) | 四ッ谷       | 139.73064 | 35.68604 | 581.273897956608 東京メトロ   | 東京メトロ丸ノ内線     | 四ツ谷       | 139.72995 | 35.68459 | 594.242941765542 JR           | JR中央・総武線         | 信濃町       | 139.72073 | 35.68003 | 904.269714847293 JR           | JR中央線(快速)         | 信濃町       | 139.72073 | 35.68003 | 904.269714847293 (9 rows) [/sql]

最後に PostGISMySQLで大きく違うところは、二点間の距離を、distance_spheroidという関数で比較的容易に算出できるところです。PostGISを利用する理由として、このようなジオメトリ関数が豊富に存在しているところかと思います。点間距離が簡単に使えるので、色々なところで応用できるのではないか?といった例でした。 今回は、PostGIS利用した位置情報検索を例に説明いたしました。実際はもっともっと細かなことができますが、もう限界・・・。詳細につきましては、PostGIS日本語マニュアル 等がありますので、そちらを参考にして、もっとすごいものを作成したほうが良いか・・・と思います。