![GRASS logo](grass_logo.png)
NAME
v.krige.py
DESCRIPTION
v.krige allows to perform kriging operations in GRASS
environment, using R software functions in background.
NOTES
v.krige is just a front-end to R. The options and parameters
are the same offered by packages automap and gstat.
Kriging, like other interpolation methods, is fully dependent on input
data features. Exploratory analysis of data is encouraged to find out
outliers, trends, anisotropies, uneven distributions and consequently
choose the kriging algorithm that will give the most acceptable
result. Good knowledge of the dataset is more valuable than hundreds
of parameters or powerful hardware. See Isaaks and Srivastava's book,
exhaustive and clear even if a bit outdated.
Dependencies
- R software >= 2.x
- rpy2
- Python binding to R. Note! rpy version 1 is not supported.
- R packages automap and gstat.
- automap is optional (provides automatic variogram fit).
Notes for Debian GNU/Linux
Install the dependiencies. Attention! python-rpy IS NOT
SUITABLE.:
aptitude install R python-rpy2
To install R packages, use either R's function (as root):
install.packages("gstat", dep=T)
install.packages("spgrass6", dep=T)
either the brand new Debian packages [5], add to repositories' list
for 32bit or 64bit (pick up the suitable line):
and get the packages via aptitude:
aptitude install r-cran-gstat r-cran-spgrass6
Notes for Windows
Compile GRASS following this
guide.
You could also use Linux in a virtual machine. Or install Linux in a
separate partition of the HD. This is not as painful as it appears,
there are lots of guides over the Internet to help you.
Computation time issues
Please note that kriging calculation is slown down both by high number
of input data points and/or high region resolution, even if they both
contribute to a better output.
EXAMPLES
Kriging example based on elevation map (Spearfish data set).
Part 1: random sampling of 2000 vector points from known
elevation map. Each point will receive the elevation value from the
elevation raster, as if it came from a point survey.
g.region rast=elevation.10m -p
v.random output=rand2k_elev n=2000
v.db.addtable map=rand2k_elev column="elevation double precision"
v.what.rast vect=rand2k_elev rast=elevation.10m column=elevation
Part 2: remove points lacking elevation attributes. Points
sampled at the border of the elevation map didn't receive any
value. v.krige has no preferred action to cope with no data values, so
the user must check for them and decide what to do (remove points,
fill with the value of the nearest point, fill with the global/local
mean...). In the following line of code, points with no data are
removed from the map.
v.extract rand2k_elev output=rand2k_elev_filt where="elevation not NULL"
Check the result of previous line ("number of NULL attributes" must be
0):
v.univar rand2k_elev_filt type=point column=elevation
Part 3: reconstruct DEM through kriging. Using automatic
variogram fit is the simplest way to run v.krige from CLI (note:
requires R's automap package). Output map name is optional, the
modules creates it automatically appending "_kriging" the the input
map name and also checks for overwrite. If output_var is specified,
the variance map is also created. Automatic variogram fit is provided
by R package automap, the variogram models tested by the fitting
functions are: exponential, spherical, Gaussian, Matern, M.Stein's
parametrisation. A wider range of models is available from gstat
package and can be tested on the GUI via the variogram plotting. If
model is specified in the CLI, also sill, nugget and range values are
to be provided, otherwise an error is raised (see second example of
v.krige command).
v.krige input=rand2k_elev_filt column=elevation output=rand2k_elev_kriging \
output_var=rand2k_elev_kriging_var
v.krige input=rand2k_elev_filt column=elevation output=rand2k_elev_kriging \
output_var=rand2k_elev_kriging_var model=Lin sill=2500 nugget=0 range=1000 \
--overwrite
Or run wxGUI, to interactively fit the variogram and explore options:
Calculate prediction error:
r.mapcalc "rand2k_elev_kriging_pe = sqrt(rand2k_elev_kriging_var)"
r.univar elevation.10m
r.univar rand2k_elev_kriging
r.univar rand2k_elev_kriging_pe
The results show high errors, as the kriging techniques (ordinary and
block kriging) are unable to handle a dataset with a trend, like the
one used in this example: elevation is higher in the southwest corner
and lower on northeast corner. Universal kriging can give far better
results in these cases as it can handle the trend. It is available in
R package gstat and will be part of a future v.krige release.
SEE ALSO
R package gstat,
mantained by Edzer J. Pebesma and others
R
package spgrass6,
mantained by Roger Bivand
The Short
Introduction to Geostatistical and Spatial Data Analysis with GRASS 6
and R statistical data language at the GRASS website. (includes
installation tips)
v.krige's wiki page
REFERENCES
Isaaks and Srivastava, 1989: "An Introduction to Applied Geostatistics"
(ISBN 0-19-505013-4)
AUTHOR
Anne Ghisla, Google Summer of Code 2009
Last changed: $Date: 2011-11-08 13:56:45 -0800 (Tue, 08 Nov 2011) $
Main index - vector index - Full index
© 2003-2014 GRASS Development Team