As an option, simultaneously with approximation, topographic parameters slope, aspect, profile curvature (measured in the direction of the steepest slope), tangential curvature (measured in the direction of a tangent to contour line) or mean curvature are computed and saved as raster maps specified by the options slope, aspect, pcurv, tcurv, mcurv respectively. If -d flag is set, v.surf.rst outputs partial derivatives fx,fy,fxx, fyy,fxy instead of slope, aspect, profile, tangential and mean curvatures respectively. If the input vector map have time stamp, the program creates time stamp for all output maps.
User can either use r.mask to set a mask or specify a raster map in mask option, which will be used as a mask. The approximation is skipped for cells which have zero or NULL value in mask. NULL values will be assigned to these cells in all output raster maps. Data points are checked for identical points and points that are closer to each other than the given dmin are removed. If sparsely digitized contours or isolines are used as input, additional points are computed between each 2 points on a line if the distance between them is greater than specified dmax. Parameter zmult allows user to rescale the values used for approximation (useful e.g. for transformation of elevations given in feet to meters, so that the proper values of slopes and curvatures can be computed).
Regularized spline with tension is used for the approximation. The tension parameter tunes the character of the resulting surface from thin plate to membrane. Smoothing parameter smooth controls the deviation between the given points and the resulting surface and it can be very effective in smoothing noisy data while preserving the geometrical properties of the surface. With the smoothing parameter set to zero (smooth=0) the resulting surface passes exactly through the data points (spatial interpolation is performed). When smoothing parameter is used, it is also possible to output a vector point map deviations containing deviations of the resulting surface from the given data.
If the number of given points is greater than segmax, segmented processing is used. The region is split into quadtree-based rectangular segments, each having less than segmax points and approximation is performed on each segment of the region. To ensure smooth connection of segments the approximation function for each segment is computed using the points in the given segment and the points in its neighborhood which are in the rectangular window surrounding the given segment. The number of points taken for approximation is controlled by npmin, the value of which must be larger than segmax. User can choose to output vector maps treeseg and overwin which represent the quad tree used for segmentation and overlapping neighborhoods from which additional points for approximation on each segment were taken.
Predictive error of surface approximation for given parameters can be computed using the -c flag. A crossvalidation procedure is then performed using the data given in the vector map input and the estimated predictive errors are stored in the vector point map cvdev. When using this flag, no raster output maps are computed. Anisotropic surfaces can be interpolated by setting anisotropy angle theta and scaling factor scalex. The program writes values of selected input and internally computed parameters to the history file of raster map elevation.
The user must run g.region before the program to set the region and resolution for approximation.
Topographic parameters are computed directly from the approximation function so that the important relationships between these parameters are preserved. The equations for computation of these parameters and their interpretation is described in Mitasova and Hofierka, 1993 or Neteler and Mitasova, 2004). Slopes and aspect are computed in degrees (0-90 and 1-360 respectively). The aspect raster map has value 0 assigned to flat areas (with slope less than 0.1%) and to singular points with undefined aspect. Aspect points downslope and is 90 to the North, 180 to the West, 270 to the South and 360 to the East, the values increase counterclockwise. Curvatures are positive for convex and negative for concave areas. Singular points with undefined curvatures have assigned zero values.
Tension and smoothing allow user to tune the surface character. For most landscape scale applications the default values should provide adequate results. The program gives warning when significant overshoots appear in the resulting surface and higher tension or smoothing should be used.
To select parameters that will produce a surface with desired properties, it is useful to know that the method is scale dependent and the tension works as a rescaling parameter (high tension "increases the distances between the points" and reduces the range of impact of each point, low tension "decreases the distance" and the points influence each other over longer range). Surface with tension set too high behaves like a membrane (rubber sheet stretched over the data points) with peak or pit ("crater") in each given point and everywhere else the surface goes rapidly to trend. If digitized contours are used as input data, high tension can cause artificial waves along contours. Lower tension and higher smoothing is suggested for such a case.
Surface with tension set too low behaves like a stiff steel plate and overshoots can appear in areas with rapid change of gradient and segmentation can be visible. Increase in tension should solve the problems.
There are two options how tension can be applied in relation to dnorm (dnorm rescales the coordinates depending on the average data density so that the size of segments with segmax=40 points is around 1 - this ensures the numerical stability of the computation):
The default is a recommended choice, however for the applications where the user needs to change density of data and preserve the approximation character the -t flag can be helpful.
Anisotropic data (e.g. geologic phenomena) can be interpolated
using theta and scalex defining orientation and ratio of
the perpendicular axes put on the longest/shortest side of the
feature, respectively. Theta is measured in degrees from East,
counterclockwise. Scalex is a ratio of axes sizes.
Setting scalex in the range 0-1, results in a pattern prolonged
in the direction defined by theta. Scalex value 0.5
means that modeled feature is approximately 2 times longer in the
direction of theta than in the perpendicular direction.
Scalex value 2 means that axes ratio is reverse resulting in a
pattern perpendicular to the previous example. Please note that
anisotropy option has not been extensively tested and may include bugs
(for example, topographic parameters may not be computed correctly) -
if there are problems, please report to GRASS bugtracker (accessible
For data with values changing over several magnitudes (sometimes the concentration or density data) it is suggested to interpolate the log of the values rather than the original ones.
v.surf.rst checks the numerical stability of the algorithm by computing the values in given points, and prints the root mean square deviation (rms) found into the history file of raster map elevation. For computation with smoothing set to 0, rms should be 0. Significant increase in tension is suggested if the rms is unexpectedly high for this case. With smoothing parameter greater than zero the surface will not pass exactly through the data points and the higher the parameter the closer the surface will be to the trend. The rms then represents a measure of smoothing effect on data. More detailed analysis of smoothing effects can be performed using the output deviations option.
v.surf.rst also writes the values of parameters used in computation into the comment part of history file elevation as well as the following values which help to evaluate the results and choose the suitable parameters: minimum and maximum z values in the data file (zmin_data, zmax_data) and in the interpolated raster map (zmin_int, zmax_int), rescaling parameter used for normalization (dnorm), which influences the tension.
If visible connection of segments appears, the program should be rerun with higher npmin to get more points from the neighborhood of given segment and/or with higher tension.
When the number of points in a vector map is not too large (less than 800), the user can skip segmentation by setting segmax to the number of data points or segmax=700.
v.surf.rst gives warning when user wants to interpolate outside the rectangle given by minimum and maximum coordinates in the vector map, zoom into the area where the given data are is suggested in this case.
When a mask is used, the program takes all points in the given region for approximation, including those in the area which is masked out, to ensure proper approximation along the border of the mask. It therefore does not mask out the data points, if this is desirable, it must be done outside v.surf.rst.
The "optimal" approximation parameters for given data can be found using a cross-validation (CV) procedure (-c flag). The CV procedure is based on removing one input data point at a time, performing the approximation for the location of the removed point using the remaining data points and calculating the difference between the actual and approximated value for the removed data point. The procedure is repeated until every data point has been, in turn, removed. This form of CV is also known as the "leave-one-out" or "jack-knife" method (Hofierka et al., 2002; Hofierka, 2005). The differences (residuals) are then stored in the cvdev output vector map. Please note that during the CV procedure no other output maps can be set, the approximation is performed only for locations defined by input data. To find "optimal parameters", the CV procedure must be iteratively performed for all reasonable combinations of the approximation parameters with small incremental steps (e.g. tension, smoothing) in order to find a combination with minimal statistical error (also called predictive error) defined by root mean squared error (RMSE), mean absolute error (MAE) or other error characteristics. A script with loops for tested RST parameters can do the job, necessary statistics can be calculated using e.g. v.univar. It should be noted that crossvalidation is a time-consuming procedure, usually reasonable for up to several thousands of points. For larger data sets, CV should be applied to a representative subset of the data. The cross-validation procedure works well only for well-sampled phenomena and when minimizing the predictive error is the goal. The parameters found by minimizing the predictive (CV) error may not not be the best for for poorly sampled phenomena (result could be strongly smoothed with lost details and fluctuations) or when significant noise is present that needs to be smoothed out.
North Carolina example (we simulate randomly distributed elevation measures which we interpolate to a gap-free elevation surface):
g.region raster=elevation -p # random elevation extraction of 500 samplings r.random elevation vector_output=elevrand n=500 v.info -c elevrand v.db.select elevrand # interpolation based on all points v.surf.rst elevrand zcol=value elevation=elev_full # apply the color table of the original raster map r.colors elev_full raster=elevation d.rast elev_full d.vect elevrand # interpolation based on subset of points (only those over 1300m/asl) v.surf.rst elevrand zcol=value elevation=elev_partial where="value > 1300" r.colors elev_partial raster=elevation d.rast elev_partial d.vect elevrand where="value > 1300"
Overview: Interpolation and Resampling in GRASS GIS
Modified program (translated to C, adapted for GRASS, new segmentation
Irina Kosinovsky, US Army CERL, Dave Gerdes, US Army CERL
Modifications for new sites format and timestamping:
Darrel McCauley, Purdue University, Bill Brown, US Army CERL
Update for GRASS5.7, GRASS6 and addition of crossvalidation:
Jaroslav Hofierka, University of Presov; Radim Blazek, ITC-irst
Parallelization using OpenMP:
Stanislav Zubal, Czech Technical University in Prague
Michal Lacko, Pavol Jozef Safarik University in Kosice
Last changed: $Date$
© 2003-2020 GRASS Development Team, GRASS GIS 7.6.2dev Reference Manual