それはたぶん計算のしすぎ

はじめに

Advent Calendarがそこここで実施されているなか、id:meketen:20121209 で「Raspberry PiでPostgreSQL(PostGIS)を動かして、現在地をしゃべらせてみた。」というエントリーに遭遇したのは、なんとなくPostGISについてぐぐってた時でした。
PostgreSQL Advent Calendar (http://atnd.org/events/34052)の一環だそう。

安価なLinuxが走るワンボードマイコンでPostGIS入れてごにょごにょというなかなか面白い話。
最後に「オチ」として、応答が非常に遅いというお話でしたので、ここまで変態なことしたんだから屋外に出てしゃべらせるところまで行って欲しいと願ってやみません。

FOSS4G Advent Calendar (http://atnd.org/events/34052) の管理人の人がその日の夕方にはスカウトしていたというので、じゃああれだ、ちょっと速度が上がるか試してみて下さい、ということで、このエントリーを書いて見ました。

傍からあれこれ言うだけで、実際にやられる方の身になっていません、すみません。

PostGISは経度・緯度の順

まずひとつチェック。PostGISでは、座標の表現は全て 東西軸 南北軸 の順にします。
POINT(lat lng) とするのは間違いで POINT(lng lat) とします。
これ重要。後述のGEOGRAPHYでは POINT(35 135) は、東経35度 北緯135度はありえん、といってエラーにしてくれます。

インデクスを効かせて絞り込むのが重要

時間がかかる根本的原因は、距離計算です。距離計算自体が、cos, sin とかを使ううえ、収束計算をさせていて、なるたけしたくないものです。

空間インデクスは id:yellow_73:20110820#p2 あたりにある通り、MBRベースで検索します。ので、絞り込みとしては粗雑ですが距離計算と比べたら計算が非常に速い。

また、2次元空間を扱う場合には、スカラ値を扱う場合よりも、絞り込みをより吟味する必要があります。スカラ値ならヒットする数は1次のオーダですが、2次元だと2次のオーダになるからです。

よって、WHEREでインデクスを効かせて十分に対象ジオメトリを絞り込んだうえで距離計算をした方が良いです。

直交座標系なら簡単

たとえば直交座標系上なら話は簡単で、1万メートルの範囲内にあるかどうかをテストするのは

SELECT ... WHERE ST_DWithin(geom, 10000);

となります。裏でインデクスが貼られているか確かめたうえで、この場合には「半径1万メートルの円に外接するMBR内にあるポイント」、言い換えると、2万メートル四方の正方形内にあるポイント、をインデクスで探してくれます。そのうえで距離計算を行います。

が、経度緯度を扱う場合にはそのようにはいきません。GEOMETRY型で経度緯度をとっても、自動的にはメートルに変換してくれません。
ST_DWithin(geom, 10000)は「半径1万度の円に外接する四角形範囲」となります。1万度て何言ってるかわかりませんね。

GEOMETRY && GEOMETRY を使って絞り込みも可能ですが、第2引数の作成が大変になります。

GEOGRAPHY型を使うといろいろ面倒が無い

今回のポイントデータは経度緯度で管理しているので、メートルから度に変換しなければなりません。てっとり早く経度緯度のデータの絞り込みにST_DWithinをメートル単位で使いたい場合には、GEOGRAPHY型を使うと便利です。

GEOGRAPHYを使うと距離や面積を経度緯度を使っていても、メートル単位で計算してくれます。

db=# SELECT ST_Distance('SRID=4326;POINT(0 0)'::GEOMETRY, 'SRID=4326;POINT(90 0)'::GEOMETRY);
st_distance

• • • • • • • • • • • • -

90
(1 行)

db=# SELECT ST_Distance('SRID=4326;POINT(0 0)'::GEOGRAPHY, 'SRID=4326;POINT(90 0)'::GEOGRAPHY);
st_distance

• • • • • • • • • • • • • • • • • -

10018754.1701286
(1 行)

db=# SELECT ST_Distance_Spheroid('SRID=4326;POINT(0 0)'::GEOMETRY, 'SRID=4326;POINT(90 0)'::GEOMETRY,'SPHEROID["GRS_1980",6378137,298.257222101]]');
st_distance_spheroid

• • • • • • • • • • • • • • • • • • • • • -

10018754.1701286
(1 行)

ST_Distance(GEOMETRY,GEOMETRY)の場合は90(度)を返しましたが、ST_Distance(GEOGRAPHY,GEOGRAPHY)とST_Distance_Spheroidは、メートル単位の距離を返してくれます。

GEOGRAPHY型を使うとインデクスもうまいことしてくれる

ST_Distance(GEOGRAPHY,GEOGRAPHY)とST_Distance_Spheroidとが同じだけなら、計算コストは全く減りません。

しかし、インデクスもうまいことしてくれます。
ST_DWithin(GEOMETRY,GEOMETRY,DOUBLE PRECISION) での空間インデクス検索条件は、「半径$3(単位は$1,$2の定義による)の円に外接するMBR」でしたが、
ST_DWithin(GEOGRAPHY,GEOGRAPHY,DOUBLE PRECISION) では、「半径$3メートルの円に外接するMBR」となります。

しぼりこみ範囲

街区レベルなら、たぶん、500m程度あれば大丈夫じゃないかと思います。

実験

• • テーブル作成
  • geog GEOGRAPHYを追加

CREATE TABLE address(
address_no bigserial PRIMARY KEY,
pref_name text,
city_name text,
street_name text,
address text,
point_no SMALLINT,
x_point decimal(8,1),
y_point decimal(8,1),
lat decimal(8,5),
lng decimal(8,5),
jyukyo_flg bool,
daihyo_flg bool,
upd_bef_flg SMALLINT,
upd_aft_flg SMALLINT,
geom GEOMETRY(POINT,4326),
geog GEOGRAPHY(POINT)
);

CREATE INDEX ix_address_geom ON address USING GiST(geom);
CREATE INDEX ix_address_geog ON address USING GiST(geog);

• • コピー

COPY address(pref_name,city_name,street_name,address,point_no,x_point,y_point,lat,lng,jyukyo_flg,daihyo_flg,upd_bef_flg,upd_aft_flg)
FROM 'C:\Users\Public\Documents\pgtest\14000-10.0a\14_2011.csv' CSV HEADER;

• • lng, lat の順でPOINT作成

UPDATE address SET geom=ST_GeomFromText('POINT('||lng||' '||lat||')',4326);

• • geogを作成

UPDATE address SET geog=('POINT(' || lng || ' ' || lat || ')')::GEOGRAPHY;

• • 本当はgeomのSRIDが4326の場合(今回の場合)には
  • UPDATE address SET geog=geom;
  • でOKで、勝手に変換してくれます。

まずは、こちらはデスクトップ機でやってるので、ラズベリーパイとの比較のために、ひとつ同じクエリを実行してみます。
ちなみに、だいたい123km×87kmの範囲を検索しています。

explain analyze SELECT
pref_name,city_name,street_name,address,
st_distance_spheroid(geom, ST_GeomFromText('POINT(139.642403 35.448575)',4326), 'SPHEROID["GRS_1980",6378137,298.257222101]')
FROM address WHERE
geom && 'BOX3D(140.0 35,138.630066 35.791083)'::BOX3D AND
st_distance_spheroid(geom, ST_GeomFromText('POINT(139.642403 35.448575)',4326),'SPHEROID["GRS_1980",6378137,298.257222101]') < 10000
ORDER BY 5 LIMIT 1
;

• • • • • • • -

Limit (cost=298497.27..298497.27 rows=1 width=173) (actual time=1720.901..1720.901 rows=1 loops=1)
-> Sort (cost=298497.27..299145.52 rows=259300 width=173) (actual time=1720.898..1720.898 rows=1 loops=1)
Sort Key: (st_distance_spheroid(geom, '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geometry, 'SPHEROID("GRS_1980",6378137,298.257222101004)'::spheroid))
Sort Method: top-N heapsort Memory: 17kB
-> Seq Scan on address (cost=0.00..297200.77 rows=259300 width=173) (actual time=84.135..1704.007 rows=58467 loops=1)
Filter: *1
Rows Removed by Filter: 359474
Total runtime: 1720.929 ms

半径10km(20km四方)を範囲とした場合は、次の通り。

explain analyze SELECT
pref_name,city_name,street_name,address,
ST_Distance(geog, ST_GeogFromText('POINT(139.642403 35.448575)'))
FROM address WHERE
ST_DWithIn(geog, ST_GeogFromText('POINT(139.642403 35.448575)'), 10000)
ORDER BY 5 LIMIT 1
;

• • • • • • • -

Limit (cost=11.80..11.81 rows=1 width=269) (actual time=578.037..578.037 rows=1 loops=1)
-> Sort (cost=11.80..11.81 rows=1 width=269) (actual time=578.035..578.035 rows=1 loops=1)
Sort Key: (_st_distance(geog, '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography, 0::double precision, true))
Sort Method: top-N heapsort Memory: 17kB
-> Index Scan using ix_address_geog on address (cost=0.00..11.79 rows=1 width=269) (actual time=3.543..559.561 rows=58467 loops=1)
Index Cond: (geog && '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography)
Filter: *2 AND _st_dwithin(geog, '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography, 10000::double precision, true))
Rows Removed by Filter: 37041
Total runtime: 579.934 ms

半径1km(2km四方)を範囲とした場合は次の通り。

explain analyze SELECT
pref_name,city_name,street_name,address,
ST_Distance(geog, ST_GeogFromText('POINT(139.642403 35.448575)'))
FROM address WHERE
ST_DWithIn(geog, ST_GeogFromText('POINT(139.642403 35.448575)'), 1000)
ORDER BY 5 LIMIT 1
;

• • • • • • • -

Limit (cost=11.80..11.81 rows=1 width=269) (actual time=19.153..19.153 rows=1 loops=1)
-> Sort (cost=11.80..11.81 rows=1 width=269) (actual time=19.150..19.150 rows=1 loops=1)
Sort Key: (_st_distance(geog, '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography, 0::double precision, true))
Sort Method: top-N heapsort Memory: 17kB
-> Index Scan using ix_address_geog on address (cost=0.00..11.79 rows=1 width=269) (actual time=0.870..18.580 rows=815 loops=1)
Index Cond: (geog && '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography)
Filter: *3 AND _st_dwithin(geog, '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography, 1000::double precision, true))
Rows Removed by Filter: 672
Total runtime: 19.203 ms

もうちょっとがんばって半径500m(1km四方)。

explain analyze SELECT
pref_name,city_name,street_name,address,
ST_Distance(geog, ST_GeogFromText('POINT(139.642403 35.448575)'))
FROM address WHERE
ST_DWithIn(geog, ST_GeogFromText('POINT(139.642403 35.448575)'), 500)
ORDER BY 5 LIMIT 1
;

• • • • • • • -

Limit (cost=11.80..11.81 rows=1 width=269) (actual time=3.206..3.206 rows=1 loops=1)
-> Sort (cost=11.80..11.81 rows=1 width=269) (actual time=3.204..3.204 rows=1 loops=1)
Sort Key: (_st_distance(geog, '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography, 0::double precision, true))
Sort Method: top-N heapsort Memory: 17kB
-> Index Scan using ix_address_geog on address (cost=0.00..11.79 rows=1 width=269) (actual time=0.474..3.108 rows=262 loops=1)
Index Cond: (geog && '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography)
Filter: *4 AND _st_dwithin(geog, '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography, 500::double precision, true))
Rows Removed by Filter: 150
Total runtime: 3.232 ms

みるみる小さくなっていきます。

まとめ

• PostGIS関数の計算はできるだけしないように絞り込むのが吉。
• GEOGRAPHYは、経度緯度ベースなのにメートル単位の計測関数が使えて楽。

ここでは余談となりますが、ある程度GISとかに慣れてくるとGEOGRAPHYは計算コストがやはりかかるので使わない方向にむいていくことになると思います。

*1:geom && '010300000001000000050000004A0A2C802954614000000000008041404A0A2C8029546140F6B52E3542E541400000000000806140F6B52E3542E54140000000000080614000000000008041404A0A2C80295461400000000000804140'::geometry) AND (st_distance_spheroid(geom, '0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geometry, 'SPHEROID("GRS_1980",6378137,298.257222101004)'::spheroid) < 10000::double precision

*2:'0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography && _st_expand(geog, 10000::double precision

*3:'0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography && _st_expand(geog, 1000::double precision

*4:'0101000020E6100000467BBC908E746140CF66D5E76AB94140'::geography && _st_expand(geog, 500::double precision