RST Spline Surfaces

RST: Regularized Spline with Tension

Overview

Thin plate splines fit a surface to a series of points by minimizing the energy it takes to "bend" the surface. You can think of it as stretching a rubber sheet between all the points and then varying the stiffness of the sheet.

Using this method you can create raster surfaces which, as compared to a simple rasterization method like a sampled TIN, better match the true topography, work better results from groundwater flow modules, and create much nicer visualization images.

A quadtree segmentation method is used to concentrate computational effort into areas of greatest input point density.

Modules

• v.surf.rst
• v.vol.rst
• r.resamp.rst
• librst
• r.fillnulls

Tuning parameters and Validation

• See the module help pages
  • [Video lecture on spline interpolation] by Prof Helena Mitasova

(TODO: write user guide)

• Cross-validation

(TODO: write How-to)

See manual page of v.surf.rst

See C-V script: http://www.grassbook.org/examples_menu3rd.php (Script to perform cross validation in GRASS/RST splines interpolation) for parametrisation

See this post to the mailing list

• To help speed things up during tuning tests use a smaller area, e.g. 100x100 cells, that is representative of the terrain.

Getting rid of those segmentation boxes

Here some work-around (avoiding segmentation at all):

Example:

      6921 input points
      ...
      WARNING: Points are more dense than specified 
      'DMIN'--ignored 5575 points (remain 1346)

Now we can calculate: 6921-5575=1346. For smooth connection of segments, npmin > segmax.

You may get errors like:

       ERROR: segmentation parameters set to invalid values: 
       npmin = 300, segmax = 850
       for smooth connection of segments, npmin > segmax (see manual)
       ERROR: input failed!

In the v.surf.rst/v.vol.rst source code, there is a MAXPOINTS value defined in the surf.h file. MAXPOINTS simply says how many input points are allowed to be computed without segmentation. I had to increase the value of MAXPOINTS in surf.h to 1346 (or a bit more), recompile, and then use segmax=1346 to avoid segmentation. Tricky but works. (n.b. actually this source code change probably has no effect other than to mask the warning message for v.surf.rst. It seems to be used in v.vol.rst though)

Update I have a patch which I currently test to replace MAXPOINTS with a new npmax parameter - Markus Neteler 12:56, 7 April 2009 (UTC)

The idea of the module authors is to change v.surf.rst/v.vol.rst to a more dynamic management of this value, the question is how to implement that.

Summary: avoid segmentation if possible.

Computational time needed

GRASS 7.svn, expirical formulae (to be stabilized or improved):

(t in minutes; points: number of vector points)

• v.surf.rst: t ~= 1.5e-11 * points^2
• v.vol.rst: t ~= 4.5e-11 * points^3

Working with contour lines

The RST modules work best with an even distribution of starting points. Contour line vertices tend to be highly packed along lines with huge gaps in between different lines. Not very evenly distributed! This will usually lead to segmentation problems (see troubleshooting section below). Some solutions are offered at the following link, demonstrating the use of the v.generalize module to remove vertices from the contour lines:

http://skagit.meas.ncsu.edu/~helena/grasswork/interpgen.html

You may prefer to try the r.surf.contour or add-on r.surf.nnbathy modules if starting with contour line data.

• See the Contour lines to DEM wiki page for a study of different interpolation methods.

Working with huge data sets

Interpolation large DEMs from large amounts of vector points requires lot of RAM. Some hints

• don't create point topology (-b flag in several modules like r.random, v.in.ascii etc)
• use segmentation (e.g., segmax=10 and npmin=120)
• interpolate in chunks with overlap

todo: explain more

Troubleshooting

• If your input data points include a dramatic change in point density across the region of interest you may get square artifacts in the output map. This is especially visible in a slope map or in contour lines created from the elev map with r.contour.
  • Solution 1: increase the npmin parameter
  • Solution 2: Create the surface in multiple passes. First isolate the low density points and use v.hull to make a polygon around them. Next run v.surf.rst for those points to make a smooth surface between them. Then you can either run v.to.rast on the hull area and create new random points with the r.random module, or sample the new surface with the v.random.cover add-on script. Finally use v.patch to create a new points file with a less dramatic change in point density.

• Spline surfaces do not handle hard breaks in the topography very well, such as found on a beach at the base of a cliff. "Ringing_artifacts" may be introduced and the module may complain about "overshoots", as the surface is bent away from an input point further than the module is happy with. Users familiar with digital electronics may recall the similar problem of fitting a sinusoid to a square wave.

ToDo

• add parallelization support to the RST library (e.g. OpenMP)