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

[strk at refractions.net]

In the head branch of postgis we already have functions able to
do this:

	-- position is a float between 0 and 1
	line_interpolate_point(line, position)

	-- return a float between 0 and 1 representing
	-- position of closest point on line
	line_locate_point(line, point)

So functionality of your function would be obtained with:

 SELECT line_interpolate_point(geom, line_locate_point(geom, 'POINT..'));

--strk;
	
On Fri, Jul 01, 2005 at 10:58:36PM -0400, Stephen Woodbridge 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