Digging Kriging
Or… Kriging with GRASS, R and OS X
Given a set of points with attributes, how do you go about interpolating between those points to create a surface? The three most common methods of doing this are Inverse Distance Weighting (IDW), Regularized Spline with Tension (RST) and Kriging. GRASS GIS can natively create IDW and RST-based interpolated surfaces, but needs a little help from its friend R to do Kriging. Here’s a quick guide to setting up R with the associated packages on OSX so that you can create Krigged surfaces. This assumes that you’re already running GRASS 6 (if not, grab a copy of OpenOSX GRASS from http://www.openosx.com/grass/freedownload.html)
Krigged surface, rendered with nviz
1. Install MacPorts: Download and install MacPorts, MacPorts is an easy to use system for building linux/unix libraries and apps on OSX
2. Install R, proj, gdal: Use MacPorts to install the statistical computing application ‘R’, cartographic projection library proj and geospatial data abstraction library gdal. Open up a terminal and type this in:
sudo port install R proj gdal
You’ll be asked for your password. Once you’ve done that MacPorts will download the source, compile and install R, proj and gdal. NB This takes a long while (> 1 hour on a MacBook 2.0GHz) so go outside and enjoy the sun, watch a film, go for a walk etc…
Install R packages R requires three additional packages installed before we get going. First off we need to install the rgdal package to get access to the gdal and proj libraries. Start up R by typing r at a terminal and hit enter. Then enter the following to download, compile and install the rgdal package.
install.packages(‘rgdal’, configure.args=’—with-proj-lib=/opt/usr/local/lib’)
Then enter the following to install the spgrass6 which gives R direct access to GRASS6 data:
install.packages(‘spgrass6’)
That’s it
To Krig a surface you’ll need to start GRASS, set the current region to the one that you wish to create a surface within (tip: use integer values for the boundaries of the region), then in a GRASS termnal start R by entering ‘r’
You can then adapt the following script, to create a krigged surface. This script is a modified version of the one found here. This script interpolates the points in the vector map name ‘points_to_krig’ with values in the field ‘dbl_1’. The kriging is carried out using a default linear variogram, and output to the GRASS raster k_surf. Save the script below to a file (eg /tmp/krig.txt) and run using r’s source() function eg source(’/tmp/krig.txt’)
##load the required libraries
library(gstat)
library(spgrass6)
##retrieve the GRASS6 environment
G <- gmeta6()
##load points
source_points <- readVECT6('points_to_krig')
#remove any duplicate points
my_points <- remove.duplicates(source_points)
## create a grid wihch will define the interpolation
grd <- gmeta2grd()
## make a new grid of (1)s
## be sure to use original data's proj data...
## doesn't work with the stuff stored in G...
new_grid <- SpatialGridDataFrame(grd, data=data.frame(k=rep(1,G$cols*G$rows)), proj4string=CRS(my_points@proj4string@projargs))
## init our gstat object, with the model formula
## note that coordinates are auto-identified from the GRASS object
g <- gstat(id="elev", formula=dbl_1 ~ 1, data=my_points)
## default linear variogram model
v.fit <- fit.variogram(variogram(g) ,vgm(model="Lin") )
## update gstat object with fitted variogram model
g <- gstat(g, id="elev", model=v.fit )
## interpolate with ord. kriging
k_surf <- predict(g, id='elev', newdata=new_data)
## write raster back to GRASS: interpolation and kriging variance:
## system('g.remove rast=elev.ok')
writeRAST6(k_surf, 'k_surf', zcol='elev.pred')
writeRAST6(k_surf, 'k_surf.var', zcol='elev.var')
#quit r
q()