(This article was first published on Bart Rogiers , and kindly contributed toR-bloggers)

Last year I wrote a short demo on variography with gstat and ggplot2 for a colleague who was planning to migrate to R . Just thought I’d share this here (with some additional stuff) as it might be useful for other people as well.

First, make sure you have the necessary packages installed and loaded:

require('gstat')

## Loading required package: gstat

require('sp')

## Loading required package: sp

require('ggplot2')

## Loading required package: ggplot2

The sp package is not really required, but it makes working with spatial data a little bit easier.

I like to work with the black & white ggplot2 theme, and legends positioned above the graphs, so I always adjust the basic settings the following way:

theme_set(theme_bw())

theme_update(legend.position='top')

Let’s create some synthetic data to work with. I’m using gstat here to simulate a random field on a 100×100 regular grid:

This is the data we will perform some variography on. Alternatively, have a look at this file to load your own data. We’ll transform our data frame into a spatial points data frame, and create a gstat object:

I commented the above line, as there is an issue with gstat 1.1-0, which basically mirrors gridded data horizontally before deriving the experimental variogram. Edzer Pebesma, the author of gstat, already solved the issue , so in the latest gstat releases this should work properly with gridded data as well.

g<-gstat(id='z',formula=z~1,data=dat)

Calculating an experimental isotropic variogram can then simply be done by:

expvar<-variogram(g)head(expvar)# show first lines of the gstatVariogram data frame

Finally, my colleague was interested in the performance of the algorithm. So here are some execution times for my personal laptop, running Linux Lite 2.6 , R 3.2.3, and having an Intel(R) Core(TM)2 Duo CPU T5800 @ 2.00GHz processor and 2969MB memory:

system.time(variogram(g))

## user system elapsed ## 5.192 0.072 7.698

system.time(variogram(g,alpha=0,tol.hor=5))

## user system elapsed ## 4.828 0.104 6.553

system.time(variogram(g,alpha=c(0,45,90,135),tol.hor=5))

## user system elapsed ## 19.240 0.440 26.115

system.time(variogram(g,alpha=c(0,45,90,135),tol.hor=5,cutoff=100))

## user system elapsed ## 29.040 0.680 40.121

That’s it! If you want to know more, have a look at

?variogram

for more options, and if you want to try this yourself, check out the source R script !