|
@@ -1,78 +1,89 @@
|
|
|
<h2>DESCRIPTION</h2>
|
|
|
-<em>v.surf.bspline</em> makes a bilinear/bicubic spline interpolation
|
|
|
-with Tykhonov regularization. The required input is an only 3d points
|
|
|
-vector map that will be used to interpolate a reference surface.
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-Interpolation is carried out by adjusting a Least-Squares (LS) system
|
|
|
-in which the parameters to estime are spline functions. The number of
|
|
|
-splines doesn't depend on the resolution region, but it depends on the
|
|
|
-spline steps values in the north-south and west-east directions. These
|
|
|
-spline steps are set by "<b><i>sin=</i></b>" and "<b><i>sie=</i></b>",
|
|
|
-respectively. If the number of splines is bigger than the number of
|
|
|
-points, the LS system is bad conditioned because there are more unkowns
|
|
|
-than observations. In that case the LS normal matrix can't be inverted.
|
|
|
-To allow the inversion of the normal matrix a Tykhonov regularization
|
|
|
-is done. The minimizing function is the gradient in the case of a bilinear
|
|
|
-interpolation, and the curvature in the bicubic interpolation. The
|
|
|
-lambda_i parameter associated with the regularization smooths the
|
|
|
-interpolation. The higher the lambda_i parameter, the smoother the
|
|
|
+<em>v.surf.bspline</em> performs a bilinear/bicubic spline interpolation with
|
|
|
+Tykhonov regularization. The input is a 2D or 3D vector points map. Values to
|
|
|
+interpolate can be the z values of 3D points or the values in a user-specified
|
|
|
+attribue column in a 2D or 3D map. Output can be a raster or vector map.
|
|
|
+Optionally, a "sparse point" vector map can be input specify vector points
|
|
|
+output.
|
|
|
+<br> <br>
|
|
|
+From a theoretical perspective, the interpolating procedure takes place in two
|
|
|
+parts: the first is an estimate of the linear coefficients of a spline function
|
|
|
+is derived from the observation points using a least squares regression; the
|
|
|
+second is the computation of the interpolated surface (or interpolated vector
|
|
|
+points). As used here, the splines are 2D piece-wise non-zero polynomial
|
|
|
+functions calculated within a limited, 2D area. The length of each spline step
|
|
|
+is defined by <b><i>sie</i></b> for the east-west direction and
|
|
|
+<b><i>sin</i></b> for the north-south direction. For optimum performance, the
|
|
|
+length of spline step should be no less than the distance between observation
|
|
|
+points. Each vector point observation is modeled as a linear function of the
|
|
|
+non-zero splines in the area around the observation. The least squares
|
|
|
+regression predicts the the coefficients of these linear functions.
|
|
|
+Regularization, avoids the need to have one one observation and one coefficient
|
|
|
+for each spline (in order to avoid instability).
|
|
|
+
|
|
|
+<p>
|
|
|
+With regularly distributed data points, a spline step corresponding to the
|
|
|
+maximum distance between two points in both the east and north directions is
|
|
|
+sufficient. But often data points are not regularly distributed and require
|
|
|
+statistial regularization or estimation. In such cases, v.surf.bspline will
|
|
|
+attempt to minimize the gradient of bilinear splines or the curvature of bicubic
|
|
|
+splines in areas lacking point observations. As a general rule, spline step
|
|
|
+length should be greater than the mean distance between observation points
|
|
|
+(twice the distance between points is a good starting point). Separate east-west
|
|
|
+and north-south spline step length arguments allows the user to account for some
|
|
|
+degree of anisotropy in the distribution of observation points. Short spline
|
|
|
+step lengths--especially spline step lengths that are less than the distance
|
|
|
+between observation points--can greatly increase processing time.
|
|
|
+
|
|
|
+<p>
|
|
|
+Moreover, the maximum number of splines for each direction at each time is
|
|
|
+fixed, regardless of the spline step length. As the total number of splines used
|
|
|
+increases (i.e., with small spline step lengths), the region is automatically
|
|
|
+into subregions for interpolation. Each subregion can contain no more than
|
|
|
+150x150 splines. To avoid subregion boundary problems, subregions are created to
|
|
|
+partially overlap each other. A weighted mean of observations, based on point
|
|
|
+locations, is calculated within each subregion.
|
|
|
+
|
|
|
+<p>
|
|
|
+The Tykhonov regularization parameter ("<b><i>lambda_i</i></b>") acts to smooth
|
|
|
+the interpolation. With a small <b><i>lambda_i</i></b>, the interpolated surface
|
|
|
+closely follows observation points; a larger value will produce a smoother
|
|
|
interpolation.
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-The number of splines has a great influence on two things, mainly. The
|
|
|
-first thing is the module's execution time. The second is the RAM use.
|
|
|
-The higher the number of splines, the longer the time of execution and
|
|
|
-the higher RAM use. A numerical example: 100 splines in each direction
|
|
|
-imply 10e4 splines in total, that is, a square LS normal matrix of 10e4
|
|
|
-size. Inverting this matrix means inverting 100 millions elements!
|
|
|
-To improve this problems a Tcholebsky method with triangulars matrixes
|
|
|
-is used in the normal matrix inversion. It has also fixed a maximum number
|
|
|
-of splines for each direction. However, it is also possible running the
|
|
|
-module with a higher number of splines. For a number of spline higher than
|
|
|
-the fixed maximum, the whole region is divided into smaller regions. Each
|
|
|
-subregion is 150x150 splines wide. To avoid contour problems, the subregions
|
|
|
-are overlaped one to each other. To estimate a single value within the
|
|
|
-overlaped zones, a weighted mean considering the point positions into each
|
|
|
-subregion is carried out.
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-The required input is a 3d points vector. If nothing is specified z-coordinates
|
|
|
-will be used in the interpolation. It could be also possible to consider
|
|
|
-an attribute value by specifying "<b><i>layer=</i></b>" and "<b><i>column=</i></b>"
|
|
|
-parameters. If a vector map with another type of features is used, only
|
|
|
-points will be considered. If the "<b><i>sparse=</i></b>" vector is
|
|
|
-used, the "<b><i>input=</i></b>" vector map will be used to create a
|
|
|
-reference surface. This surface will be used to make an estimation on the
|
|
|
-points within the "<b><i>sparse=</i></b>". In this case a vector output
|
|
|
-("<b><i>output=</i></b>") must be specify. If the "<b><i>sparse=</i></b>"
|
|
|
-is not supplied, the final interpolation output will be the interpolated
|
|
|
-reference surface from the "<b><i>input=</i></b>" vector map. In this case,
|
|
|
-one of both the raster or vector output format can be choosen. For raster
|
|
|
-format ("<b><i>raster=</i></b>"), the point estimation will be done
|
|
|
-on a regular grid with a resolution equal to the GRASS region. For vector
|
|
|
-format, the estimation will be done on the sparse points of the
|
|
|
-"<b><i>input=</i></b>" vector supplied. Both, vector and raster output,
|
|
|
-are not allowed simultaneously.
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-A cross validation method has been implemented. It helps to find the optimal
|
|
|
-lambda_i value that fits the data. It shows the <i>mean</i> and <i>rms</i>
|
|
|
-of the residuals from the true point value and the estimated from the
|
|
|
-interpolation made with all the data without the point itself. This procedure
|
|
|
-is done for fixed lambda_i values. The results of the cross validation will
|
|
|
-appear in the stdout and no vector nor raster output will be created. The
|
|
|
-external input ("<b><i>sparse=</i></b>") will be not considered. Due to
|
|
|
-the nature of the algorithm, it is advised the user no to try the cross-
|
|
|
-validation with more than 100 points at a time because it will take too long.
|
|
|
-The execution time could be reduced by considering a lower number of splines.
|
|
|
-Although, as seen, it is possible to use a high number of splines, more than
|
|
|
-150x150 splines is not recommended.
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-In a raster map output ("<b><i>raster=</i></b>"), region resolution implying
|
|
|
-more than 2000x2000 (4 mill) cells are not allowed. If the user tries with a
|
|
|
-more than those cells an error message will ask for a lower region resolution.
|
|
|
+
|
|
|
+<p>
|
|
|
+The input can be a 2D pr 3D vector points map. If "<b><i>layer =</i></b>" 0 the
|
|
|
+z-value of a 3D map is used for interpolation. If layer > 0, the user must
|
|
|
+specify an attribute column to used for interpolation using the
|
|
|
+"<b><i>column=</i></b>" argument (2D or 3D map).
|
|
|
+
|
|
|
+<p>
|
|
|
+v.surf.bspline can produce a raster OR a vector output (NOT simultaneously).
|
|
|
+
|
|
|
+<p>
|
|
|
+If output is a vector points map and a "<b><i>sparse=</i></b>" vector points map
|
|
|
+is not specified, the output vector map will contain points at the same
|
|
|
+locations as observation points in the input map, but the values of the output
|
|
|
+points are interpolated values. If a "<b><i>sparse=</i></b>" vector points map
|
|
|
+is specified, the output vector map will contain points at the same locations as
|
|
|
+the sparse vector map points, and values will be those of the interpolated
|
|
|
+raster surface at those points.
|
|
|
+
|
|
|
+<p>
|
|
|
+A cross validation "leave-one-out" analysis is available to help to determine
|
|
|
+the optimal <b><i>lambda_i</i></b> value that produces an interpolation that
|
|
|
+best fits the original observation data. The more points used for
|
|
|
+cross-validation, the longer the time needed for computation. Empirical testing
|
|
|
+indicates a threshold of a maximum of 100 points is recommended. The
|
|
|
+cross-validation output reports <i>mean</i> and <i>rms</i> of the residuals from
|
|
|
+the true point value and the estimated from the interpolation for a fixed series
|
|
|
+of <b><i>lambda_i</i></b> values. No vector nor raster output will be created
|
|
|
+when cross-validation is selected.
|
|
|
+
|
|
|
+<p>
|
|
|
+A raster output map ("<b><i>raster=</i></b>") of more than 2000x2000 (4 mill)
|
|
|
+cells is not allowed. If an output map would exceed this size, an error message
|
|
|
+is generated.
|
|
|
+
|
|
|
|
|
|
<h2>EXAMPLES</h2>
|
|
|
|
|
@@ -82,85 +93,90 @@ more than those cells an error message will ask for a lower region resolution.
|
|
|
v.surf.bspline input=point_vector output=interpolate_surface type=bicubic
|
|
|
</pre></div>
|
|
|
|
|
|
-In this case, a bicubic spline interpolation will be done and an
|
|
|
-estimation on the points of point_vector will be the output.
|
|
|
+A bicubic spline interpolation will be done and a vector points map with estimated
|
|
|
+(i.e., interpolated) values will be created.
|
|
|
|
|
|
-<h4>Basic interpolation and raster output with a long spline step</h4>
|
|
|
+<h4>Basic interpolation and raster output with a longer spline step</h4>
|
|
|
|
|
|
<div class="code"><pre>
|
|
|
v.surf.bspline input=point_vector raster=interpolate_surface sie=25 sin=25
|
|
|
</pre></div>
|
|
|
|
|
|
-Now, a bilinear spline interpolation will be done on a grid. The spline steps
|
|
|
-are set to 25. It doesn't mean that the grid will have a resolution equal to 25,
|
|
|
-but that each 25 units there will be a spline.
|
|
|
+A bilinear spline interpolation will be done with a spline step length of 25 map
|
|
|
+units. An interpolated raster map will be created at the current region resolution.
|
|
|
|
|
|
-<h4> Estimation of lambda_i parameter with a cross validation proccess</h4>
|
|
|
+<h4>Estimation of <b><i>lambda_i</i></b> parameter with a cross validation proccess</h4>
|
|
|
|
|
|
<div class="code"><pre>
|
|
|
v.surf.bspline -c input=point_vector
|
|
|
</pre></div>
|
|
|
|
|
|
-
|
|
|
<h4>Estimation on sparse points</h4>
|
|
|
|
|
|
<div class="code"><pre>
|
|
|
v.surf.bspline input=point_vector sparse=sparse_points output=interpolate_surface
|
|
|
</pre></div>
|
|
|
|
|
|
-In this last case, an estimation on the points of the sparse_points vector
|
|
|
-will be done. The reference surface used for this estimation will be that
|
|
|
-interpolated using the point_vector vector.
|
|
|
+An output map of vector points will be created, corresponding to the sparse vector map, with interpolated values.
|
|
|
|
|
|
<h4>Using attribute values instead Z-coordinates</h4>
|
|
|
<div class="code"><pre>
|
|
|
v.surf.bspline input=point_vector raster=interpolate_surface layer=1 column=attrib_column
|
|
|
</pre></div>
|
|
|
|
|
|
-This last case, the module uses the attribute values in attrib_column
|
|
|
-in the table associated to layer 1.
|
|
|
+The interpolation will be done using the values in attrib_column, in the
|
|
|
+table associated with layer 1.
|
|
|
|
|
|
<h2>BUGS</h2>
|
|
|
Known issues:
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-In order to avoid RAM memory problems, an auxiliar table will be needed for
|
|
|
-recording some intermediate calculi. Since the "<b>GROUP BY</b>" SQL function is used,
|
|
|
-which is not supported by the "<b>dbf</b>" driver, this driver is not
|
|
|
-allowed with the vector map output "<b><i>output=</i></b>". There is no problem
|
|
|
-with the raster map output.
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-At this time, using the external vector input ("<b><i>sparse=</i></b>") implies
|
|
|
-interpoling with Z-coordinates. Updates to allow using attribute values
|
|
|
-will be done in a near future (I hope).
|
|
|
-<br>
|
|
|
-<br>
|
|
|
+<p>
|
|
|
+In order to avoid RAM memory problems, an auxiliary table is needed for
|
|
|
+recording some intermediate calculations. This requires the "<b>GROUP BY</b>"
|
|
|
+SQL function is used, which is not supported by the "<b>dbf</b>" driver. For
|
|
|
+this reason, vector map output "<b><i>output=</i></b>" is not permitted with the
|
|
|
+DBF driver. There are no problems with the raster map output from the DBF
|
|
|
+driver.
|
|
|
+
|
|
|
+<p>
|
|
|
+At this time, sparse vector input ("<b><i>sparse=</i></b>") can only be used
|
|
|
+with interpolation from 3D vector z-coordinates (<a href="http://trac.osgeo.org/grass/ticket/96">trac #96</a>).
|
|
|
+<p>
|
|
|
|
|
|
|
|
|
<h2>SEE ALSO</h2>
|
|
|
-<em><a HREF="v.surf.rst.html">v.surf.rst</a></em>
|
|
|
+
|
|
|
+<em><a
|
|
|
+href="http://grass.osgeo.org/grass64/manuals/html64_user/v.surf.idw.html">v.surf.idw</a>,
|
|
|
+<a
|
|
|
+href="http://grass.osgeo.org/grass64/manuals/html64_user/v.surf.rst.html">v.surf.rst</a></em>
|
|
|
|
|
|
<h2>AUTHORS</h2>
|
|
|
+
|
|
|
Original version in GRASS 5.4: (s.bspline.reg)
|
|
|
-<BR>
|
|
|
+<br>
|
|
|
Maria Antonia Brovelli, Massimiliano Cannata, Ulisse Longoni, Mirko Reguzzoni
|
|
|
-<BR><BR>
|
|
|
+<p>
|
|
|
Update for GRASS 6.X and improvements:
|
|
|
-<BR>
|
|
|
+<br>
|
|
|
Roberto Antolin
|
|
|
|
|
|
<h2>REFERENCES</h2>
|
|
|
-Brovelli M. A., Cannata M., and Longoni U.M., 2004, LIDAR Data Filtering and DTM Interpolation Within GRASS, Transactions in GIS, April 2004, vol. 8, iss. 2, pp. 155-174(20), Blackwell Publishing Ltd
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-Brovelli M. A. and Cannata M., 2004, Digital Terrain model reconstruction in urban areas from airborne laser scanning data: the method and an example for Pavia (Northern Italy). Computers and Geosciences 30, pp.325-331
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-Brovelli M. A e Longoni U.M., 2003, Software per il filtraggio di dati LIDAR, Rivista dell'Agenzia del Territorio, n. 3-2003, pp. 11-22 (ISSN 1593-2192)
|
|
|
-<br>
|
|
|
-<br>
|
|
|
-Brovelli M. A., Cannata M. and Longoni U.M., 2002, DTM LIDAR in area urbana, Bollettino SIFET N.2, 2002, pp. 7-26
|
|
|
-<br>
|
|
|
+
|
|
|
+Brovelli M. A., Cannata M., and Longoni U.M., 2004, LIDAR Data
|
|
|
+Filtering and DTM Interpolation Within GRASS, Transactions in GIS,
|
|
|
+April 2004, vol. 8, iss. 2, pp. 155-174(20), Blackwell Publishing Ltd
|
|
|
+<p>
|
|
|
+Brovelli M. A. and Cannata M., 2004, Digital Terrain model
|
|
|
+reconstruction in urban areas from airborne laser scanning data: the
|
|
|
+method and an example for Pavia (Northern Italy). Computers and
|
|
|
+Geosciences 30, pp.325-331
|
|
|
+<p>
|
|
|
+Brovelli M. A e Longoni U.M., 2003, Software per il filtraggio di
|
|
|
+dati LIDAR, Rivista dell'Agenzia del Territorio, n. 3-2003, pp. 11-22
|
|
|
+(ISSN 1593-2192)
|
|
|
+<p>
|
|
|
+Antolin R. and Brovelli M.A., 2007, LiDAR data Filtering with GRASS GIS for the
|
|
|
+Determination of Digital Terrain Models. Proceedings of Jornadas de SIG Libre,
|
|
|
+Girona, España. CD ISBN: 978-84-690-3886-9 <br>
|
|
|
|
|
|
<p><i>Last changed: $Date$</i>
|