Geografic-Information-System/raster/r.resamp.bspline/r.resamp.bspline.html

<h2>DESCRIPTION</h2>

<em>r.resamp.bspline</em> performs a bilinear/bicubic spline interpolation with
Tykhonov regularization. The input is a raster surface map, e.g. elevation,
temperature, precipitation etc. Output is a raster map. Optionally, only
input NULL cells are interpolated, useful to fill NULL cells, an alternative
to <em><a href="r.fillnulls.html">r.fillnulls</a></em>. Using the <b>-n</b> flag to only
interpolate NULL cells will considerably speed up the module.
<p>
The input raster map is read at its native resolution, the output raster
map will be produced for the current computational region set with
<em><a href="g.region.html">g.region</a></em>. Any MASK will be respected, masked
values will be treated as NULL cells in both the input and the output map.
<p>Spline step values <b>ew_step</b> for the east-west direction and
<b>ns_step</b> for the north-south direction should not be smaller than
the east-west and north-south resolutions of the input map. For a raster
map without NULL cells, 1 * resolution can be used, but check for
undershoots and overshoots. For very large areas with missing values
(NULL cells), larger spline step values may be required, but most of the
time the defaults (1.5 x resolution) should be fine.
<p>
The Tykhonov regularization parameter (<b>lambda</b>) acts to
smooth the interpolation. With a small <b>lambda</b>, the
interpolated surface closely follows observation points; a larger value
will produce a smoother interpolation. Reasonable values are 0.0001,
0.001, 0.005, 0.01, 0.02, 0.05, 0.1 (needs more testing). For seamless
NULL cell interpolation, a small value is required and default is set to 0.005.
<p>
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;
these are 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>ew_step</b> for the east-west direction and
<b>ns_step</b> for the north-south direction. For optimal performance, the
spline step values should be no less than the east-west and north-south
resolutions of the input map. Each non-NULL cell 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>A cross validation "leave-one-out" analysis is available to help to determine
the optimal <b>lambda</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. Note that cross
validation can run very slowly if more than 100 observations are used. 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>lambda</b> values. No vector nor raster output will be created
when cross-validation is selected. 

<h2>EXAMPLES</h2>

<h3>Basic interpolation</h3>

<div class="code"><pre>
r.resamp.bspline input=raster_surface output=interpolated_surface method=bicubic
</pre></div>

A bicubic spline interpolation will be done and a raster map with estimated
(i.e., interpolated) values will be created. 

<h3>Interpolation of NULL cells and patching</h3>

General procedure:
<div class="code"><pre>
# set region to area with NULL cells, align region to input map
g.region n=north s=south e=east w=west align=input -p
# interpolate NULL cells
r.resamp.bspline -n input=input_raster output=interpolated_nulls method=bicubic
# set region to area with NULL cells, align region to input map
g.region raster=input -p
# patch original map and interpolated NULLs
r.patch input=input_raster,interpolated_nulls output=input_raster_gapfilled
</pre></div>

<h3>Interpolation of NULL cells and patching (NC data)</h3>

In this example, the SRTM elevation map in the
North Carolina sample dataset location is filtered for outlier
elevation values; missing pixels are then re-interpolated to obtain
a complete elevation map:

<div class="code"><pre>
g.region raster=elev_srtm_30m -p
d.mon wx0
d.histogram elev_srtm_30m

r.univar -e elev_srtm_30m

# remove too low elevations (esp. lakes)
# Threshold: thresh = Q1 - 1.5 * (Q3 - Q1)
r.mapcalc "elev_srtm_30m_filt = if(elev_srtm_30m &lt; 50.0, null(), elev_srtm_30m)"

# verify
d.histogram elev_srtm_30m_filt
d.erase
d.rast elev_srtm_30m_filt

r.resamp.bspline -n input=elev_srtm_30m_filt output=elev_srtm_30m_complete \
  method=bicubic

d.histogram elev_srtm_30m_complete
d.rast elev_srtm_30m_complete
</pre></div>

<h3>Estimation of <b>lambda</b> parameter with a cross validation process</h3>

A random sample of points should be generated first with
<em><a href="r.random.html">r.random</a></em>, and the current region should not
include more than 100 non-NULL random cells. 

<div class="code"><pre>
r.resamp.bspline -c input=input_raster 
</pre></div>

<h2>REFERENCES</h2>

<ul>
<li>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</li>
<li>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</li>
<li>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)</li>
<li>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&ntilde;a. CD ISBN: 978-84-690-3886-9</li>
</ul>

<h2>SEE ALSO</h2>

<em>
<a href="r.fillnulls.html">r.fillnulls</a>,
<a href="r.resamp.rst.html">r.resamp.rst</a>,
<a href="r.resamp.interp.html">r.resamp.interp</a>,
<a href="v.surf.bspline.html">v.surf.bspline</a>
</em>

<p>
Overview: <a href="https://grasswiki.osgeo.org/wiki/Interpolation">Interpolation and Resampling</a> in GRASS GIS

<h2>AUTHORS</h2>
Markus Metz<br>
<br>
based on <em><a href="v.surf.bspline.html">v.surf.bspline</a></em> by
<br>
Maria Antonia Brovelli, Massimiliano Cannata, Ulisse Longoni, Mirko Reguzzoni, Roberto Antolin

<p>
<i>Last changed: $Date$</i>