R universal kriging with autoKrige()

Question

I'm trying to use the autoKrige() function in the automap package for a simple application of universal kriging. I have an irregularly spaced grid of measurements, and I want to interpolate between them on a fine spatial scale. Example code:

    library('automap')

    # create an irregularly spaced grid
    y <-x <-c(-5,-4,-2,-1,-0.5,0,0.5,1,2,4,5)
    grid <-expand.grid(x,y)
    names(grid) <-c('x', 'y')

    # create some measurements, greatest in the centre, with some noise
    vals <-apply(grid,1, function(x) {12/(0.1+sqrt(x[1]^2 + x[2]^2))+rnorm(1,2,1.5)})

    # get data into sp format
    s <-SpatialPointsDataFrame(grid, data.frame(vals))

    # make some prediction locations and get them into sp format
    pred <-expand.grid(seq(-5,5,by=0.5), seq(-5,5,by=0.5))
    pred <-cbind(pred[,1], pred[,2])    # this seems to be needed, not sure why
    pred <-SpatialPoints(pred)

    # try universal kriging
    surf <-autoKrige(vals~x+y, s, new_data=pred)

This results in the error:

    Error in gstat.formula.predict(d$formula, newdata, na.action = na.action,  : 
      NROW(locs) != NROW(X): this should not occur

I have tried making new_data have the same number of rows as the original data, and have even tried making the co-ordinates in new_data exactly the same as the original data, but I still get this error. I'm new to geostatistics techniques so apologies if I'm making a basic mistake. Can anyone advise where I'm going wrong? Thanks.

The `pred <-cbind(pred[,1], pred[,2]) ` step can be left out for me, `autoKrige` works even without it. — Paul Hiemstra, Sep 08 '15 at 14:03

score 6 · Accepted Answer · answered Sep 08 '15 at 13:57

The problem is that you have the syntax of the autoKrige function wrong. The formula input to autoKrige specifies the linear model you want to use, e.g.:

log(zinc) ~ dist

from the meuse dataset. In this case, you model log(zinc) versus dist using a linear model, and the residuals to this model are interpolated using the variogram. Essentially, universal kriging is linear regression with spatially correlated residuals.

In your case, you specify:

val ~ x+y

so autoKrige (gstat actually) will try to first model the linear model of vals versus x and y (multivariate regression), and interpolate the residuals using the variogram model. However, the x and y variables are not present in the SpatialPointsDataFrame.

What I think you want to do is to only interpolate spatially using the variogram model. In that case, the linear model is very simple, actually just fitting a mean value:

vals ~ 1

where the mean of vals is determined and the residuals are interpolated using the variogram model. This is actually known as Ordinary Kriging. Your call to autoKrige would be something like:

surf <-autoKrige(vals ~ 1, s, new_data=pred)

R universal kriging with autoKrige()

1 Answers1