1v.surf.rst(1) GRASS GIS User's Manual v.surf.rst(1)
2
6 v.surf.rst - Performs surface interpolation from vector points map by
7 splines.
8 Spatial approximation and topographic analysis from given point or iso‐
9 line data in vector format to floating point raster format using regu‐
10 larized spline with tension.
11
13 vector, surface, interpolation, splines, RST, 3D, no-data filling
14
16 v.surf.rst
17 v.surf.rst --help
18 v.surf.rst [-ctd] input=name [layer=string] [zcolumn=name]
19 [where=sql_query] [elevation=name] [slope=name] [aspect=name]
20 [pcurvature=name] [tcurvature=name] [mcurvature=name] [devia‐
21 tions=name] [cvdev=name] [treeseg=name] [overwin=name]
22 [nprocs=integer] [mask=name] [tension=float] [smooth=float]
23 [smooth_column=string] [segmax=integer] [npmin=integer]
24 [dmin=float] [dmax=float] [zscale=float] [theta=float]
25 [scalex=float] [--overwrite] [--help] [--verbose] [--quiet]
26 [--ui]
27 Flags:
29 -c
30 Perform cross-validation procedure without raster approximation
31 -t
33 Use scale dependent tension
34 -d
36 Output partial derivatives instead of topographic parameters
37 --overwrite
39 Allow output files to overwrite existing files
40 --help
42 Print usage summary
43 --verbose
45 Verbose module output
46 --quiet
48 Quiet module output
49 --ui
51 Force launching GUI dialog
52 Parameters:
54 input=name [required]
55 Name of input vector map
56 Or data source for direct OGR access
57 layer=string
59 Layer number or name
60 Vector features can have category values in different layers. This
61 number determines which layer to use. When used with direct OGR
62 access this is the layer name.
63 Default: 1
64 zcolumn=name
66 Name of the attribute column with values to be used for approxima‐
67 tion
68 If not given and input is 2D vector map then category values are
69 used. If input is 3D vector map then z-coordinates are used.
70 where=sql_query
72 WHERE conditions of SQL statement without ’where’ keyword
73 Example: income < 1000 and population >= 10000
74 elevation=name
76 Name for output surface elevation raster map
77 slope=name
79 Name for output slope raster map
80 aspect=name
82 Name for output aspect raster map
83 pcurvature=name
85 Name for output profile curvature raster map
86 tcurvature=name
88 Name for output tangential curvature raster map
89 mcurvature=name
91 Name for output mean curvature raster map
92 deviations=name
94 Name for output deviations vector point map
95 cvdev=name
97 Name for output cross-validation errors vector point map
98 treeseg=name
100 Name for output vector map showing quadtree segmentation
101 overwin=name
103 Name for output vector map showing overlapping windows
104 nprocs=integer
106 Number of threads for parallel computing
107 Default: 1
108 mask=name
110 Name of raster map used as mask
111 tension=float
113 Tension parameter
114 Default: 40.
115 smooth=float
117 Smoothing parameter
118 Smoothing is by default 0.5 unless smooth_column is specified
119 smooth_column=string
121 Name of the attribute column with smoothing parameters
122 segmax=integer
124 Maximum number of points in a segment
125 Default: 40
126 npmin=integer
128 Minimum number of points for approximation in a segment (>segmax)
129 Default: 300
130 dmin=float
132 Minimum distance between points (to remove almost identical points)
133 dmax=float
135 Maximum distance between points on isoline (to insert additional
136 points)
137 zscale=float
139 Conversion factor for values used for approximation
140 Default: 1.0
141 theta=float
143 Anisotropy angle (in degrees counterclockwise from East)
144 scalex=float
146 Anisotropy scaling factor
147
149 v.surf.rst program performs spatial approximation based on z-values
150 (input vector map is 3D and zcolumn parameter is not given), categories
151 (input vector map is 2D and zcolumn parameter is not given), or
152 attributes (zcolumn parameter is given) of point or isoline data given
153 in a vector map named input to grid cells in the output raster map ele‐
154 vation representing a surface.
155 As an option, simultaneously with approximation, topographic parameters
157 slope, aspect, profile curvature (measured in the direction of the
158 steepest slope), tangential curvature (measured in the direction of a
159 tangent to contour line) or mean curvature are computed and saved as
160 raster maps specified by the options slope, aspect, pcurv, tcurv, mcurv
161 respectively. If -d flag is set, v.surf.rst outputs partial derivatives
162 fx,fy,fxx, fyy,fxy instead of slope, aspect, profile, tangential and
163 mean curvatures respectively. If the input vector map have time stamp,
164 the program creates time stamp for all output maps.
165 User can either use r.mask to set a mask or specify a raster map in
167 mask option, which will be used as a mask. The approximation is skipped
168 for cells which have zero or NULL value in mask. NULL values will be
169 assigned to these cells in all output raster maps. Data points are
170 checked for identical points and points that are closer to each other
171 than the given dmin are removed. If sparsely digitized contours or
172 isolines are used as input, additional points are computed between each
173 2 points on a line if the distance between them is greater than speci‐
174 fied dmax. Parameter zmult allows user to rescale the values used for
175 approximation (useful e.g. for transformation of elevations given in
176 feet to meters, so that the proper values of slopes and curvatures can
177 be computed).
178 Regularized spline with tension is used for the approximation. The ten‐
180 sion parameter tunes the character of the resulting surface from thin
181 plate to membrane. Smoothing parameter smooth controls the deviation
182 between the given points and the resulting surface and it can be very
183 effective in smoothing noisy data while preserving the geometrical
184 properties of the surface. With the smoothing parameter set to zero
185 (smooth=0) the resulting surface passes exactly through the data points
186 (spatial interpolation is performed). When smoothing parameter is used,
187 it is also possible to output a vector point map deviations containing
188 deviations of the resulting surface from the given data.
189 If the number of given points is greater than segmax, segmented pro‐
191 cessing is used. The region is split into quadtree-based rectangular
192 segments, each having less than segmax points and approximation is per‐
193 formed on each segment of the region. To ensure smooth connection of
194 segments the approximation function for each segment is computed using
195 the points in the given segment and the points in its neighborhood
196 which are in the rectangular window surrounding the given segment. The
197 number of points taken for approximation is controlled by npmin, the
198 value of which must be larger than segmax. User can choose to output
199 vector maps treeseg and overwin which represent the quad tree used for
200 segmentation and overlapping neighborhoods from which additional points
201 for approximation on each segment were taken.
202 Predictive error of surface approximation for given parameters can be
204 computed using the -c flag. A crossvalidation procedure is then per‐
205 formed using the data given in the vector map input and the estimated
206 predictive errors are stored in the vector point map cvdev. When using
207 this flag, no raster output maps are computed. Anisotropic surfaces
208 can be interpolated by setting anisotropy angle theta and scaling fac‐
209 tor scalex. The program writes values of selected input and internally
210 computed parameters to the history file of raster map elevation.
211 The user must run g.region before the program to set the region and
213 resolution for approximation.
214
216 v.surf.rst uses regularized spline with tension for approximation from
217 vector data. The module does not require input data with topology,
218 therefore both level1 (no topology) and level2 (with topology) vector
219 point data are supported. Additional points are used for approximation
220 between each 2 points on a line if the distance between them is greater
221 than specified dmax. If dmax is small (less than cell size) the number
222 of added data points can be vary large and slow down approximation sig‐
223 nificantly. The implementation has a segmentation procedure based on
224 quadtrees which enhances the efficiency for large data sets. Special
225 color tables are created by the program for output raster maps.
226 Topographic parameters are computed directly from the approximation
228 function so that the important relationships between these parameters
229 are preserved. The equations for computation of these parameters and
230 their interpretation is described in Mitasova and Hofierka, 1993 or
231 Neteler and Mitasova, 2004). Slopes and aspect are computed in degrees
232 (0-90 and 1-360 respectively). The aspect raster map has value 0
233 assigned to flat areas (with slope less than 0.1%) and to singular
234 points with undefined aspect. Aspect points downslope and is 90 to the
235 North, 180 to the West, 270 to the South and 360 to the East, the val‐
236 ues increase counterclockwise. Curvatures are positive for convex and
237 negative for concave areas. Singular points with undefined curvatures
238 have assigned zero values.
239 Tension and smoothing allow user to tune the surface character. For
241 most landscape scale applications the default values should provide
242 adequate results. The program gives warning when significant over‐
243 shoots appear in the resulting surface and higher tension or smoothing
244 should be used.
245 To select parameters that will produce a surface with desired proper‐
247 ties, it is useful to know that the method is scale dependent and the
248 tension works as a rescaling parameter (high tension "increases the
249 distances between the points" and reduces the range of impact of each
250 point, low tension "decreases the distance" and the points influence
251 each other over longer range). Surface with tension set too high
252 behaves like a membrane (rubber sheet stretched over the data points)
253 with peak or pit ("crater") in each given point and everywhere else the
254 surface goes rapidly to trend. If digitized contours are used as input
255 data, high tension can cause artificial waves along contours. Lower
256 tension and higher smoothing is suggested for such a case.
257 Surface with tension set too low behaves like a stiff steel plate and
259 overshoots can appear in areas with rapid change of gradient and seg‐
260 mentation can be visible. Increase in tension should solve the prob‐
261 lems.
262 There are two options how tension can be applied in relation to dnorm
264 (dnorm rescales the coordinates depending on the average data density
265 so that the size of segments with segmax=40 points is around 1 - this
266 ensures the numerical stability of the computation):
267 1 Default: the given tension is applied to normalized data
269 (x/dnorm), that means that the distances are multiplied
270 (rescaled) by tension/dnorm. If density of points is changed,
271 e.g., by using higher dmin, the dnorm changes and tension needs
272 to be changed too to get the same result. Because the tension
273 is applied to normalized data its suitable value is usually
274 within the 10-100 range and does not depend on the actual scale
275 (distances) of the original data (which can be km for regional
276 applications or cm for field experiments).
277 2 Flag-t: The given tension is applied to un-normalized data
279 (rescaled tension = tension*dnorm/1000 is applied to normalized
280 data (x/dnorm) and therefore dnorm cancels out) so here tension
281 truly works as a rescaling parameter. For regional applications
282 with distances between points in km the suitable tension can be
283 500 or higher, for detailed field scale analysis it can be 0.1.
284 To help select how much the data need to be rescaled the program
285 writes dnorm and rescaled tension fi=tension*dnorm/1000 at the
286 beginning of the program run. This rescaled tension should be
287 around 20-30. If it is lower or higher, the given tension param‐
288 eter should be changed accordingly.
289 The default is a recommended choice, however for the applications where
291 the user needs to change density of data and preserve the approximation
292 character the -t flag can be helpful.
293 Anisotropic data (e.g. geologic phenomena) can be interpolated using
295 theta and scalex defining orientation and ratio of the perpendicular
296 axes put on the longest/shortest side of the feature, respectively.
297 Theta is measured in degrees from East, counterclockwise. Scalex is a
298 ratio of axes sizes. Setting scalex in the range 0-1, results in a
299 pattern prolonged in the direction defined by theta. Scalex value 0.5
300 means that modeled feature is approximately 2 times longer in the
301 direction of theta than in the perpendicular direction. Scalex value 2
302 means that axes ratio is reverse resulting in a pattern perpendicular
303 to the previous example. Please note that anisotropy option has not
304 been extensively tested and may include bugs (for example, topographic
305 parameters may not be computed correctly) - if there are problems,
306 please report to GRASS bugtracker (accessible from
307 https://grass.osgeo.org/).
308 For data with values changing over several magnitudes (sometimes the
310 concentration or density data) it is suggested to interpolate the log
311 of the values rather than the original ones.
312 v.surf.rst checks the numerical stability of the algorithm by computing
314 the values in given points, and prints the root mean square deviation
315 (rms) found into the history file of raster map elevation. For computa‐
316 tion with smoothing set to 0, rms should be 0. Significant increase in
317 tension is suggested if the rms is unexpectedly high for this case.
318 With smoothing parameter greater than zero the surface will not pass
319 exactly through the data points and the higher the parameter the closer
320 the surface will be to the trend. The rms then represents a measure of
321 smoothing effect on data. More detailed analysis of smoothing effects
322 can be performed using the output deviations option.
323 v.surf.rst also writes the values of parameters used in computation
325 into the comment part of history file elevation as well as the follow‐
326 ing values which help to evaluate the results and choose the suitable
327 parameters: minimum and maximum z values in the data file (zmin_data,
328 zmax_data) and in the interpolated raster map (zmin_int, zmax_int),
329 rescaling parameter used for normalization (dnorm), which influences
330 the tension.
331 If visible connection of segments appears, the program should be rerun
333 with higher npmin to get more points from the neighborhood of given
334 segment and/or with higher tension.
335 When the number of points in a vector map is not too large (less than
337 800), the user can skip segmentation by setting segmax to the number of
338 data points or segmax=700.
339 v.surf.rst gives warning when user wants to interpolate outside the
341 rectangle given by minimum and maximum coordinates in the vector map,
342 zoom into the area where the given data are is suggested in this case.
343 When a mask is used, the program takes all points in the given region
345 for approximation, including those in the area which is masked out, to
346 ensure proper approximation along the border of the mask. It therefore
347 does not mask out the data points, if this is desirable, it must be
348 done outside v.surf.rst.
349 Cross validation procedure
351 The "optimal" approximation parameters for given data can be found
352 using a cross-validation (CV) procedure (-c flag). The CV procedure is
353 based on removing one input data point at a time, performing the
354 approximation for the location of the removed point using the remaining
355 data points and calculating the difference between the actual and
356 approximated value for the removed data point. The procedure is
357 repeated until every data point has been, in turn, removed. This form
358 of CV is also known as the "leave-one-out" or "jack-knife" method
359 (Hofierka et al., 2002; Hofierka, 2005). The differences (residuals)
360 are then stored in the cvdev output vector map. Please note that during
361 the CV procedure no other output maps can be set, the approximation is
362 performed only for locations defined by input data. To find "optimal
363 parameters", the CV procedure must be iteratively performed for all
364 reasonable combinations of the approximation parameters with small
365 incremental steps (e.g. tension, smoothing) in order to find a combina‐
366 tion with minimal statistical error (also called predictive error)
367 defined by root mean squared error (RMSE), mean absolute error (MAE) or
368 other error characteristics. A script with loops for tested RST param‐
369 eters can do the job, necessary statistics can be calculated using e.g.
370 v.univar. It should be noted that crossvalidation is a time-consuming
371 procedure, usually reasonable for up to several thousands of points.
372 For larger data sets, CV should be applied to a representative subset
373 of the data. The cross-validation procedure works well only for
374 well-sampled phenomena and when minimizing the predictive error is the
375 goal. The parameters found by minimizing the predictive (CV) error may
376 not not be the best for for poorly sampled phenomena (result could be
377 strongly smoothed with lost details and fluctuations) or when signifi‐
378 cant noise is present that needs to be smoothed out.
379
381 Setting for lidar point cloud
382 Lidar point clouds as well as UAS SfM-based (phodar) point clouds tend
383 to be dense in relation to the desired raster resolution and thus a
384 different set of parameters is more advantageous, e.g. in comparison to
385 a typical temperature data interpolation.
386 v.surf.rst input=points elevation=elevation npmin=100
387 Usage of the where parameter
389 Using the where parameter, the interpolation can be limited to use only
390 a subset of the input vectors.
391 North Carolina example (we simulate randomly distributed elevation mea‐
393 sures which we interpolate to a gap-free elevation surface):
394 g.region raster=elevation -p
395 # random elevation extraction of 500 samplings
396 r.random elevation vector_output=elevrand n=500
397 v.info -c elevrand
398 v.db.select elevrand
399 # interpolation based on all points
400 v.surf.rst elevrand zcol=value elevation=elev_full
401 # apply the color table of the original raster map
402 r.colors elev_full raster=elevation
403 d.rast elev_full
404 d.vect elevrand
405 # interpolation based on subset of points (only those over 1300m/asl)
406 v.surf.rst elevrand zcol=value elevation=elev_partial where="value > 1300"
407 r.colors elev_partial raster=elevation
408 d.rast elev_partial
409 d.vect elevrand where="value > 1300"
410
412 · Mitasova, H., Mitas, L. and Harmon, R.S., 2005, Simultaneous
413 spline approximation and topographic analysis for lidar eleva‐
414 tion data in open source GIS, IEEE GRSL 2 (4), 375- 379.
415 · Hofierka, J., 2005, Interpolation of Radioactivity Data Using
417 Regularized Spline with Tension. Applied GIS, Vol. 1, No. 2,
418 pp. 16-01 to 16-13. DOI: 10.2104/ag050016
419 · Hofierka J., Parajka J., Mitasova H., Mitas L., 2002, Multi‐
421 variate Interpolation of Precipitation Using Regularized Spline
422 with Tension. Transactions in GIS 6(2), pp. 135-150.
423 · H. Mitasova, L. Mitas, B.M. Brown, D.P. Gerdes, I. Kosinovsky,
425 1995, Modeling spatially and temporally distributed phenomena:
426 New methods and tools for GRASS GIS. International Journal of
427 GIS, 9 (4), special issue on Integrating GIS and Environmental
428 modeling, 433-446.
429 · Mitasova, H. and Mitas, L., 1993: Interpolation by Regularized
431 Spline with Tension: I. Theory and Implementation, Mathematical
432 Geology ,25, 641-655.
433 · Mitasova, H. and Hofierka, J., 1993: Interpolation by Regular‐
435 ized Spline with Tension: II. Application to Terrain Modeling
436 and Surface Geometry Analysis, Mathematical Geology 25,
437 657-667.
438 · Mitas, L., and Mitasova H., 1988, General variational approach
440 to the approximation problem, Computers and Mathematics with
441 Applications, v.16, p. 983-992.
442 · Neteler, M. and Mitasova, H., 2008, Open Source GIS: A GRASS
444 GIS Approach, 3rd Edition, Springer, New York, 406 pages.
445 · Talmi, A. and Gilat, G., 1977 : Method for Smooth Approximation
447 of Data, Journal of Computational Physics, 23, p.93-123.
448 · Wahba, G., 1990, : Spline Models for Observational Data,
450 CNMS-NSF Regional Conference series in applied mathematics, 59,
451 SIAM, Philadelphia, Pennsylvania.
452
454 v.vol.rst, v.surf.idw, v.surf.bspline, r.fillnulls, g.region
455 Overview: Interpolation and Resampling in GRASS GIS
457 For examples of applications see GRASS4 implementation and GRASS5 and
459 GRASS6 implementation.
460
462 Original version of program (in FORTRAN) and GRASS enhancements:
463 Lubos Mitas, NCSA, University of Illinois at Urbana Champaign, Illi‐
464 nois, USA (1990-2000); Department of Physics, North Carolina State Uni‐
465 versity, Raleigh
466 Helena Mitasova, USA CERL, Department of Geography, University of Illi‐
467 nois at Urbana-Champaign, USA (1990-2001); MEAS, North Carolina State
468 University, Raleigh
469 Modified program (translated to C, adapted for GRASS, new segmentation
471 procedure):
472 Irina Kosinovsky, US Army CERL, Dave Gerdes, US Army CERL
473 Modifications for new sites format and timestamping:
475 Darrel McCauley, Purdue University, Bill Brown, US Army CERL
476 Update for GRASS5.7, GRASS6 and addition of crossvalidation:
478 Jaroslav Hofierka, University of Presov; Radim Blazek, ITC-irst
479 Parallelization using OpenMP:
481 Stanislav Zubal, Czech Technical University in Prague
482 Michal Lacko, Pavol Jozef Safarik University in Kosice
483
485 Available at: v.surf.rst source code (history)
486 Main index | Vector index | Topics index | Keywords index | Graphical
488 index | Full index
489 © 2003-2020 GRASS Development Team, GRASS GIS 7.8.5 Reference Manual
491GRASS 7.8.5 v.surf.rst(1)