I have a very sparse set of measurements of a quantity, around 10 measurements at different locations across a continent. The distributions of locations is highly irregular. I would like to smooth these into values on a regularly spaced grid over the continent.

We actually have a much denser set of observations of a related quantity and techniques are being developed to create a map of the first quantity using both sets of measurements. The intention is to estimate the accuracy of this hybrid approach using leave one out cross validation (testing each of the 10 sites in turn having created the map using all the others plus the additional measurements). As a baseline, I first need to see how accurate a map is, using LOOCV, created just from the original network of 10 stations. I should also point out this quantity varies over timescale of minutes and we want to produce a map in near real time.

I am not an expert in geospatial techniques, however I comprehend the basics of Kriging and after checking the task view of CRAN, it seems autoKrige in the automap package would be a suitable function to employ. The results I am getting are nothing like what I expected so any help would be appreciated.

Here is an example of one set of input training data (all available data for one particular time):

   lat   lon      value
  -27.53 152.92   98
  -35.32 149.00   79
  -34.05 150.67   81
  -12.45 130.95   92
  -42.92 147.32   73
  -22.25 114.08   91
  -29.03 167.97  108
  -31.94 115.95   89

As you can see, the target quantity has a range of values over a widely spaced area.

I convert from a data frame containing lat, lon and the value into a spatialPointsDataFrame by doing

coordinates(xtrain) = ~lon + lat
coordinates(xtest) = ~lon + lat
proj4string(xtrain) <-CRS("")
proj4string(xtest) <-CRS("")

where xtest contains the grid of points I want to sample the map over. I then Krige this using

result <- autoKrige(value~1,xtrain,xtest) 

I plot the results


and get the following image

Kriging result

As you can see, the result is essentially the mean value at all points. My understanding is that autoKrige should have determined the best values of the Kriging parameters and produced something more realistic than this. Just to illustrate, the range of the predicted values on the grid is 0.22, which is clearly much less than the 30 point scatter seen in the input values.

So, where am I going wrong? Am I misusing the R package, or is my problem more fundamental in the way I am attempting to employ Kriging?

  • 1
    $\begingroup$ Attempting to krige a dataset of just $10$ points is a lost cause unless you have strong prior information about the variogram. (Rules of thumb in the literature range from $30$ to $100$ points as reasonably minima.) Using an automatic kriging routine for any dataset is just foolhardy. Nevertheless, this one seems to have done a good job: the variogram is reasonable. As a result, the map is a complete artifact of whatever neighborhood search method the program chose. One moral to draw from this experience is that you tend to get out of statistical procedures what you put into them. $\endgroup$ – whuber Jul 25 '12 at 15:53

From ten points you are going to have 45 (10*(10-1)/2) points in your variogram cloud from the distances between each pair of points. Once the system has binned that, or even without binning, its going to be dominated by noise, I reckon. Get a plot of the variogram cloud to see what I mean.

If autokrige can't fit a nice smooth variogram then it will do what it did, and just go 'heck, I can't work out the correlation with distance with just 10 points, my best guess is just the mean'. It really can't do better.

If you want something to look 'realistic', then you could feed it variogram parameters with a bigger range, that would over-smooth the output. But then you may as well just do inverse-distance weighting if all you want is a pretty picture. The advantage of kriging is that it is realistic. But it rejects your reality and replaces it with its own...


  • Get a plot of the variogram cloud
  • Get more data :)
  • Look into bivariate kriging for your case with the two different data sets. I think the theory exists, there may even be code for it...
  • 2
    $\begingroup$ (+1) "Bivariate kriging" is called "cokriging." There's lots of commercial code. For an open source solution geoRglm will do it, I believe. $\endgroup$ – whuber Jul 25 '12 at 15:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.