[postgis-users] Bezier or Spline smoothing implementation

[Michael.Scholz at dlr.de]

Hey!
First of all I can recommend two tutorials about how to use PL/R in PostGIS [1][2]. One possibility to interpolate splines with R is to use R's xspline package [3]. Some theoretical background 'bout those x-splines you'll find in the related paper [4]. I tried some other spline-methods in R but with dubious results. xspline works fine for me. A wrapping PL/R-function could look like this:

/**
* Interpolates a LineString using R's X-spline-method
*
 * Parameters:
*   1st - A *single* PostGIS-geometry (POINT, LINESTRING or POLYGON) *as text*: Use ST_AsEWKT() to pass a geometry as text.
*   2nd - Shape parameter used in R's X-spline-method [-1.0,...,1.0]. A shape of 0.0 does not interpolate anything. Default: 0.5.
*
 * Returns:
*   A PostGIS-LINESTRING-geometry covering the interpolated points.
*
 * Usage:
*   > SELECT interpolateXSpline(ST_AsEWKT(geom), 1.0) FROM your_table
*/
create or replace function interpolateXSpline(text, float DEFAULT 0.5) returns geometry as $$
  # Select the points (also from line or polygon-boundary) into a data frame
  points <- pg.spi.exec(
    sprintf("WITH pts AS (SELECT (ST_DumpPoints(ST_GeomFromEWKT(%1$s))).geom AS geom)
      SELECT ST_X(geom) AS x, ST_Y(geom) AS y FROM pts;",
        pg.quoteliteral(arg1)
    )
  )

  # Interpolate spline
  plot(points$x, points$y, type="n") # Dummy-call of plot() for xspline to work
  spline_pts <- xspline(points$x, points$y, shape=arg2, draw=FALSE)
  # Here you could simply output spline_pts and return.
  # Otherwise proceed to build a WKT linestring as shown below.

  # LineString-creation: Build beginning of WKT string for output to PostGIS
  out_line = "LINESTRING("
    for (i in 1:(length(spline_pts$x)-1)) {
      out_line = sprintf(
        "%1$s %2$f %3$f,",
        out_line, spline_pts$x[i], spline_pts$y[i]
      )
    }
  # LineString-creation: Append middle part
  out_line =
    sprintf("%1$s %2$f %3$f)",
      out_line,
      spline_pts$x[length(spline_pts$x)],
      spline_pts$y[length(spline_pts$y)]
    )
  # LineString-creation: Finalise (close)
  oline <- pg.spi.exec(
    sprintf("SELECT ST_GeomFromText('%1$s', ST_SRID(ST_GeomFromEWKT(%2$s))) AS ln;",
      out_line,
      pg.quoteliteral(arg1)
    )
 )
  return(oline$ln[1])
$$ language 'plr' IMMUTABLE;

I tried to get it to work with the 1st parameter as bytea to be able to pass binary geometry but got unserialize-errors in R. So for now you have to convert your input geometry to EWKT first and pass it as text.

Have fun! Michi

[1] http://www.bostongis.com/?content_name=postgresql_plr_tut01#87
[2] http://www.bostongis.com/?content_name=postgresql_plr_tut02#98
[3] http://stat.ethz.ch/R-manual/R-devel/library/graphics/html/xspline.html
[4] http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.44.5770&rep=rep1&type=pdf
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20130918/d3b5476e/attachment.html>