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

[Simon Greener]

Adrian,

That's excellent.

regards
Simon

On Thu, 11 Apr 2013 16:51:14 +1000, <adrian.kitchingman at dse.vic.gov.au> wrote:

> 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.



-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20130411/e6789359/attachment.html>