rkriging

Kriging at a single spatial point?

This question may have to do with my poor knowledge of kriging - is it possible to compute kriged value at a single spatial location? As I understand it, a typical kriging method uses the spatial correlation embedded in a sparse data to interpolate values at all points on a regular grid that has the same spatial extent as the data. What I'd like to know is if I can do this interpolation on a specific point (that may not fall on the grid). As an example, the code below does kriging on copper concentration in the meuse data and overlays the values on the Meuse map. I'd like to know how to compute kriged value at "Maasband", near the center of the map.

Thanks.

# transform meuse data to SpatialPointsDataFrame
suppressMessages(library(sp))
data(meuse)
coordinates(meuse) <- ~ x + y
proj4string(meuse) <- CRS("+proj=stere
                          +lat_0=52.15616055555555 +lon_0=5.38763888888889
                          +k=0.999908 +x_0=155000 +y_0=463000
                          +ellps=bessel +units=m +no_defs
                          +towgs84=565.2369,50.0087,465.658,
                          -0.406857330322398,0.350732676542563,-1.8703473836068, 4.0812")

# define a regular grid for kriging
xrange <- range(as.integer(meuse@coords[, 1])) + c(0,1)
yrange <- range(as.integer(meuse@coords[, 2]))
grid <- expand.grid(x = seq(xrange[1], xrange[2], by = 40),
                    y = seq(yrange[1], yrange[2], by = 40))
coordinates(grid) <- ~ x + y
gridded(grid) <- T

# do kriging
suppressMessages(library(automap))
krg <- autoKrige(formula = copper ~ 1,
                 input_data = meuse,
                 new_data = grid)

# extract kriged data
krg_df <- data.frame(krg$krige_output@coords,
                     pred = krg$krige_output@data$var1.pred)

# transform to SpatialPointsDF & assign original (meuse) projection
krg_spdf <- krg_df
coordinates(krg_spdf) <- ~ x + y 
proj4string(krg_spdf) <- proj4string(meuse)

# transform again to longlat coordinates (for overlaying on google map below)
krg_spdf <- spTransform(krg_spdf, CRS("+init=epsg:4326"))
krg_df <- data.frame(krg_spdf@coords, pred = krg_spdf@data$pred)

# get meuse map and overlay kriged data
suppressMessages(library(ggmap))
suppressMessages(library(RColorBrewer))
lon <- range(krg_df$x)
lat <- range(krg_df$y)
meuse_map <- get_map(location = c(lon = mean(lon), lat = mean(lat)),
                     zoom = 13)
print(ggmap(meuse_map, extent = "normal", maprange = F) +
        stat_summary_2d(aes(x = x, y = y, z = pred),
                        binwidth = c(0.001,0.001),
                        alpha = 0.5,
                        data = krg_df) +
        scale_fill_gradientn(name = "Copper",
                             colours = brewer.pal(6, "YlOrRd")) +
        coord_cartesian(xlim = lon, ylim = lat, expand = 0) +
        theme(aspect.ratio = 1))

# geocode for Maasband
longlat <- as.numeric(geocode("Maasband"))
x0 <- longlat[1]
y0 <- longlat[2]

enter image description here

Solution

  • You can predict a single point and the point can be outside the training grid, but statistically it is not correct and if you get too far away you will only get a average.

    Example of prediction of a point outside training (newdata):

    > range(coordinates(meuse)[,1])
[1] 178605 181390
> range(coordinates(meuse)[,2])
[1] 329714 333611
> newdata <- data.frame(x = 178600,
+                     y = 329710)
> coordinates(newdata) <- ~ x + y
> krg <- autoKrige(formula = copper ~ 1,
+                  input_data = meuse,
+                  new_data = newdata)
[using ordinary kriging]
> krg
$krige_output
       coordinates var1.pred var1.var var1.stdev
1 (178600, 329710)  43.10343 538.8824   23.21384

$exp_var
     np       dist    gamma dir.hor dir.ver   id
1    17   59.33470 142.2353       0       0 var1
2    36   86.01449 288.4722       0       0 var1
3   114  131.02870 259.1184       0       0 var1
4   149  176.18845 417.4295       0       0 var1
5   184  226.75652 322.5353       0       0 var1
6   711  337.60359 473.3298       0       0 var1
7   830  502.04773 529.3259       0       0 var1
8  1349  713.21485 594.9974       0       0 var1
9  1314  961.27179 628.1739       0       0 var1
10 1139 1213.41157 612.3446       0       0 var1
11 1355 1506.55052 583.7055       0       0 var1

$var_model
  model     psill    range kappa
1   Nug  29.45239   0.0000   0.0
2   Ste 592.45663 365.6081   0.4

$sserr
[1] 78.06413

attr(,"class")
[1] "autoKrige" "list"