[postgis-users] Finding closet intersect point along direction of line

[adrian.kitchingman at dse.vic.gov.au]

Thanks for the code Simon. After a few tweaks I got it working perfectly 
in PostGIS2.

Here's the updated code.

Adrian

CREATE OR REPLACE FUNCTION ST_Extend(eje_ geometry, bound_ geometry)
RETURNS geometry
AS $$
-- Use Case: I need to "extend" a linestring in another one. The new 
linestring must be
-- obtained such as their new extremes are the intersection of it with a
-- polygon that contains it. The original linestring is composed of only a
-- linear segment.
-- version: alfa , by Julio A. Galindo, April 17, 2007: juliogalindoq at 
gmail.com
DECLARE
      b_ geometry = st_boundary(bound_);
      dist float;
      max_dist float = 0;
      n_points int;
      pto_1 geometry;
      pto_2 geometry;
      first_pto geometry;
      last_pto geometry;
      u_1 float;
      u_2 float;
      norm float;
      result text = 'LINESTRING(';
BEGIN
      IF  st_GeometryType(eje_)   NOT LIKE 'ST_LineString' OR
           st_GeometryType(bound_) NOT LIKE 'ST_Polygon' THEN
         RETURN NULL;
      END IF;
      -- First Search how far is the boundary: (worst case)
      pto_1 := st_StartPoint(eje_);
      pto_2 := st_EndPoint(eje_);
      FOR i IN 1..st_NumPoints(b_)-1 LOOP
         dist := st_distance(st_PointN(b_,i),pto_1);
         IF dist > max_dist THEN max_dist := dist; --END IF;
            max_dist := dist;
         END IF;
         dist := st_distance(st_PointN(b_,i),pto_2);
         IF dist > max_dist THEN max_dist := dist; --END IF;
            max_dist := dist;
         END IF;
      END LOOP;
      -- Now extent the linestring:
      pto_2 := st_PointN(eje_,2);
      u_1 := st_X(pto_2)-st_X(pto_1);
      u_2 := st_Y(pto_2)-st_Y(pto_1);
      norm := sqrt(u_1^2 + u_2^2);
      first_pto := 
st_MakePoint(st_X(pto_1)-u_1/norm*dist,st_Y(pto_1)-u_2/norm*dist);
      n_points := st_nPoints(eje_);
      IF n_points > 2 THEN
         pto_1 := st_PointN(eje_,n_points-1);
         pto_2 := st_PointN(eje_,n_points);
         u_1   := st_X(pto_2)-st_X(pto_1);
         u_2   := st_Y(pto_2)-st_Y(pto_1);
         norm  := sqrt(u_1^2 + u_2^2);
      END IF;
      last_pto := 
st_MakePoint(st_X(pto_2)+u_1/norm*dist,st_Y(pto_2)+u_2/norm*dist);
      result   := result || st_X(first_pto) || ' ' || st_Y(first_pto) || 
',';
      FOR i IN 1..st_NumPoints(eje_) LOOP
          result := result || st_X(st_PointN(eje_,i)) || ' ' || 
st_Y(st_PointN(eje_,i)) || ',';
      END LOOP;
      result := result || st_X(last_pto) || ' ' || st_Y(last_pto) || ')';
      -- Find the final Linestring:
      b_ := st_intersection(st_GeomFromText(result,st_SRID(eje_)),bound_);
      RETURN b_;
  END $$
  LANGUAGE plpgsql
  STABLE
  RETURNS NULL ON NULL INPUT;





Notice:
This email and any attachments may contain information that is personal, confidential,
legally privileged and/or copyright. No part of it should be reproduced, adapted or communicated without the prior written consent of the copyright owner. 

It is the responsibility of the recipient to check for and remove viruses.

If you have received this email in error, please notify the sender by return email, delete it from your system and destroy any copies. You are not authorised to use, communicate or rely on the information contained in this email.

Please consider the environment before printing this email.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20130411/7b9313c8/attachment.html>