Web屋のネタ帳

結論はPostgreSQLの幾何データ機能を有効に使おうということである（MySQLでも同様の機能はあるけど）。 とにかくあいかわらず地理情報データとGoogleMaps APIの組み合わせが アツイようなので、その基本的な要件となるであろう、 データベース上の位置情報を効率的に検索する方法を紹介したい。

たとえばおいしいケーキ屋さんの位置情報がデータベース上にあるとしよう。 お店の名前とその緯度経度の情報が数百件あるとする。 GoogleMapsなどである範囲の地図を表示したとして、 お店の位置を地図上にマーキングさせたい場合には、 その地図の範囲の情報をキーにしてデータベース上の緯度経度を検索する必要がある。

普通に考えうるSELECT文は次のとおり。なお、x1,y1は地図の左上の位置、 x2y2は地図の右下の位置とする。

select 店名,緯度,経度 from shop
where 緯度<=y1 and 緯度>=y2 and 経度<=x2 and 経度>=x1

ひとまずこれで要件は満たせるだろう。 だが、ひとたび

• ある1点から半径rの円内に該当するデータを検索したい
• さらにその検索結果を、中心点からの距離でソートしたい

といったことになると、とたんに難しくなる。 しかし、PostgreSQLにもう5年以上前から実装されている 幾何データ型、幾何関数、幾何演算子を使えば、SELECT一発でできることだ。

まずは実験環境を用意する。PostgreSQL8.1.4をインストールし、 2次元座標データ型を持つテーブルを作成。こんな感じ。

create table geodata (id int4, geo point);

pointという聞きなれないデータ型が指定されていることにご注目。 VARCHAR型とかDATE型といった一般的な型しかないOracleなどと違って、 PostgreSQLはこういう多彩なデータ型を利用することができる。

実際にデータを入れる場合はこんな感じになる。

insert into geodata (id, geo) values (1, point(3,4))

x座標で3、y座標で4というデータがはいる。 point()は通常の数値を幾何型に変換する関数。

テストのために、アトランダムな座標を持つデータを10万件ほどつっこんでおく。

testdb=> select * from geodata;
  id   |     geo
-------+--------------
     0 | (726,3654)
     1 | (3804,8386)
     2 | (3644,2121)
     3 | (3079,1080)
(以下省略。全10万レコード)

ちなみにx,y座標を明確に取り出したい場合には

testdb=> select geo[0] as x, geo[1] as y from geodata;
   x   |   y
-------+-------
   726 |  3654
  3804 |  8386

といった感じで可能。

手始めに、ある矩形（長方形のこと）の範囲内にあるレコードを検索してみる。

testdb=> select * from geodata where geo @ box(point(100,100),point(200,200));
  id   |    geo
-------+-----------
 41282 | (198,151)
 44787 | (105,103)
 91065 | (179,104)
 96581 | (125,100)
 97615 | (112,142)
(5 rows)

成功。@は「含む、または境界上」という演算子。 boxは二つのpointを対角とした長方形を示す。

次に、ある円の範囲内にあるレコードを検索してみる。

testdb=> select * from geodata where geo @ circle(point(300,300),40);
  id   |    geo
-------+-----------
  5501 | (289,336)
 20403 | (270,306)
 43715 | (319,265)
 74237 | (274,293)
 84113 | (326,327)
(5 rows)

circle(point(300,300),40) がミソ。これはx=300,y=300の座標を中心とする半径40の円を表す。

さらに、上の例で、中心からの距離で結果をソートしてみる。

testdb=> select * from geodata where geo @ circle(point(300,300),40) order by geo <-> point(300,300);
  id   |    geo
-------+-----------
 74237 | (274,293)
 20403 | (270,306)
 84113 | (326,327)
  5501 | (289,336)
 43715 | (319,265)
(5 rows)

疑り深い人は↓

testdb=> select id, geo, geo <-> point(300,300) as dist from geodata 
where geo @ circle(point(300,300),40) order by geo <-> point(300,300);
  id   |    geo    |       dist
-------+-----------+------------------
 74237 | (274,293) | 26.9258240356725
 20403 | (270,306) | 30.5941170815567
 84113 | (326,327) | 37.4833296279826
  5501 | (289,336) | 37.6430604494374
 43715 | (319,265) | 39.8246155034798
(5 rows)

さて気になるのは性能だ。インデックスは効くのだろうか？

testdb=> create index testindex on geodata (geo);
ERROR:  data type point has no default operator class for access method "btree"
HINT:  You must specify an operator class for the index or define a default operator class for the data type.

