v.surf.rst  Performs surface interpolation from vector points map
by splines.
Spatial approximation and topographic analysis from given point or isoline data
in vector format to floating point raster format using regularized spline with
tension.
vector, surface, interpolation, RST, 3D, nodata filling
v.surf.rst
v.surf.rst help
v.surf.rst [
ctd]
input=
name
[
layer=
string] [
zcolumn=
name]
[
where=
sql_query] [
elevation=
name]
[
slope=
name] [
aspect=
name]
[
pcurvature=
name] [
tcurvature=
name]
[
mcurvature=
name] [
deviations=
name]
[
cvdev=
name] [
treeseg=
name]
[
overwin=
name] [
nprocs=
integer]
[
mask=
name] [
tension=
float]
[
smooth=
float] [
smooth_column=
string]
[
segmax=
integer] [
npmin=
integer]
[
dmin=
float] [
dmax=
float]
[
zscale=
float] [
theta=
float]
[
scalex=
float] [
overwrite] [
help]
[
verbose] [
quiet] [
ui]
 c

Perform crossvalidation procedure without raster approximation
 t

Use scale dependent tension
 d

Output partial derivatives instead of topographic parameters
 overwrite

Allow output files to overwrite existing files
 help

Print usage summary
 verbose

Verbose module output
 quiet

Quiet module output
 ui

Force launching GUI dialog
 input=name [required]

Name of input vector map
Or data source for direct OGR access
 layer=string

Layer number or name
Vector features can have category values in different layers. This number
determines which layer to use. When used with direct OGR access this is
the layer name.
Default: 1
 zcolumn=name

Name of the attribute column with values to be used for approximation
If not given and input is 2D vector map then category values are used. If
input is 3D vector map then zcoordinates are used.
 where=sql_query

WHERE conditions of SQL statement without ’where’ keyword
Example: income < 1000 and population >= 10000
 elevation=name

Name for output surface elevation raster map
 slope=name

Name for output slope raster map
 aspect=name

Name for output aspect raster map
 pcurvature=name

Name for output profile curvature raster map
 tcurvature=name

Name for output tangential curvature raster map
 mcurvature=name

Name for output mean curvature raster map
 deviations=name

Name for output deviations vector point map
 cvdev=name

Name for output crossvalidation errors vector point map
 treeseg=name

Name for output vector map showing quadtree segmentation
 overwin=name

Name for output vector map showing overlapping windows
 nprocs=integer

Number of threads for parallel computing
Default: 1
 mask=name

Name of raster map used as mask
 tension=float

Tension parameter
Default: 40.
 smooth=float

Smoothing parameter
Smoothing is by default 0.5 unless smooth_column is specified
 smooth_column=string

Name of the attribute column with smoothing parameters
 segmax=integer

Maximum number of points in a segment
Default: 40
 npmin=integer

Minimum number of points for approximation in a segment (>segmax)
Default: 300
 dmin=float

Minimum distance between points (to remove almost identical points)
 dmax=float

Maximum distance between points on isoline (to insert additional
points)
 zscale=float

Conversion factor for values used for approximation
Default: 1.0
 theta=float

Anisotropy angle (in degrees counterclockwise from East)
 scalex=float

Anisotropy scaling factor
v.surf.rst program performs spatial approximation based on
zvalues (input vector map is 3D and
zcolumn parameter is not
given),
categories (input vector map is 2D and
zcolumn parameter
is not given), or
attributes (
zcolumn parameter is given) of
point or isoline data given in a vector map named
input to grid cells
in the output raster map
elevation representing a surface.
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 quadtreebased 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.
v.surf.rst uses regularized spline with tension for approximation from
vector data. The module does not require input data with topology, therefore
both level1 (no topology) and level2 (with topology) vector point data are
supported. Additional points are used for approximation between each 2 points
on a line if the distance between them is greater than specified
dmax.
If
dmax is small (less than cell size) the number of added data points
can be vary large and slow down approximation significantly. The
implementation has a segmentation procedure based on quadtrees which enhances
the efficiency for large data sets. Special color tables are created by the
program for output raster maps.
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 (090 and 1360 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):
 1
 Default: the given tension is applied to normalized data (
x/dnorm), that means that the distances are multiplied (rescaled)
by tension/dnorm. If density of points is changed, e.g., by using
higher dmin, the dnorm changes and tension needs to
be changed too to get the same result. Because the tension is
applied to normalized data its suitable value is usually within the 10100
range and does not depend on the actual scale (distances) of the original
data (which can be km for regional applications or cm for field
experiments).
 2
 Flagt: The given tension is applied to unnormalized data
