PostGIS の領域判定

　ここギコさんの記事からパクらせていただきました。
PostgreSQL には、GIST という領域型のインデックスがあり、これを利用できるとクエリを劇的に速くすることができます。さて、地図で利用する座標系は、いろいろあって、広い領域を扱おうと思うとJGD2000などの緯度経度の座標系を扱う事になると思います。ところが、点から半径１０ｋｍ圏内といった指定をするためには、距離を緯度経度に変換しなければならず、面倒です。しかも、長さは緯度により変化するので、厄介です。
　以下の関数は、指定図形を指定座標系(SRID)と、その単位により拡張した領域を返すものです。

CREATE OR REPLACE FUNCTION expand_pseudo(geometry, integer, double precision)
  RETURNS geometry AS
$BODY$
DECLARE
  sid integer; 
  geo geometry;
BEGIN
  -- SRIDを退避
  sid := SRID($1);
  -- 指定のSRIDにGeometryを変換する
  geo := Transform( $1, $2 );
  -- 指定SRIDの単位系で図形を拡張する
  geo := expand( geo, $3 );
  -- 元のSRIDへ変換する
  geo := Transform( geo, sid );
  -- バウンディング・ボックスを返す
  RETURN Envelope( geo );
END
$BODY$
  LANGUAGE 'plpgsql' IMMUTABLE
  COST 100;

　いちいち、geo という変数に代入していますが、連結した方が効率は良いかもしれません。何をやっているか、コメントを書いて、効率よりもわかり易さを重点に書いてます。

　使い方は、

DECLARE
  geo geometry;
  rec record;
BEGIN
  -- JGD2000: 4612
  -- 世界測地日本公共座標１２系: 2454
  geo := GeomFromString('POINT(140.223412,42.123456)', 4612);
  geo := expand_pseudo( geo, 2454, 10000 );
  FOR rec IN select * from hoge_table where geom && expand_pseudo( geo, 2454, 10000 ) LOOP
    --- 何か
  END LOOP;
END

という感じです。