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

[Stephen Woodbridge]

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