おっと。インデックスを張れない！ インデックスの種類の問題かと思ったが

testdb=> create index testindex on geodata using gist (geo);
ERROR:  data type point has no default operator class for access method "gist"
HINT:  You must specify an operator class for the index or define a default operator class for the data type.
testdb=> create index testindex on geodata using rtree (geo);
ERROR:  data type point has no default operator class for access method "rtree"
HINT:  You must specify an operator class for the index or define a default operator class for the data type.

デフォルトのbtreeはもちろんgistもrtreeもだめ。 こいつは困ったなんでだろう？

実は、現在のPostgreSQLでは、インデックスを張ることができるのは box型(長方形)、circle型（円）、polygon型(多角形)の3種類だけ。 point型にはインデックスを張れない。 マニュアルのここの下のほうにあるSELECTを投げて調べてみたが、 幾何データ型と思われるものは上の3つの演算子についてしか見当たらない。 つまり、「boxをboxでインデックス検索」「円を円でインデックス検索」 「多角形を多角形でインデックス検索」ということはできても、 「円で点をインデックス検索」といったことができないのだ。 なんでそんなことになっているのかは知らないがとにかくそうなっている。 もちろん将来のバージョンで改善されることを期待したいところだ。 DBにおけるインデックスの内部構造とC言語に自信のある人はpoint型用のgist演算子クラスの実装にチャレンジしてみてはいかがだろう。

しかし、現在のPostgreSQLであってもちょっとした方法で解決できる。 PostgreSQLには関数インデックスという便利なインデックスの張り方がある。 値に対してそのままインデックス化するのではなく、 何らかの一定の関数を通した結果をインデックスしておくというものだ。 この機能と、「点は、半径ゼロの円である」「点は、4点全てが同じ長方形である」 といった幾何学的トリックを組み合わせればいい。

testdb=> create index testindex on geodata using gist(circle(geo,0));
CREATE INDEX
testdb=>

これでgeo列（point型）に対してそれを半径ゼロのcircleとする関数をかませて circle型とし、それに対してインデックスが張れたことになる。 検索する場合には次のようにする。

testdb=> select * from geodata where circle(geo,0) @ circle(point(300,300),40);
  id   |    geo
-------+-----------
  5501 | (289,336)
 20403 | (270,306)
 43715 | (319,265)
 74237 | (274,293)
 84113 | (326,327)
(5 rows)

OK。さっきと結果は同じ。

では、インデックスを張ったことでどの程度の効果が出たのだろうか。 先ほどのインデックスを一旦dropして、explainで速度を量ってみる。

testdb=> drop index testindex;
DROP INDEX
testdb=> explain analyze select * from geodata where circle(geo,0) @ circle(point(300,300),40);
                                                QUERY PLAN
----------------------------------------------------------------------------------------------------------
 Seq Scan on geodata  (cost=0.00..2137.00 rows=100 width=20) (actual time=21.102..381.205 rows=5 loops=1)
   Filter: (circle(geo, 0::double precision) @ '<(300,300),40>'::circle)
 Total runtime: 404.447 ms
(3 rows)

当然ながらSeq Scan（全走査）されている。インデックスを張った状態で同じSELECTをするとどうなるか。

testdb=> create index testindex on geodata using gist(circle(geo,0));
CREATE INDEX
testdb=> explain analyze select * from geodata where circle(geo,0) @ circle(point(300,300),40);
                                                     QUERY PLAN
--------------------------------------------------------------------------------------------------------------------
 Bitmap Heap Scan on geodata  (cost=2.35..268.61 rows=100 width=20) (actual time=0.720..0.781 rows=5 loops=1)
   Filter: (circle(geo, 0::double precision) @ '<(300,300),40>'::circle)
   ->  Bitmap Index Scan on testindex  (cost=0.00..2.35 rows=100 width=0) (actual time=0.654..0.654 rows=5 loops=1)
         Index Cond: (circle(geo, 0::double precision) @ '<(300,300),40>'::circle)
 Total runtime: 0.978 ms
(5 rows)

インデックススキャンが走るようになった。実行時間も、404ミリ秒から0.9ミリ秒へ。 実に400倍のスピードアップである。

このように、PostgreSQLの幾何データ機能は有効に使えるものだ。 詳しくはマニュアルを。

• 幾何データ型
• 幾何関数と演算子

see also:
Google Maps とIEの「開けません。 操作は中断されました」(2006.4)
Google Maps APIと文字コード(2005.9)
地図情報と郵便番号(2005.6)