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;
