さて、先日（と言っても結構前）、地理情報をPostgreSQLで扱う例として、巡回セールスマン問題をPostgreSQLで解きつつGoogle Mapsで可視化するエントリを書きました。
- 巡回セールスマン問題における最短経路をpgRoutingで探索する
http://pgsqldeepdive.blogspot.jp/2016/09/pgrouting.html
実現したいことは、
- 観光に関連する情報をPostgreSQLに取り込んで、
- Google Mapsで地図上にマッピングして可視化しつつ、
- 地図上をブラウジングしながら、
- 興味のある場所があったらそのままGoogle検索に飛ぶ
年末年始のお出かけの検討に、または雑談のお供にご活用いただければと思います。
■オープンデータ「国土数値情報 観光資源データ」とは
まず、今回使うデータですが、国土交通省が公開している「国土数値情報」の中から「観光資源」のデータを使います。
- 国土数値情報ダウンロードサービス
http://nlftp.mlit.go.jp/ksj/index.html - 国土数値情報 観光資源データの詳細
http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P12-v2_2.html
- 観光資源_ID
- 観光資源名
- 都道府県コード
- 行政コード
- 種別名称
- 所在地住所
- 観光資源分類コード
- 観光資源（地理情報）
このデータは地理情報のデータフォーマットとして広く使われている「シェープファイル（Shape File）」と呼ばれる形式で配布されています。
PostGISには、このシェープファイルをPostgreSQLに取り込むためのコマンドラインツールが含まれており、比較的容易にPostgreSQLのデータベースに取り込むことができます。
- シェープファイルのデータをインポートしてみよう コマンドライン編 - Qiita
http://qiita.com/yellow_73/items/f2d8388e88d534066bc0
■データを準備する
データセットをダウンロードすると分かりますが、47都道府県×3種類で、100以上のシェープファイルが含まれています。
シェープファイルを取り込む shp2pgsql コマンドは1シェープファイル1テーブルとして作成しますので、普通に取り込むと100個以上のテーブルができることになります。前処理としてそれらのテーブルを統合しておいた方が良いでしょう。
というわけで、取り込むスクリプトは以下です。
各シェープファイルの取り込みと、それらを統合して「p12_14」というテーブルを作成してくれます。
実行すると以下のようになります。
[snaga@localhost postgis_test]$ ./004_P12-14_GML.sh Archive: ./data/P12-14_GML.zip creating: P12-14_GML/ inflating: P12-14_GML/KS-META-P12_14-01.xml inflating: P12-14_GML/KS-META-P12_14-02.xml inflating: P12-14_GML/KS-META-P12_14-03.xml ... DROP TABLE DROP TABLE DROP TABLE COMMIT 都道府県コード | count ----------------+------- 01 | 53 02 | 524 03 | 513 04 | 8 05 | 679 ... 45 | 265 46 | 13 47 | 198 (47 rows) [snaga@localhost postgis_test]$ psql -c 'select count(*) from p12_14' gistest count ------- 19140 (1 row) [snaga@localhost postgis_test]$
ついでに都道府県コードの変換テーブルも作成しておきます。
[snaga@localhost postgis_test]$ psql -f 201_PREFCODE.SQL gistest
psql:201_PREFCODE.SQL:1: ERROR: table "prefcode" does not exist
CREATE TABLE
COPY 47
[snaga@localhost postgis_test]$ psql -c 'select count(*) from prefcode' gistest
count
-------
47
(1 row)
[snaga@localhost postgis_test]$
データをロードし終わると、以下のような状態になります。
gistest=> \d
List of relations
Schema | Name | Type | Owner
--------+-------------------+-------+----------
public | geography_columns | view | postgres
public | geometry_columns | view | postgres
public | p12_14 | table | snaga
public | prefcode | table | snaga
public | raster_columns | view | postgres
public | raster_overviews | view | postgres
public | spatial_ref_sys | table | postgres
(7 rows)
gistest=> \d p12_14
Table "public.p12_14"
Column | Type | Modifiers
--------------------+----------------+-----------
観光資源ID | integer |
観光資源名 | text |
都道府県コード | character(2) |
行政コード | character(5)[] |
種別名称 | text |
所在地住所 | text |
観光資源分類コード | numeric(2,0) |
geom | geometry |
Indexes:
"p12_14_idx" btree ("都道府県コード", "観光資源ID")
gistest=> \d prefcode
Table "public.prefcode"
Column | Type | Modifiers
----------------+--------------+-----------
都道府県コード | character(2) | not null
都道府県名 | text | not null
Indexes:
"prefcode_pkey" PRIMARY KEY, btree ("都道府県コード")
gistest=>
■ポリゴンデータを点データに変換する
次に、geomカラムに含まれているポリゴンデータを点のデータに変換します。
ポリゴンデータのままだと、地理情報を使った演算に時間がかかってしまいます。今回はそこまで厳密な距離の計算などは必要ないため、計算量を減らすためにポリゴンから「点」、具体的には「ポリゴンの重心」に変換します。
重心を求めるには PostGIS の ST_Centroid() 関数を使います。
gistest=> select pg_total_relation_size('p12_14');
pg_total_relation_size
------------------------
10821632
(1 row)
gistest=> update p12_14 set geom = ST_Centroid(geom);
UPDATE 19140
gistest=> vacuum full p12_14 ;
VACUUM
gistest=> select pg_total_relation_size('p12_14');
pg_total_relation_size
------------------------
3571712
(1 row)
gistest=>
ポリゴンから点に変換したことによって、テーブルサイズが1/3程度になりました。
■観光資源データを検索する
ここまでできたら観光資源データを検索してみます。
今回の検索は基本的に「緯度経度で矩形を指定し、その中に存在している観光資源の情報を取得する」というものです。というのは、最終的にはGoogle Mapsで表示、検索できるようにしますので、Google Mapsの表示領域である矩形に切り取れる必要があるのです。
クエリとしては、まず検索したい対象の区域の緯度経度からポリゴンを作成し、そのポリゴンに含まれる観光資源を取得します。
ポリゴンを作成するには ST_MakePolygon() 関数を、作成したポリゴンに座標系を設定するには ST_SetSRID() 関数を使います。
そして、ST_Contains() 関数を使って、作成したポリゴンに観光資源が含まれるかどうかを判定します。
例えば、都庁から浜離宮恩賜庭園にかけての一帯を検索したいのであれば、それぞれの座標は
- 都庁: 緯度 35.689634 経度 139.692101
- 浜離宮恩賜庭園: 緯度 35.6597374 経度 139.7634925
この辺。
矩形を指定する方法ですが、「4つの頂点を指定する」のではなく、「頂点を結ぶ線を定義する」必要があります。つまり、最後には「開始点に戻ってくる」必要がある、ということです。開始点と終了点が一致して閉じられていないとエラーが出ますので注意してください。
具体的なクエリとしては以下のようになります。
gistest=> SELECT
観光資源名,
ST_Y(geom) lat,
ST_X(geom) lon,
所在地住所
FROM
p12_14
WHERE
ST_Contains(
ST_SetSRID(ST_MakePolygon('LINESTRING(139.692101 35.689634, 139.692101 35.6597374,
139.7634925 35.6597374, 139.7634925 35.689634,
139.692101 35.689634)'::geometry), 4612),
geom
);
観光資源名 | lat | lon | 所在地住所
----------------------------------+------------------+------------------+-----------------------
江戸城跡 | 35.6848450646745 | 139.753409059839 | 千代田区千代田
江戸城跡 | 35.6848545616798 | 139.753620448968 | 千代田区千代田
原宿 | 35.6677931758621 | 139.706963201124 | 渋谷区神宮前
明治神宮 | 35.6759794856927 | 139.699423688458 | 渋谷区代々木神園町1-1
新宿御苑 | 35.6852930986455 | 139.71005110951 | 新宿区
国会議事堂 | 35.6758880055466 | 139.744858001005 | 千代田区永田町1-7-1
根津美術館 | 35.6622430012588 | 139.717093001577 | 港区南青山6-6-5-1
江戸前の寿司 | 35.6896339956514 | 139.692101001444 | 新宿区
国立能楽堂で上演される能・狂言 | 35.6804039978788 | 139.708209998586 | 渋谷区千駄ヶ谷4-18-1
国立劇場で上演される歌舞伎・文楽 | 35.6815600003142 | 139.74327700092 | 千代田区隼町4-1
明治神宮 | 35.6757916030472 | 139.699511838812 | 渋谷区代々木神園町1-1
新宿御苑 | 35.6850607466651 | 139.709971187167 | 新宿区
(12 rows)
gistest=>
東京にこれしか観光資源ないのかよというツッコミはあろうかと思いますが、PostgreSQLのせいではないのでここでは不問とします。
■検索をREST API化する
矩形で検索できるようになったら、これをブラウザのJavaScriptから呼び出せるようにREST API化します。
Flaskを使ってさくっとREST API化して手元なりどこかのPaaSなりで動作させておきます。
ソースコードは以下です。
■Google Mapsと連携する
最後にGoogle Mapsと連携させます。
（力尽きたので詳細は割愛。ソース見てちょ） 特長、というか、こんな感じで動いてます。
- Google Maps JavaScript API使ってます。
- 表示している領域の四隅の緯度経度を取得して、その領域内にあるアイテムを検索・表示しています。
- アイテムをクリックして情報ウィンドウを表示させ、そこからGoogle検索することができます。
- 表示領域が変わると、新たな座標をパラメータにしてREST APIを呼び、表示すべき項目を取得します。
- 広域表示をすると表示する項目が多くなりすぎるので、最大100件に間引いてます。
- 雑に間引いているので、ズームアウトすると見えなくなったり見えたりします。
- 観光マップ feat. オープンデータ
http://snaga-lab.mybluemix.net/static/maps/index.html
日本全域表示。最大100件に間引いて表示してます。（北海道に何も表示されていませんが、これは雑に間引いて表示している偶然の産物であり他意はございません。ズームしていくと見えてきます）
■まとめ
というわけで、今回はオープンデータをPostgreSQL/PostGISに投入して、クエリをREST API化することによってGoogle Mapsと連携できることを示しました。
先日の巡回セールスマン問題のエントリではMy MapsにKMLファイルをインポートする方式を取りましたので、扱っているデータはあくまでもstaticなデータでしたが、今回はREST APIを使って位置情報を動的に取得して可視化することができることを示しました。
フロントエンドと動的に連携できるようになりましたので、キーワード検索や条件検索などなど、「データベースとつながっていることの価値」が出せるようになるのではないかと思います。
また、今回は時間の都合上確認できませんでしたが、近隣にあるデータをクラスタリングして代表点を求める、みたいなこともできるのではないかと考えています（AVG() GROUP BYのようなノリで）。それができると、広域表示の時にもう少し分かりやすい表示になるように思います。
PostgreSQLの強みの一つにはGISデータの扱いであり、かつ、GISデータを扱えるというだけではなく、先のエントリで紹介したような巡回セールスマン問題のソルバーのようなライブラリが存在している、ということも非常にユニークなところだと思います。
地理情報は、まだまだいろんな使い方ができるのではないかなーと感じています。
興味のある方は、ぜひこれを機会に何かチャレンジしてみていただければと思います。
では、また。
0 件のコメント:
コメントを投稿