[postgis-users] project a point onto a line and get x, y, offset of the point of projection

[Alex Smith]

Thanks for posting this.  Let me make sure that I have this correct - as we may need something similar to this in an upcoming project.  It sounds like this function essentially "snaps" a point to the nearest line.  Correct?  If so, does it snap the point to the nearest VERTEX on the line?  Or will it place the point at the appropriate location between two vertex's on the line?  I just ask because an Avenue script ArcView 3.2 would place the point at the nearest vertex and therefore not do what we need.  
 
Thanks again for sharing!
 
AS

Stephen Woodbridge <woodbri at swoodbridge.com> wrote:
Hi all,

Someone on the list asked if they could get the point location of where 
on a line object the shortest distance to that line was from a point. 
The distance(line, point) will tell you the shortest distance, but it 
does not project the point onto the line at the point where the shortest 
distance is located.

Then I needed the same thing for a project at Where2GetIt.com and it 
seemed like this functionality should be a part of PostGIS. So attached 
is the code I wrote for our needs. There are some examples of its use. 
It has an assumption that you are working in decimal degrees, but it 
could be modified easily to work in any space. I will leave that up to 
others.

The algorithm is very straight forward. Given a linestring and a point, 
project the point onto each segment of the linestring and save the 
projection of the closest point. Then compute the distance from the 
start of the linestring to the projection point and return an array of 
pointx, pointy, distance.

Paul, please feel free to add this or something like it to PostGIS.

Enjoy, I would be open to any feedback or improvements any might suggest.

-Steve
CREATE OR REPLACE FUNCTION x_y_offset_from_line_pnt (geometry, geometry) RETURNS float[] AS '
/*
* {x, y, offset} = x_y_offset_from_line_pnt(LINESTRING, POINT($lon $lat))
*/
DECLARE
line ALIAS FOR $1;
pnt ALIAS FOR $2;
nline geometry;
ipnt geometry;
bestp geometry;
pa geometry;
pb geometry;
npts integer;
besti integer;
seglen float;
dist float;
best float;
offset float;
U float;
--
BEGIN
IF NOT (GeometryType(line) = ''LINESTRING'' ||
GeometryType(line) = ''MULTILINESTRING'' ) AND
GeometryType(pnt) != ''POINT'' THEN
RETURN NULL;
END IF;
/*
* iterate over the linestring segments and find the segment and 
* point on the linestring that are closest to the pnt
*/
npts := npoints(line);
best := 360.0;
FOR i IN 1 .. npts-1 LOOP
pa := PointN(line,i);
pb := PointN(line,i+1);
seglen := distance(pa, pb);
IF (seglen > 0.0) THEN
U := ( (X(pnt)-X(pa)) * (X(pb)-X(pa)) +
(Y(pnt)-Y(pa)) * (Y(pb)-Y(pb)) ) /
(seglen*seglen);
IF U < 0.0 THEN
ipnt := pa;
ELSIF U > 1.0 THEN
ipnt := pb;
ELSE
ipnt := MakePoint(
X(pa) + U*(X(pb) - X(pa)),
Y(pa) + U*(Y(pb) - Y(pa)) );
END IF;
dist := distance(pnt, ipnt);
IF dist < best THEN
best := dist;
besti := i;
bestp := ipnt;
END IF;
END IF;
END LOOP;
/*
* we should now have the best point so calcuate the offset
* and return the results.
*/
IF besti = 1 THEN
nline = MakeLine(PointN(line,1), bestp);
ELSE
nline = MakeLine(PointN(line,1), PointN(line,2));
IF besti > 2 THEN
FOR i IN 3 .. besti LOOP
pa := PointN(line,i);
nline := AddPoint(nline, pa, -1);
END LOOP;
END IF;
nline := AddPoint(nline, bestp, -1);
END IF;
offset := length_spheroid(nline,''SPHEROID["GRS_1980", 6378137, 298.257222101]'');
RETURN ARRAY[X(bestp), Y(bestp), offset];
END;
' LANGUAGE plpgsql;

-- here are some test calls to the function

select x_y_offset_from_line_pnt(GeomFromText('LINESTRING(0 0,1 0)',-1), GeomFromText('POINT(0.5 0.1)', -1));

select x_y_offset_from_line_pnt(GeomFromText('LINESTRING(0 0,0.1 0,0.3 .1,.4 .15,.6 .2)',-1), GeomFromText('POINT(0.5 0.1)', -1));

select link_id, the_geom, min(distance(the_geom, GeomFromText('POINT(-71.38900904888 42.619402684661)', -1))) as dist
from 
rgeo 
where
expand(GeomFromText('POINT(-71.38900904888 42.619402684661)', -1), 0.0013) &&
the_geom
group by the_geom order by dist limit 1;

select (x_y_offset_from_line_pnt(
(select the_geom from rgeo where link_id=$segment_id), 
GeomFromText('POINT(-71.38900904888 42.619402684661)', -1))
)[3] as offset;
_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

		
---------------------------------
Yahoo! Sports
 Rekindle the Rivalries. Sign up for Fantasy Football
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20050703/6a86e8fa/attachment.html>