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