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