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

[Stephen Woodbridge]

Hi Strk,

This is great news, you guys really work fast :) Well I guess my 
function will help people that aren't yet on HEAD CVS. ;)

-Steve


strk at refractions.net wrote:
> 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
> 
> 
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>