Provided by: grass-doc_7.6.0-1_all bug

NAME  - 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, no-data filling

SYNOPSIS --help  [-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]

           Perform cross-validation procedure without raster approximation

           Use scale dependent tension

           Output partial derivatives instead of topographic parameters

           Allow output files to overwrite existing files

           Print usage summary

           Verbose module output

           Quiet module output

           Force launching GUI dialog

       input=name [required]
           Name of input vector map
           Or data source for direct OGR access

           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

           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 z-coordinates are used.

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

           Name for output surface elevation raster map

           Name for output slope raster map

           Name for output aspect raster map

           Name for output profile curvature raster map

           Name for output tangential curvature raster map

           Name for output mean curvature raster map

           Name for output deviations vector point map

           Name for output cross-validation errors vector point map

           Name for output vector map showing quadtree segmentation

           Name for output vector map showing overlapping windows

           Number of threads for parallel computing
           Default: 1

           Name of raster map used as mask

           Tension parameter
           Default: 40.

           Smoothing parameter
           Smoothing is by default 0.5 unless smooth_column is specified

           Name of the attribute column with smoothing parameters

           Maximum number of points in a segment
           Default: 40

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

           Minimum distance between points (to remove almost identical points)

           Maximum distance between points on isoline (to insert additional points)

           Conversion factor for values used for approximation
           Default: 1.0

           Anisotropy angle (in degrees counterclockwise from East)

           Anisotropy scaling factor

DESCRIPTION  program  performs spatial approximation based on z-values (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, 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

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

       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 10-100 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      Flag-t: The given tension is applied to  un-normalized  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  20-30.  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

       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. 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.  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.  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

   Cross validation procedure
       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.


       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 gap-free elevation surface):
       g.region raster=elevation -p
       # random elevation extraction of 500 samplings
       r.random elevation vector_output=elevrand n=500 -c elevrand elevrand
       # interpolation based on all points 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) 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.  16-01  to  16-13.  DOI:

           ·   Hofierka J., Parajka J., Mitasova H., Mitas L., 2002,  Multivariate  Interpolation
               of Precipitation Using Regularized Spline with Tension.  Transactions in GIS 6(2),
               pp. 135-150.

           ·   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, 433-446.

           ·   Mitasova,  H.  and  Mitas,  L.,  1993:  Interpolation  by  Regularized Spline with
               Tension: I. Theory and Implementation, Mathematical Geology ,25, 641-655.

           ·   Mitasova, H. and Hofierka, J., 1993:  Interpolation  by  Regularized  Spline  with
               Tension:  II.  Application  to  Terrain  Modeling  and  Surface Geometry Analysis,
               Mathematical Geology 25, 657-667.

           ·   Mitas,  L.,  and  Mitasova  H.,  1988,   General  variational  approach   to   the
               approximation  problem,  Computers  and  Mathematics  with  Applications, v.16, p.

           ·   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.93-123.

           ·   Wahba, G., 1990,  :  Spline  Models  for  Observational  Data,  CNMS-NSF  Regional
               Conference series in applied mathematics, 59, SIAM, Philadelphia, Pennsylvania.


        v.vol.rst,,, r.fillnulls, g.region

       Overview: Interpolation and Resampling in GRASS GIS

       For   examples   of   applications   see  GRASS4  implementation  and  GRASS5  and  GRASS6


       Original version of program (in FORTRAN) and GRASS enhancements:
       Lubos Mitas, NCSA, University of Illinois at Urbana Champaign, Illinois, USA  (1990-2000);
       Department of Physics, North Carolina State University, Raleigh
       Helena   Mitasova,   USA   CERL,  Department  of  Geography,  University  of  Illinois  at
       Urbana-Champaign, USA (1990-2001); 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, ITC-irst

       Parallelization using OpenMP:
       Stanislav Zubal, Czech Technical University in Prague
       Michal Lacko, Pavol Jozef Safarik University in Kosice

       Last changed: $Date: 2018-08-27 17:17:20 +0200 (Mon, 27 Aug 2018) $


       Available at: source code (history)

       Main index | Vector index | Topics index | Keywords index | Graphical index | Full index

       © 2003-2019 GRASS Development Team, GRASS GIS 7.6.0 Reference Manual