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

[Stephen Woodbridge]

Alex Smith wrote:
> 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? 

given a point(P) and a line(S-E) it projects the point onto the line 
where the snapped point(X) is at the shortest distance from the point to 
the line. The shortest distance it always the one closest to the 
beginning of the line if there are multiple points all with the same 
shortest distance. The point may be a vertex, but it may be along the 
segment if that is the shortest distance. The offset is the distance 
along the line from the start to the snap point.

     *----X----*
    /           \
   /      P      \
  /               E
S

If S-E is the linestring, P is the point, then X is the snap point and 
the function returns ARRAY[Xx, Xy, offset] where POINT(Xx, Xy) is the 
snap point and offset = length of the segment S-X

HTH,
   -Steve

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