I am getting acquainted with gstat
and kriging more generally in R. The end goal is to interpolate temperature across my spatial grid, and I want to be able to do so for different time points (average daily temperature over several years). So, I've been imagining an interpolation that estimates temperature and considers something like the season or year as a random effect. I've been exploring universal kriging, but I'm struggling to understand the formula set up.
I've included some example code using the meuse
data (from sp
). When I assign a variogram, I can state that my dependent variable is predicted by other variables and include random effects (rather than using something like zinc~1
).
However, when I do this with krige()
, the only alternative formula to zinc~1
seems to be zinc~x+y
. I believe the formula zinc~x+y
is what allows anisotropy and avoids the stationarity assumption, but why can't I include other data here like I can with the variogram? Is it because it is already accounted for in the variogram, or am I not understanding something about kriging more generally?
library(sp)
library(gstat)
data(meuse)
coordinates(meuse) <- ~ x+y
data(meuse.grid)
coordinates(meuse.grid) <- ~ x+y
# assigning and fitting the variogram
sample.vgm <- variogram(zinc~copper+(1|dist.m), # the formula - copper is a fixed effect, dist.m is a random effect
meuse) # SPDF dataset
vgm.fit <- fit.variogram(sample.vgm, # the sample variogram
model=vgm(model="Sph")) # assigning a model type
# kriging
z.krige <- krige(zinc~x+y, # the formula
meuse, # the SPDF that has data
meuse.grid, # the SPDF to krige over
model=vgm.fit)
# this doesn't work if the formula is anything besides ~1 or ~x+y
z.krige2 <- krige(zinc~copper+(1|dist.m),
meuse,
meuse.grid,
model=vgm.fit)
I am using R version 4.1.3 on Windows 10 (x64).
Based on this post, the explanatory variables apparently need to be loaded into the grid for Kriging to find them.
dist
should work since it is already in the meuse
grid:
# dist is in the grid
z.krige2 <- krige(zinc~dist,
meuse,
meuse.grid,
model=vgm.fit)
For copper
and dist.m
, I've just added random values here:
# add random values for copper and dist.m into the grid
meuse.grid$copper <- runif(nrow(meuse.grid))
meuse.grid$dist.m <- runif(nrow(meuse.grid))
# now copper and dist.m variables are available
z.krige2 <- krige(zinc~copper+(1|dist.m),
meuse,
meuse.grid,
model=vgm.fit)