RST Spline Surfaces: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
m (→‎Overview: quadtree)
(+https://www.youtube.com/watch?v=eLYsx-ADOGI spline interpolation lecture)
 
(12 intermediate revisions by 2 users not shown)
Line 16: Line 16:
* {{cmd|r.fillnulls}}
* {{cmd|r.fillnulls}}


=== Tuning parameters ===
=== Tuning parameters and Validation ===


* See the module help pages
* See the module help pages
** [[https://www.youtube.com/watch?v=eLYsx-ADOGI Video lecture on spline interpolation]] by Prof Helena Mitasova
: '''''(TODO: write user guide)'''''
: '''''(TODO: write user guide)'''''


=== Validation ===
* Cross-validation
* Cross-validation
: '''''(TODO: write How-to)'''''
: '''''(TODO: write How-to)'''''


:: See manual page of v.surf.rst
:: See manual page of {{cmd|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 C-V script: http://www.grassbook.org/examples_menu3rd.php (Script to perform cross validation in GRASS/RST splines interpolation) for parametrisation
:: See this [http://lists.osgeo.org/pipermail/grass-user/2007-May/039530.html 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 ===
=== Getting rid of those segmentation boxes ===
Line 47: Line 49:
         ERROR: input failed!
         ERROR: input failed!


In the v.surf.rst/v.vol.rst source code, there is a MAXPOINTS value defined in the <tt>surf.h</tt> file. MAXPOINTS simply says how many input points are allowed to be computed without segmentation. I had to increase the value of MAXPOINTS in <tt>surf.h</tt> to 1346 (or a bit more), recompile, and then use segmax=1346 to avoid segmentation.  
In the {{cmd|v.surf.rst}}/{{cmd|v.vol.rst}} source code, there is a MAXPOINTS value defined in the <tt>surf.h</tt> file. MAXPOINTS simply says how many input points are allowed to be computed without segmentation. I had to increase the value of MAXPOINTS in <tt>surf.h</tt> 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'')
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 [[User:Neteler|Neteler]] 12:56, 7 April 2009 (UTC)
: Update I have a patch which I currently test to replace MAXPOINTS with a new npmax parameter - Markus [[User:Neteler|Neteler]] 12:56, 7 April 2009 (UTC)
Line 58: Line 60:


GRASS 7.svn, expirical formulae (to be stabilized or improved):
GRASS 7.svn, expirical formulae (to be stabilized or improved):
* v.surf.rst: t ~ 1.5e-11 * 7838946^2 (t in minutes; points: number of vector points)
: ''('''t''' in minutes; '''points''': number of vector points)''
* v.vol.rst: t ~ 4.5e-11 * points^3   (t in minutes; points: number of vector points)
* {{cmd|v.surf.rst}}: t ~= 1.5e-11 * points^2
* {{cmd|v.vol.rst}}: t ~= 4.5e-11 * points^3


=== Working with contour lines ===
=== 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:
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 {{cmd|v.generalize}} module to remove vertices from the contour lines:
: http://skagit.meas.ncsu.edu/~helena/grasswork/interpgen.html
: 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.
You may prefer to try the {{cmd|r.surf.contour}} or add-on {{AddonCmd|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 ===
=== Working with huge data sets ===
Line 72: Line 77:
Interpolation large DEMs from large amounts of vector points requires lot of RAM. Some hints
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)
* don't create point topology (-b flag in several modules like {{cmd|r.random}}, {{cmd|v.in.ascii}} etc)
* use segmentation (e.g., segmax=10 and npmin=120)
* use segmentation (e.g., segmax=10 and npmin=120)
* interpolate in chunks with overlap
* interpolate in chunks with overlap
Line 80: Line 85:
=== Troubleshooting ===
=== 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.
* 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 {{cmd|r.contour}}.
** Solution 1: increase the ''npmin'' parameter
** 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 addon script. Finally use v.patch to create a new points file with a less dramatic change in point density.
** Solution 2: Create the surface in multiple passes. First isolate the low density points and use {{cmd|v.hull}} to make a polygon around them. Next run {{cmd|v.surf.rst}} for those points to make a smooth surface between them. Then you can either run {{cmd|v.to.rast}} on the hull area and create new random points with the {{cmd|r.random}} module, or sample the new surface with the {{AddonCmd|v.random.cover}} add-on script. Finally use {{cmd|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" 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.
* Spline surfaces do not handle hard breaks in the topography very well, such as found on a beach at the base of a cliff. "{{wikipedia|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 ===
=== ToDo ===
Line 92: Line 97:


[[Category:FAQ]]
[[Category:FAQ]]
[[Category:Interpolation]]
[[Category: Parallelization]]
[[Category: Parallelization]]

Latest revision as of 05:50, 30 June 2017

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

Tuning parameters and Validation

(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)

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.

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)