[postgis-users] Loess Data Smoothing With R

[David William Bitner]

Attached is an example of using Loess prediction in PL/R for smoothing xyzm
data (flight tracks) with respect to the m (time) value. Compared to
previous work that I had posted that used a time-weighted moving average
(similar to McMaster method) with the windowing functions, using the Loess
prediction is far less influenced by outliers. In the attached image the red
track is the raw track (data from multilateration sensors). You can see that
there is a lot of noise to this data. The blue line is the result of the
following functions to smooth the data.

This is the workhorse PL/R function. It takes two arrays and a span value.
The first array is the data you want to smooth (in the case of the 4d data
this is an array of the x, y, or z values), the second array is the array of
time values that are the steps used for the smoothing. The final value is
the span value that controls what portion of the data is considered for each
prediction.

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



This function wraps the PL/R call and wraps it to call the loess function on
each dimension. This function unwraps the data to and from the single
dimension arrays to pass into the above PL/R function. The span value that
this wrapper takes is an absolute unit in the units of the m value that get
converted into the proportion that is used in the above function.

CREATE OR REPLACE FUNCTION loess(geometry, double precision)
>   RETURNS geometry AS
> $BODY$
> SELECT
> makeline(setsrid(makepoint(x,y,z,m),srid($1))) from (
>  select unnest(x) as x,unnest(y) as y,unnest(z) as z,unnest(m) as m from (
> select
> loess(array_agg(st_x(point)),array_agg(st_m(point)),$2/max(lastm)) as x,
>  loess(array_agg(st_y(point)),array_agg(st_m(point)),$2/max(lastm)) as y,
> loess(array_agg(st_z(point)),array_agg(st_m(point)),$2/max(lastm)) as z,
>  array_agg(st_m(point)) as m
> from (
> select dump_linestring_points($1) as point ,st_m(st_endpoint($1)) as lastm
>  ) as foo
> ) as bar ) as baz
> $BODY$
>   LANGUAGE sql VOLATILE
>   COST 100;


dump_linestring_points is a hack to get around the performance limitations
of st_dumppoints that simply parses the EWKT to create the list of points
rather than creating a loop through the points.

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;



Hopefully this can be of use to someone else who has noisy track type data!

David

-- 
************************************
David William Bitner
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20101015/b8efb3e7/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smooth.png
Type: image/png
Size: 18201 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20101015/b8efb3e7/attachment.png>