2007年10月17日

PostGIS応用例：オンザフライの経緯度データが特定ポリゴンに含まれるかを調べる方法

マッシュアップアワードで賞を取られていた、ふむふむソフトさんの手書き地図検索：CHIZTEKというサービスがあります。
地図画面上でマウスで手書きポリゴンを描くと、そのポリゴン内に含まれるお店等を検索してくれるという、直感的で面白いサービスです。

▲ マウスで書いたポリゴン内のお店検索 ▲

ポリゴンに内包されるPOIの検索、というのは極めてGIS的な発想なので、このサービスを見た時はすごく興奮しました。
どうやって実現しているのか、と聞いてみると、POI（Point of Interesting：興味の対象となるスポットの意）自体はポリゴンの最小内接矩形（MBR）内に含まれるものを全部取ってきて、クライアントサイドのJavaScriptでポリゴン内に入っているか否かを判定するアルゴリズムを組み、フィルタリングをしているそうです。
JavaScriptが使いこなせていない事もありますし、その手の本格的なアルゴリズムの実装が苦手な事もあって、私ではとてもできない技術力です。
すごい。

でも、サーバサイドで、かつ自力実装ではなくツールを使ってよければ、与えられたポリゴンとPOI群に対し、ポリゴン内に含まれるPOIをフィルタリングすることは比較的簡単に可能です。
PostGISが、その問題を解決してくれます。

PostGISが、DB内に蓄積された位置情報データについて、ポリゴンへの内包等の計算をして結果を出してくれることについてはよく知られていると思いますが、PostGISのコマンドもしょせんはSQLなので、DB内への問い合わせだけではなく、

SELECT 1+2;

みたいな計算もSQLはオンザフライで結果を返してくれるのと同様、地理情報計算もオンザフライでやってくれます。

例えば、どこかのWeb APIを通じて、

経度136度緯度33度
経度136度緯度34度
経度140度緯度50度
経度135度緯度40度

の4つのPOIを得たとします。
このPOIのうち、

経度135度緯度35度、経度137度緯度30度、経度139度緯度35度

の3点からなる三角形ポリゴンの中に含まれるPOIだけを抜き出すには、以下のようにすればOKです。

SELECT AsText(Intersection(GeometryFromText('POLYGON((135 35,137 30,139 35,135 35))',4326),
GeometryFromText('MULTIPOINT(136 33,136 34,140 50,135 40)',4326)) );

結果：MULTIPOINT(136 33,136 34)

で、２つのPOIだけが三角形の中に含まれることがわかります。
やっている事は、POLYGON（多角形）とMULTIPOINT（複数の点の集合体）の交わる図形を求める関数（Intersection）を求めているのですが、その結果は当然MULTIPOINTになるはずなので、その戻ってきたMULTIPOINT内の座標をピックアップしてやれば、ポリゴン内に含まれる点が判る、ということです。

ただしこのやり方で注意して欲しいのは、PostGISは座標の扱いは飽くまで曲率のない直角直交座標として扱われるので、ここで内包判定に用いられるポリゴンの辺は、頂点間を結ぶ大円上の弧とはならない（つまり楕円体上での最短経路とはならない）、ということです。
なので、一つの街の中でのポリゴン内の内包計算なんかだと実質上問題になるような誤差は生じないでしょうが、上の例のような緯度経度が1度以上変わるようなでっかい領域での計算だと、望むような結果は得られません。

また別の問題として、この方法だとSQLの中に含められるPOIのデータは経緯度データだけとなるため、このSQLの実行結果とPOIのその他の属性を紐付けるのは経緯度データだけ、ということになります。
なので、SQLに与えた経緯度データとSQLの結果で出てくる経緯度データが、有効数字の差等で変わってしまうと、うまく紐付けられないことになってしまいます。

例えば、

SELECT AsText(Intersection(GeometryFromText('POLYGON((135 35,137 30,139 35,135 35))',4326),
GeometryFromText('MULTIPOINT(136.00 33.00,136.00 34.00,140.00 50.00,135.00 40.00)',4326)) );

結果：MULTIPOINT(136 33,136 34)

みたいな形になってしまうと、うまく紐付けられないですよね。

そんな時は、以前の記事で紹介した、拙作のPostGISでLocapointを使う関数を使ってみて下さい。

SELECT AsLocapoint(Intersection(GeometryFromText('POLYGON((135 35,137 30,139 35,135 35))',4326),
GeometryFromLocapoint('MULTIPOINT(RT9.WV3.IR4.UF8,RX6.WV3.XC9.UF8,UF7.XC8.UF8.XC9,SU2.WT5.FU3.AA0)')) );

結果：MULTIPOINT(RT9.WV3.IR4.UF8,RX6.WV3.XC9.UF8)