[OSGeoJapan-discuss] PostGIS用の便利なファンクション

[Hiroo Imaki]

リストの皆様、

先日、pgRoutingを使ったダムから河口までの距離計算について質問した今木です。しばらくこの問題と格闘している間に、いくつか便利なファンクションを書いたので（といっても誰かの書いたものを取り込んだり、改善したりしたものですが）皆さんのお役に立てばと思いポストします。

一つ目は、複数の点（POINT)を線（LINESTRING）にスナップするファンクションです。これは、Paul
Ramseyのブログからメインのコードを拝借してファンクションにしたものです。使い方は、
SELECT * FROM points_snap2lines
   ('SELECT gid AS id, the_geom AS geom FROM point_table',
    'SELECT gid AS id, the_geom AS geom FROM line_table',
     0.0001);
のようにして、最後に地図の単位でスナップの許容度を入力します。このファンクションを実行するとポイント、ラインのそれぞれのIDとスナップ後のジオメトリーが帰ってきます。

2つ目は、複数の点を使って複数の線をクリップするファンクションです。これは、
http://postgis.refractions.net/pipermail/postgis-users/2007-September/017159.htmlの歩スティングにあったファンクションを書き換えたものです。このポスティングのままでは私のpostgresql8.3.3ではうまく動かなかったので、手直ししたのと、ST_DWithin()を使って多少、点が線からずれていてもクリップできるようにしたのと、基本的なループ上のエラーを直したものです。使い方は、
SELECT * FROM split_lines2('SELECT gid AS id, the_geom AS geom FROM
line_table',
                           'SELECT the_geom AS geom FROM point_table',
0.0001);
のようにして、 ST_DWithinの距離パラメターを最後に入力します。

私の環境（WindowsXP,
PostgreSQL8.3.3)ではうまく動きました。改良、バグなどありましたら連絡ください。使った感想などもお聞かせください。

いまき

CREATE OR REPLACE FUNCTION points_snap2lines (in pt_q text, in ln_q text, in
torelance float4, out p_id int, out l_id int, out p_geom geometry)
RETURNS SETOF RECORD AS
$$
DECLARE
  snap_q text;
  pointrec record;
BEGIN
 EXECUTE 'CREATE TEMP TABLE line_tmp as '|| ln_q;
 EXECUTE 'CREATE TEMP TABLE point_tmp as '|| pt_q;
snap_q :='
SELECT
  pt_id,
  ln_id,
  ST_line_interpolate_point(
    ln_geom,
    ST_line_locate_point(ln_geom, pt_geom)
  ) AS the_geom
FROM
  (
    SELECT DISTINCT ON (pt.id)
      ln.the_geom AS ln_geom,
      pt.the_geom AS pt_geom,
      ln.id AS ln_id,
      pt.id AS pt_id
    FROM
      point_tmp pt INNER JOIN
      line_tmp ln
    ON
      ST_DWithin(pt.the_geom, ln.the_geom, '||cast(torelance as text) ||')
    ORDER BY
      pt.id,ST_Distance(ln.the_geom, pt.the_geom)
  ) as foo' ;
 FOR pointrec in EXECUTE snap_q LOOP
   p_id   := pointrec.pt_id;
   l_id   := pointrec.ln_id;
   p_geom := pointrec.the_geom;
   RETURN NEXT;
 END LOOP;
 DROP TABLE line_tmp;
 DROP TABLE point_tmp;
 RETURN;
END;
$$
LANGUAGE plpgsql;



CREATE OR REPLACE FUNCTION split_lines2(in lineq text, in pointq text, in
torelance float4, out lineid int, out the_geom geometry)
RETURNS SETOF RECORD AS
$$
DECLARE
  linerec record;
  pointrec record;
  linepos float;
  start_ float;
  end_ float;
  loopqry text;
BEGIN
  EXECUTE 'CREATE TEMP TABLE line_tmp as '|| lineq;
  EXECUTE 'CREATE TEMP TABLE point_tmp as '|| pointq;
  FOR linerec in EXECUTE 'SELECT * FROM line_tmp ORDER BY id' LOOP
  start_ := 0;
  end_   := 0;
  loopqry := '
    SELECT
        *, ST_line_locate_point('''||cast(linerec.geom as text)||''',geom)
AS frac
    FROM point_tmp
    WHERE ST_DWithin(geom,'''||cast(linerec.geom as text)||''',
'||torelance||')
    ORDER BY ST_line_locate_point('''||cast(linerec.geom as
text)||''',geom)';
      FOR pointrec in EXECUTE loopqry LOOP
          end_   := pointrec.frac;
          lineid := linerec.id;
          the_geom   := ST_line_substring(linerec.geom, start_, end_);
          start_ := end_;
        RETURN NEXT;
      END LOOP;
      lineid := linerec.id;
      the_geom:= ST_line_substring(linerec.geom, end_,1.0);
      RETURN NEXT;
  END LOOP;
  DROP TABLE line_tmp;
  DROP TABLE point_tmp;
  RETURN;
END;
$$
LANGUAGE plpgsql;



-- 
Hiroo Imaki
hiroo ＠ angeli.org
http://www.geopacific.org
-------------- next part --------------
HTML¤ÎźÉÕ¥Õ¥¡¥¤¥ë¤¬½üµî¤µ¤ì¤Þ¤·¤¿.
URL: http://lists.osgeo.org/pipermail/osgeojapan-discuss/attachments/20100721/f97c05af/attachment.html