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

strk at refractions.net strk at refractions.net
Sun Jul 3 10:45:03 PDT 2005

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

```