(rescaled tension = tension*dnorm/1000 is applied to normalized
data ( x/dnorm) and therefore dnorm cancels out) so here
tension truly works as a rescaling parameter. For regional
applications with distances between points in km the suitable tension can
be 500 or higher, for detailed field scale analysis it can be 0.1. To help
select how much the data need to be rescaled the program writes
dnorm and rescaled tension fi=tension*dnorm/1000 at the
beginning of the program run. This rescaled tension should be
around 2030. If it is lower or higher, the given tension parameter
should be changed accordingly.
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
01, 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 from http://grass.osgeo.org/).
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 crossvalidation (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
"leaveoneout" or "jackknife" 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
timeconsuming 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 crossvalidation procedure works well only for wellsampled
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.
Using the
where parameter, the interpolation can be limited to use only a
subset of the input vectors.
North Carolina example (we simulate randomly distributed elevation measures
which we interpolate to a gapfree 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"
 •
 Mitasova, H., Mitas, L. and Harmon, R.S., 2005, Simultaneous spline
approximation and topographic analysis for lidar elevation data in open
source GIS, IEEE GRSL 2 (4), 375 379.
 •
 Hofierka, J., 2005, Interpolation of Radioactivity Data Using Regularized
Spline with Tension. Applied GIS, Vol. 1, No. 2, pp. 1601 to 1613. DOI:
10.2104/ag050016
 •
 Hofierka J., Parajka J., Mitasova H., Mitas L., 2002, Multivariate
Interpolation of Precipitation Using Regularized Spline with Tension.
Transactions in GIS 6(2), pp. 135150.
 •
 H. Mitasova, L. Mitas, B.M. Brown, D.P. Gerdes, I. Kosinovsky, 1995,
Modeling spatially and temporally distributed phenomena: New methods and
tools for GRASS GIS. International Journal of GIS, 9 (4), special issue on
Integrating GIS and Environmental modeling, 433446.
 •
 Mitasova, H. and Mitas, L., 1993: Interpolation by Regularized Spline with
Tension: I. Theory and Implementation, Mathematical Geology ,25,
641655.
 •
 Mitasova, H. and Hofierka, J., 1993: Interpolation by Regularized Spline
with Tension: II. Application to Terrain Modeling and Surface Geometry
Analysis, Mathematical Geology 25, 657667.
 •
 Mitas, L., and Mitasova H., 1988, General variational approach to the
approximation problem, Computers and Mathematics with Applications, v.16,
p. 983992.
 •
 Neteler, M. and Mitasova, H., 2008, Open Source GIS: A GRASS GIS Approach,
3rd Edition, Springer, New York, 406 pages.
 •
 Talmi, A. and Gilat, G., 1977 : Method for Smooth Approximation of Data,
Journal of Computational Physics, 23, p.93123.
 •
 Wahba, G., 1990, : Spline Models for Observational Data, CNMSNSF Regional
Conference series in applied mathematics, 59, SIAM, Philadelphia,
Pennsylvania.
v.vol.rst, v.surf.idw, v.surf.bspline,
r.fillnulls, g.region
Overview: Interpolation and Resampling in GRASS GIS
For examples of applications see GRASS4 implementation and GRASS5 and GRASS6
implementation.
Original version of program (in FORTRAN) and GRASS enhancements:
Lubos Mitas, NCSA, University of Illinois at Urbana Champaign, Illinois, USA
(19902000); Department of Physics, North Carolina State University, Raleigh
Helena Mitasova, USA CERL, Department of Geography, University of Illinois at
UrbanaChampaign, USA (19902001); MEAS, North Carolina State University,
Raleigh
Modified program (translated to C, adapted for GRASS, new segmentation
procedure):
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, ITCirst
Parallelization using OpenMP:
Stanislav Zubal, Czech Technical University in Prague
Michal Lacko, Pavol Jozef Safarik University in Kosice
Last changed: $Date: 20180827 17:17:54 +0200 (Mon, 27 Aug 2018) $
Available at: v.surf.rst source code (history)
Main index  Vector index  Topics index  Keywords index  Graphical index 
Full index
© 20032018 GRASS Development Team, GRASS GIS 7.4.3 Reference
Manual