Ordinary Kriging Example: R via text file

Submitted by dylan on Mon, 2007-06-11 01:58.

 
Overview:
A simple example of how to use R to perform interpolation with ordinary kriging, using data from a text file. This example makes use of the gstat library for R. Additional examples of how to use the following gstat functions are included:

  • variogram maps
  • directional variogram plots
  • ploting the interpolated surface directly from R

Note that this example is not meant to be an authoritative guide on variogram selection, or proper modeling of anisotropy-- just an example. The Kansas Geological Survey has an interesting set of reports that illustrate selection of a directional variogram in the presence of a strong, regional trend.

Elevation Data and Sample Points: 300 randomly placed points where elevation data was sampled.Original elevation data and sample points: 300 randomly placed points where elevation data was sampled.

 
Data Prep:
Export the coordinates and elevation values from the previous example. See attached file elev.txt.

# two new columns to the random point vector from the previous example
v.db.addcol map=rs columns="x double, y double"
# upload coordinates
v.to.db option=coor column=x,y map=rs
# export to text file
db.select rs fs="," > elev.csv

 
Start R:
Load in the text file, and coerce to format that gstat can use.

## load some libraries first:
library(gstat)
## load data
d <- read.csv('elev.csv')

## gstat does not like missing data, subset original data:
e <- na.omit(d)

## convert simple data frame into a spatial data frame object:
coordinates(e) <- ~ x+y

## test result with simple bubble plot:
bubble(e, zcol='elev', fill=FALSE, do.sqrt=FALSE, maxsize=2)

## create a grid onto which we will interpolate:
## first get the range in data
x.range <- as.integer(range(e@coords[,1]))
y.range <- as.integer(range(e@coords[,2]))

## now expand to a grid with 500 meter spacing:
grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=500), y=seq(from=y.range[1], to=y.range[2], by=500) )

## convert to SpatialPixel class
coordinates(grd) <- ~ x+y
gridded(grd) <- TRUE

## test it out:
plot(grd, cex=0.5)
points(e, pch=1, col='red', cex=0.7)
title("Interpolation Grid and Sample Points")

Interpolation GridInterpolation Grid

 
Create GSTAT Objects:
Make some diagnostic plots, model variogram, check for anisotropy, etc.

## make gstat object:
g <- gstat(id="elev", formula=elev ~ 1, data=e)

## the original data had a large north-south trend, check with a variogram map
plot(variogram(g, map=TRUE, cutoff=4000, width=200), threshold=10)

## another approach:
# create directional variograms at 0, 45, 90, 135 degrees from north (y-axis)
v <- variogram(g, alpha=c(0,45,90,135))

## 0 and 45 deg. look good. lets fit a linear variogram model:
## an un-bounded variogram suggests additional source of anisotropy... oh well.
v.fit <- fit.variogram(v, model=vgm(model='Lin' , anis=c(0, 0.5)))

## plot results:
plot(v, model=v.fit, as.table=TRUE)

## update the gstat object:
g <- gstat(g, id="elev", model=v.fit )

Variogram MapVariogram Map

Directional Variogram PlotsDirectional Variogram Plots

 
Perform OK and View Results:
Examples using standard and lattice graphics.

## perform ordinary kriging prediction:
p <- predict(g, model=v.fit, newdata=grd)

## visualize it:

## base graphics
par(mar=c(2,2,2,2))
image(p, col=terrain.colors(20))
contour(p, add=TRUE, drawlabels=FALSE, col='brown')
points(e, pch=4, cex=0.5)
title('OK Prediction')

## lattice graphics: thanks for R. Bivand's advice on this
##
## alternatively plot quantiles with
## ... col.regions=terrain.colors(6), cuts=quantile(p$elev.pred) ...
##
pts <- list("sp.points", e, pch = 4, col = "black", cex=0.5)
spplot(p, zcol="elev.pred", col.regions=terrain.colors(20), cuts=19, sp.layout=list(pts), contour=TRUE, labels=FALSE, pretty=TRUE, col='brown', main='OK Prediction')

## plot the kriging variance as well
spplot(p, zcol='elev.var', col.regions=heat.colors(100), cuts=99, main='OK Variance',sp.layout=list(pts) )

## quit and convert saved EPS files to PNG:
## for i in *.eps ; do convert $i `basename $i .eps`.png ; done

OK Prediction: created with the spplot() functionOK Prediction: created with the spplot() function

OK VarianceOK Variance

AttachmentSize
elev.csv_.txt10.96 KB