GRASS GIS: v.surf.bspline

NAME

v.surf.bspline - Bicubic or bilinear spline interpolation with Tykhonov regularization.

KEYWORDS

vector, interpolation

SYNOPSIS

v.surf.bspline
v.surf.bspline help
v.surf.bspline [-c] input=name [sparse=name] [output=name] [raster=name] [sie=float] [sin=float] [type=string] [lambda_i=float] [layer=integer] [column=string] [--overwrite] [--verbose] [--quiet]

Flags:

-c: Find best parameters using a cross validation method
--overwrite: Allow output files to overwrite existing files
--verbose: Verbose module output
--quiet: Quiet module output

Parameters:

input=name: Name of input vector map
sparse=name: Name of input vector map of sparse points
output=name: Name for output vector map
raster=name: Name for output raster map
sie=float: Interpolation spline step value in east direction; Default: 4
sin=float: Interpolation spline step value in north direction; Default: 4
type=string: Spline type of interpolation; Options: bilinear,bicubic; Default: bilinear
lambda_i=float: Thychonov regularization weigth; Default: 1
layer=integer: Layer number; Field value. If set to 0, z coordinates are used. (3D vector only); Default: 0
column=string: Attribute table column with values to interpolate (if layer>0)

DESCRIPTION

v.surf.bspline 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.

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 "sin=" and "sie=", 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 interpolation.

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.

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 "layer=" and "column=" parameters. If a vector map with another type of features is used, only points will be considered. If the "sparse=" vector is used, the "input=" vector map will be used to create a reference surface. This surface will be used to make an estimation on the points within the "sparse=". In this case a vector output ("output=") must be specify. If the "sparse=" is not supplied, the final interpolation output will be the interpolated reference surface from the "input=" vector map. In this case, one of both the raster or vector output format can be choosen. For raster format ("raster="), 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 "input=" vector supplied. Both, vector and raster output, are not allowed simultaneously.

A cross validation method has been implemented. It helps to find the optimal lambda_i value that fits the data. It shows the mean and rms 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 ("sparse=") 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.

In a raster map output ("raster="), 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.

EXAMPLES

Basic interpolation

v.surf.bspline input=point_vector output=interpolate_surface type=bicubic

In this case, a bicubic spline interpolation will be done and an estimation on the points of point_vector will be the output.

Basic interpolation and raster output with a long spline step

v.surf.bspline input=point_vector raster=interpolate_surface sie=25 sin=25

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.

Estimation of lambda_i parameter with a cross validation proccess

v.surf.bspline -c input=point_vector

Estimation on sparse points

v.surf.bspline input=point_vector sparse=sparse_points output=interpolate_surface

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.

Using attribute values instead Z-coordinates

v.surf.bspline input=point_vector raster=interpolate_surface layer=1 column=attrib_column

This last case, the module uses the attribute values in attrib_column in the table associated to layer 1.

BUGS

Known issues:

In order to avoid RAM memory problems, an auxiliar table will be needed for recording some intermediate calculi. Since the "GROUP BY" SQL function is used, which is not supported by the "dbf" driver, this driver is not allowed with the vector map output "output=". There is no problem with the raster map output.

At this time, using the external vector input ("sparse=") implies interpoling with Z-coordinates. Updates to allow using attribute values will be done in a near future (I hope).

AUTHORS

Original version in GRASS 5.4: (s.bspline.reg)
Maria Antonia Brovelli, Massimiliano Cannata, Ulisse Longoni, Mirko Reguzzoni

Update for GRASS 6.X and improvements:
Roberto Antolin

REFERENCES

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

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

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)

Brovelli M. A., Cannata M. and Longoni U.M., 2002, DTM LIDAR in area urbana, Bollettino SIFET N.2, 2002, pp. 7-26

Last changed: $Date: 2007-02-26 03:08:52 -0800 (Mon, 26 Feb 2007) $

Main index - vector index - Full index