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

[Simon Greener]

Adrian,

Below is a function that I found on the web that should do what you want.

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 = 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 ( GeometryType(eje_)   NOT LIKE 'LINESTRING' OR
           GeometryType(bound_) NOT LIKE 'POLYGON' THEN
         RETURN NULL;
      END IF;
      -- First Search how far is the boundary: (worst case)
      pto_1 := StartPoint(eje_);
      pto_2 := EndPoint(eje_);
      FOR i IN 1..NumPoints(b_)-1 LOOP
         dist := distance(PointN(b_,i),pto_1);
         IF dist > max_dist THEN max_dist := dist; END IF;
            max_dist := dist;
         END IF;
         dist := distance(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 := PointN(eje_,2);
      u_1 := X(pto_2)-X(pto_1);
      u_2 := Y(pto_2)-Y(pto_1);
      norm := sqrt(u_1^2 + u_2^2);
      first_pto := MakePoint(X(pto_1)-u_1/norm*dist,Y(pto_1)-u_2/norm*dist);
      n_points := nPoints(eje_);
      IF n_points > 2 THEN
         pto_1 := PointN(eje_,n_points-1);
         pto_2 := PointN(eje_,n_points);
         u_1   := X(pto_2)-X(pto_1);
         u_2   := Y(pto_2)-Y(pto_1);
         norm  := sqrt(u_1^2 + u_2^2);
      END IF;
      last_pto := MakePoint(X(pto_2)+u_1/norm*dist,Y(pto_2)+u_2/norm*dist);
      result   := result || X(first_pto) || ' ' || Y(first_pto) || ',';
      FOR i IN 1..NumPoints(eje_) LOOP
          result := result || X(PointN(eje_,i)) || ' ' || Y(PointN(eje_,i)) || ',';
      END LOOP;
      result := result || X(last_pto) || ' ' || Y(last_pto) || ')';
      -- Find the final Linestring:
      b_ := intersection(GeomFromText(result,SRID(eje_)),bound_);
      RETURN b_;
  END $$
  LANGUAGE plpgsql
  STABLE
  RETURNS NULL ON NULL INPUT;

-- Test 1:
select asText(extendLine(geomFromText('LINESTRING(5 5, 12 12)',-1),GeomFromText('POLYGON((0 0,15 0,15 20,0 20,0 0))',-1)));

select asText(extendLine(geomFromText('LINESTRING(5 5, 12 12, 13 12)',-1),GeomFromText('POLYGON((0 0,15 0,15 20,0 20,0 0))',-1)));

regards
Simon


On Wed, 10 Apr 2013 15:34:04 +1000, <adrian.kitchingman at dse.vic.gov.au> wrote:

> Hi All
>
> I'm looking to see if anyone has a PostGIS
> solution to finding the closest point on another polyline to a line end
> point along the line's azimuth. The picture attached shows what I'm looking
> for. For end point A the usual closest point is B but I want C. The picture
> shows a polygon but I expect a conversion to polyline to get points on
> the edge.
>
> Cheers
>
> Adrian
>
> Using: PostGIS 2.0 on PostgreSQL 9.2,
> Windows 7
>
>
>
>
>
> Notice:
> This email and any attachments may contain information that ispersonal, confidential,
> legally privileged and/or copyright. No part of itshould be reproduced, adapted or communicated without the prior written consentof the copyright owner.
>
>
> It is the responsibility of the recipient to check for and removeviruses.
>
>
>
> If you have received this email in error, please notify the sender by returnemail, delete it from your system and destroy any copies. You are not authorisedto use, communicate or rely on the information contained in this email.
>
>
>
> Please consider the environment before printing this email.
>
>



-- 
Holder of "2011 Oracle Spatial Excellence Award for Education and Research."
SpatialDB Advice and Design, Solutions Architecture and Programming,
Oracle Database 10g Administrator Certified Associate; Oracle Database 10g SQL Certified Professional
Oracle Spatial, SQL Server, PostGIS, MySQL, ArcSDE, Manifold GIS, FME, Radius Topology and Studio Specialist.
39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia.
Website: www.spatialdbadvisor.com
   Email: simon at spatialdbadvisor.com
   Voice: +61 362 396397
Mobile: +61 418 396391
Skype: sggreener
Longitude: 147.20515 (147° 12' 18" E)
Latitude: -43.01530 (43° 00' 55" S)
GeoHash: r22em9r98wg
NAC:W80CK 7SWP3