[postgis-users] Help please with PL/R

[Gery .]

Hello,

Since a while I've been searching a way of smoothing long/lat data coming from ship tracks using Loess from PL/R. I found David Bitner's method (link below), here the links:


post: http://postgis.refractions.net/pipermail/postgis-users/2010-October/030264.html

image output: http://postgis.refractions.net/pipermail/postgis-users/attachments/20101015/b8efb3e7/attachment.png


I've been trying to reproduce what he did but had have no success in getting work with my data. My table of points looks like:

 id   | fldr |    sx    |    sy    |      profile       
-------+------+----------+----------+--------------------
     1 |    1 | -70.7033 | -26.5362 | LF
     2 |    2 | -70.7033 | -26.5362 | LF
     3 |    3 | -70.7033 | -26.5362 | LF
     4 |    4 | -70.7033 | -26.5362 | LF
     5 |    5 | -70.7033 | -26.5362 | LF
     6 |    6 | -70.7033 | -26.5362 | LF
     7 |    7 | -70.7033 | -26.5362 | LF
     8 |    8 | -70.7033 | -26.5362 | LF
     9 |    9 | -70.7033 | -26.5362 | LF
    10 |   10 | -70.7033 | -26.5362 | LF
    11 |   11 | -70.7033 | -26.5362 | LF
    12 |   12 | -70.7033 | -26.5362 | LF
    13 |   13 | -70.7033 | -26.5362 | LF
    14 |   14 | -70.7033 | -26.5362 | LF
    15 |   15 | -70.7033 | -26.5362 | LF
    16 |   16 | -70.7033 | -26.5362 | LF
    17 |   17 | -70.7033 | -26.5362 | LF
    18 |   18 | -70.7033 | -26.5362 | LF
    19 |   19 | -70.7033 | -26.5362 | LF
    20 |   20 | -70.7033 | -26.5362 | LF
    21 |   21 | -70.7033 | -26.5362 | LF
    22 |   22 | -70.7033 | -26.5362 | LF
    23 |   23 | -70.7033 | -26.5361 | LF
    24 |   24 | -70.7033 | -26.5361 | LF
    25 |   25 | -70.7033 | -26.5361 | LF
    26 |   26 | -70.7033 | -26.5361 | LF
    27 |   27 | -70.7032 | -26.5361 | LF
    28 |   28 | -70.7032 | -26.5361 | LF
    29 |   29 | -70.7032 | -26.5361 | LF
    30 |   30 | -70.7032 | -26.5361 | LF
    31 |   31 | -70.7032 | -26.5361 | LF
    32 |   32 | -70.7032 | -26.5361 | LF
    33 |   33 | -70.7032 | -26.5361 | LF
    34 |   34 | -70.7032 | -26.5361 | LF
    35 |   35 | -70.7032 | -26.5361 | LF
    36 |   36 | -70.7032 | -26.5361 | LF
    37 |   37 | -70.7032 | -26.5361 | LF
    38 |   38 | -70.7032 | -26.5361 | LF
    39 |   39 | -70.7032 | -26.5361 | LF
    40 |   40 | -70.7032 | -26.5361 | LF
...(continue up to 245 rows)

as you can see, the first 22 rows in sx/sy (ie. long/lat) are the same, something similar happens up to the end of the table (245 records), I attached an image of this track (yellow) with the smoothed version (red) using Bezier interpolation from here (http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=smooth_line_%28data_management%29).

I've been trying the code below with the following command, using my table of points:

mydb=# select astext(loess(geom)) from rawnav_point_wgs84;

however, it gives this error message:
CREATE FUNCTION
CREATE FUNCTION
psql:smoothingtest.sql:51: ERROR:  function loess(double precision[]) does not exist
LINE 8:    SELECT loess(array_agg(st_x(point))) AS x, loess(array_ag...
                  ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
psql:smoothingtest.sql:53: ERROR:  function loess(geometry) does not exist
LINE 1: select astext(loess(geom)) from chanaral_rawnav_point_wgs84;
                      ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

I think the code below works with a table of points as input, however I tried the table of lines made from these points but didn't work either.

I'd appreciate any support, thanks.

Please find the 'smoothingtest.sql' (modified version of David Bitner's method) here below:

*******************************************
CREATE OR REPLACE FUNCTION loess(double precision [], double precision[])
RETURNS double precision[] AS
$BODY$
        x <- arg1
        y <- arg2
        x.loess <- loess(x ~ y, span=0.75 data.frame(x=x, y=y))
        x.predict <- predict(x.loess, data.frame(y=y))
$BODY$
LANGUAGE plr VOLATILE STRICT
COST 100;

CREATE OR REPLACE FUNCTION dump_linestring_points(geometry)
RETURNS SETOF geometry AS
$BODY$
        SELECT (
                'SRID=' ||
                bar ||
                ';POINT(' ||
                foo ||
                ')'
                )::geometry
        FROM
                regexp_split_to_table(
                rtrim(
                substring(asewkt($1) from 23)
                ,')'
                ),
                ','
                ) as foo , srid($1) as bar
        WHERE st_geometrytype($1)='ST_LineString';
$BODY$
LANGUAGE sql IMMUTABLE STRICT
COST 1
ROWS 100;

CREATE OR REPLACE FUNCTION loess(geometry, double precision)
RETURNS geometry AS
$BODY$
        SELECT makeline(setsrid(makepoint(x,y),srid($1)))
        FROM (
                SELECT unnest(x) AS x, unnest(y) AS y
                FROM (
                        SELECT loess(array_agg(st_x(point))) AS x, loess(array_agg(st_y(point))) AS y
                        FROM (
                                SELECT dump_linestring_points($1) AS point
                        ) AS foo
                ) AS bar
        ) AS baz
$BODY$
LANGUAGE sql VOLATILE
COST 100;
*******************************************
 		 	   		  
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screenshot.png
Type: image/png
Size: 861697 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20120409/0f46a29d/attachment.png>