# Ordinary Kriging Example: GRASS-R Bindings

**Update: 2012-02-13**

Many of the examples used in this demonstration are now somewhat dated, probably inefficient, and in need of revision. I'll spend some time on an updated version for the GRASS wiki soon.

**Overview:**

A simple example of how to use GRASS+R to perform interpolation with ordinary kriging, using data from the spearfish sample dataset. This example makes use of the gstat library for R.

**Helpful References:**

- Issaks, E.H. & Srivastava, R.M. An Introduction to Applied Geostatistics Oxford University Press, 1989
- GSTAT Manual
- GSTAT Examples

**Data Prep:**

As a contrived example, we will generate 300 random points within the current region, and sample an elevation raster at each of these points.

**Load GRASS Data into R:**

Remember that you will need to install these R packages onto your computer.

**Variogram Modeling:**

A very simple example, using default parameters for a *non-directional variogram* is presented below. Modeling the variogram for an actual spatial problem requires knowlege of both your dataset (distribution, collection methods, etc.), the natural processes involved (stationarity vs. anisotropy ?), and a bit about the assumptions of geostatistics.

**Interpolation by Ordinary Kriging:**

The prediction is done for every instance of a '1' in the object passed to the `newdata=` argument.

**Send Results Back to GRASS:**

**Viewing Results in GRASS:**

**Simple Comparison with RST:**

RST (regularized splines with tension) and OK (ordinary kriging) are two common interpolation methods. Computing the RMSE (root-mean-square-error) between the interpolated raster and the original raster provides a simple quantitative measure of how well the interpolation performed, at least in terms mean *magnitude* of error. A spatial description of interpolation error can be generated by subtracting the new raster from the original. Note that the steps involve cell-wise computation of the square-error (SE), region-wise computation of the mean-square-error (MSE); the square root of MSE gives the root-mean-square-error or RMSE.

**Root-mean-square-error Comparison:**

Looks like the RSME are pretty close...

Method | OK | RST |

RMSE | 61 meters | 64 meters |

## Software

- General Purpose Programming with Scripting Languages
- LaTeX Tips and Tricks
- PostGIS: Spatially enabled Relational Database Sytem
- PROJ: forward and reverse geographic projections
- GDAL and OGR: geodata conversion and re-projection tools
- R: advanced statistical package
- Access Data Stored in a Postgresql Database
- Additive Time Series Decomposition in R: Soil Moisture and Temperature Data
- Aggregating SSURGO Data in R
- Cluster Analysis 1: finding groups in a randomly generated 2-dimensional dataset
- Color Functions
- Comparison of Slope and Intercept Terms for Multi-Level Model
- Comparison of Slope and Intercept Terms for Multi-Level Model II: Using Contrasts
- Creating a Custom Panel Function (R - Lattice Graphics)
- Customized Scatterplot Ideas
- Estimating Missing Data with aregImpute() {R}
- Exploration of Multivariate Data
- Interactive 3D plots with the rgl package
- Making Soil Property vs. Depth Plots
- Numerical Integration/Differentiation in R: FTIR Spectra
- Plotting XRD (X-Ray Diffraction) Data
- Using lm() and predict() to apply a standard curve to Analytical Data
- Working with Spatial Data
- Customizing Maps in R: spplot() and latticeExtra functions
- Converting Alpha-Shapes into SP Objects
- Some Ideas on Interpolation of Categorical Data
- Visual Interpretation of Principal Coordinates (of) Neighbor Matrices (PCNM)
- Visualizing Random Fields and Select Components of Spatial Autocorrelation
- Generation of Sample Site Locations [sp package for R]
- Target Practice and Spatial Point Process Models
- Ordinary Kriging Example: GRASS-R Bindings
- Point-process modelling with the sp and spatstat packages
- Simple Map Creation

- Comparison of PSA Results: Pipette vs. Laser Granulometer

- GRASS GIS: raster, vector, and imagery analysis
- Generic Mapping Tools: high quality map production