PostGISのGeospatial化大作戦 その３（デタラメ編）

PostGISは地球を球体(sphere)または回転楕円体(spheroid)とみなした上での、地表面距離を計算して距離を計算するための関数（distance_sphere、distance_spheroid）を持っています。
なので、

SELECT * FROM geo_table WHERE distance_sphere(geo_column, GeomFromText('POINT(140 89)',4326)) < 10000;

みたいなSQLで、経緯度で位置情報データが登録されたテーブルから周辺10kmのデータだけを検索する事ができます。
ですがこの検索、空間インデックスが利用されない（総当り検索になる）のが悩みの種で。

もっとも、そもそも地表面距離を計算する関数でなくとも、格納している座標系上での距離を計算する関数distance（主に、経緯度座標系ではなく、メートル単位での座標系となる横メルカトル図法とか、そういった地図投影済みの地理データを格納している場合に使います）でも、空間インデックスは利用されません。
基本的に、空間インデックスというのは、2次元または3次元座標系での、矩形どうしの重なり合いまでしか対象にできず、それ以上詳細な条件に対する絞込みは、まず矩形を使った空間インデックスでの絞り込みの後に、詳細条件でさらに絞り込む、という手順になります。
よって、例えばメートル単位の横メルカトル図法での座標系で座標(100,100)の周辺10kmのデータを検索する場合、こんな感じになると思います。

SELECT * FROM tm_table WHERE expand(GeomFromText('POINT(100 100)',2447),10000) && tm_column AND distance(tm_column, GeomFromText('POINT(100 100)',2447) ) < 10000;

※　expand関数は指定した座標系単位での周辺第2引数分の矩形を自動的に作ってくれる関数、2447は日本の平面直角座標系第５系のEPSGコード。また&&演算子は矩形が重なり合いを持つかの条件判断を空間インデックスを利用して行ってくれる演算子。

説明を加えると、まずexpand関数を用いて発生させた周辺10kmの矩形と座標カラムの重なりを空間インデックスを使って抽出後、抽出されたレコードに対して、総当りで距離を測定し10km以内の条件をチェックしている、という手順になっています。

よって、distanceに対して絞込みのためにexpand関数があるように、distance_sphereやdistance_spheroidに対しても、expand_sphereやexpand_spheroidのような関数を準備してやれば、空間インデックスを活用できるのではないか？というのが、本エントリのその１の主旨。
でも球面や回転楕円体面で一定距離離れた円の接する経緯度の導出が思いの外ムズカシスなので、とりあえずインデックスの効果だけ確かめたのが本エントリその２の内容でした。

んでその後、全くその方面追求してなかったんだけど、Web2.0ワークショップで講演することになったので、流石にインデックスも使えない周辺検索の方法を涼しい顔して紹介するわけにもいかんだろうと。
かといって厳密な球面や回転楕円体面で一定距離離れた円の接する経緯度の導出もそう簡単にはできんので、どうしようかと悩んだ挙句、日本近辺での利用だけでお茶を濁そうと、そのような結論にいたりました。

理科年表(第76冊：2003年版，p.562)によると、経度方向及び緯度方向に1秒移動する弧の長さは、緯度によって異なるがそれぞれだいたい20m以上、30m以上であると。
ということは、日本国内に限れば、中心点から一定距離：x m離れた円は、経緯度で経度方向に±x/20秒、緯度方向に±x/30秒離れた経緯度座標系上での矩形の中に含まれていると考えてよいのではないかと。
（正確には、緯度方向はともかく経度方向の弧は、中心点で正射投影した場合の線が直線ではなく弧になるので、本当にこの矩形の中に絶対望む円が含まれるかは判らないけど...その辺はマージンとってるし、まあ誤差ありありの手法という事で、いいのではないかと。）
というわけで、そういう考え方で作成した関数・expand_sphere_pseudoのplpgsqlソースをさらします。

CREATE FUNCTION expand_sphere_pseudo(geometry,double precision)
RETURNS geometry
AS '
    DECLARE
        fx double precision;
        fy double precision;
        dx double precision;
        dy double precision;
        se varchar(11);
        sw varchar(11);
        ss varchar(11);
        sn varchar(11);
        sid integer;
        wkt text;
        geo geometry;
    BEGIN
        fx = x($1);
        fy = y($1);
        sid = SRID($1);
        dx = $2/20.0/3600.0;
        dy = $2/30.0/3600.0;
        se = to_char(fx-dx,''S999D999999'');
        sw = to_char(fx+dx,''S999D999999'');
        ss = to_char(fy-dy,''S999D999999'');
        sn = to_char(fy+dy,''S999D999999'');
        wkt = ''LINESTRING('' || se || '' '' || ss || '',
              '' || sw || '' '' || sn || '')'';
        geo = GeomFromText(wkt,sid);
        RETURN geo;
    END
' IMMUTABLE LANGUAGE plpgsql;

ここでIMMUTABLEを忘れると、関数の実行結果がキャッシュされず総当りで評価されるので、かえって遅くなってしまうので注意（ちょっとハマった）。

上記関数を用いて、その２記事でも使った10万件の位置情報データに対し、東経135度、北緯35度からの半径100kmの円内検索を実行したところ、

SELECT * FROM geo_table WHERE distance_sphere(geo_column, GeomFromText('POINT(135.0 35.0)',4326)) < 100000;

では525msの検索時間がかかったのに対し、

SELECT * FROM geo_table WHERE expand_sphere_pseudo(GeomFromText('POINT(135.0 35.0)',4326),100000) && geo_column AND distance_sphere(geo_column, GeomFromText('POINT(135.0 35.0)',4326)) < 100000;

では35msと、10倍以上の高速化が実現できました。

タイトルにも書いたとおり、所詮日本近傍でしか使えないデタラメテクニックですが、一応結果を出せたので報告をば。