postgisを使ってみる

postgis はインストール済み前提です。

テーブル作成

CREATE TABLE geodata
(
id integer NOT NULL,
"name" character varying(256) NOT NULL,
geo geometry NOT NULL,
CONSTRAINT geodata_pkey PRIMARY KEY (id)
)

データ投入

INSERT INTO geodata (id, name, geo) VALUES
(1,'上野駅', GeomFromText('POINT(139.777254 35.713768)', 4326)),
(2,'西郷隆盛像', GeomFromText('POINT(139.774029 35.711846)', 4326)),
(3,'上野の森美術館', GeomFromText('POINT(139.774744 35.712737)', 4326)),
(4,'不忍池弁財天', GeomFromText('POINT(139.770872 35.712351)', 4326)),
(5,'野口英世博士像', GeomFromText('POINT(139.775696 35.716293)', 4326)),
(6,'国立西洋美術館', GeomFromText('POINT(139.775803 35.71542)', 4326)),
(7,'国立科学博物館', GeomFromText('POINT(139.776544 35.716319)', 4326)),
(8,'東京都美術館', GeomFromText('POINT(139.772776 35.717186)', 4326)),
(9,'東京国立博物館', GeomFromText('POINT(139.776462 35.718883)', 4326)),
(10,'花やしき', GeomFromText('POINT(139.794547 35.71528)', 4326)),
(11,'雷門', GeomFromText('POINT(139.792692 35.710635)', 4326));

function作成

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;

検索してみる

SELECT * FROM geodata WHERE expand_sphere_pseudo(GeomFromText('POINT(139.774744 35.712737)', 4326),1000) && geo AND distance_sphere(geo, GeomFromText('POINT(139.774744 35.712737)', 4326)) < 1000;

おぉ！検索できますね。
多少の誤差はあるでしょうが、とりあえずいけてるっぽい！

ソートする場合

SELECT id, name, distance_sphere(geo, GeometryFromText('POINT(139.774744 35.712737)', 4326)) as distance FROM geodata WHERE expand_sphere_pseudo(GeomFromText('POINT(139.774744 35.712737)', 4326), 1000) && geo AND distance_sphere(geo, GeometryFromText('POINT(139.774744 35.712737)', 4326)) < 1000 ORDER BY distance;

これで1000m以内の近い順。

geometryのインデックス

CREATE INDEX geodata_geo_idx
ON geodata
USING gist
(geo);

こんな感じでgistというインデックスがはれるらしい。
本当にこれで早くなるのか・・・
テストデータ作るのが面倒なので、とりあえずここまで。