Search code examples
rgeospatialspatialkriginggstat

Universal kriging using lat long gstat R


I'm new at R and I'm having some trouble to perform a universal kriging with gstat R.

As Hengl et al. (2004) say "Universal kriging should be reserved for the case where the drift (or trend) is modelled as a function of the coordinates only". So, I want to use only the coordinates and not the dist in universal kriging.

Anyone can tell me how? I'm proceeding like this:

library(sp)
library(gstat)

data(meuse)
coordinates(meuse) <- c("x", "y")
data(meuse.grid)
coordinates(meuse.grid) <- c("x", "y")
meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")

plot(variogram(log(zinc) ~ meuse@coords, meuse),pch=19)
v1 <- variogram(log(zinc) ~ meuse@coords, meuse)
p1 <- vgm(psill = 0.42, model = 'Cir', range = 1000, nugget = 0.08)
fit1 <- fit.variogram(v1, p1)

# Trying to use the coordinates with meuse@coords
uk1 <- krige(log(zinc) ~ meuse@coords, meuse, meuse.grid, fit1)

# Trying to get coordinates as data column
xy <- as.data.frame(meuse@coords)
meuse$long <- xy$x
meuse$lat <- xy$y

uk2 <- krige(log(zinc) ~ meuse$long + meuse$lat, meuse, meuse.grid, fit1)

Thank you!


Solution

  • I am actually not quite sure what T. Hengl meant in his 2004 paper. That paper is about the distinction between UK, RK and KED, which are actually all equivalent. First you fit a trend, a multiple linear model, and then a variogram model is fitted on those residuals.

    If you want to use the coordinates as drift in your UK model in R you can simply do:

    library(sp)
    library(gstat)
    
    data(meuse)
    coordinates(meuse) <- c("x", "y")
    data(meuse.grid)
    coordinates(meuse.grid) <- c("x", "y")
    meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")
    
    meuse$X <- as.numeric(meuse@coords[,1])
    meuse$Y <- as.numeric(meuse@coords[,2])
    
    meuse.grid$X <-as.numeric(meuse.grid@coords[,1])
    meuse.grid$Y <-as.numeric(meuse.grid@coords[,2])
    
    plot(variogram(log(zinc) ~ X+Y, meuse),pch=19)
    v1 <- variogram(log(zinc) ~  X+Y, meuse)
    p1 <- vgm(psill = 0.42, model = 'Cir', range = 1000, nugget = 0.08)
    fit1 <- fit.variogram(v1, p1)
    
    
    
    # Trying to use the coordinates with meuse@coords
    uk1 <- krige(log(zinc) ~ X+Y, meuse, meuse.grid, fit1)
    
    # Trying to get coordinates as data column
    uk2 <- krige(log(zinc) ~  X+Y , meuse, meuse.grid, fit1)