CRAN package gstat comes with functions for spatio-temporal variogram modelling and prediction. The vignettes on this topic might be a good starting point.
demo(localKrigeST)
demonstrates ST kriging using a local neighbourhood.
You will find all you need in the excellent (and didactic) technical note from Rossiter (2012)*:
Technical Note: Co-kriging with the gstat package of the R environment for statistical computing.
Co-kriging will use different functions from those with univariate kriging (for example, ordinary kriging).
The datasets (target and co-variables) should remain in separate data frames, but within the same object of class gstat. And predictions (interpolation) are carried out with predict.gstat.
In Rossiter (2012), chapters 6 and 7 explain in details how to do it:
- (6) Modelling a bivariate co-regionalisation.
- (7) Co-kriging with one co-variable.
Below is the main code from Rossiter (2012) which addresses the question. It uses the meuse dataset as example:
g <- gstat(NULL, id = "ltpb", form = ltpb ~ 1, data=meuse.pb) #target variable lead (pb).
g <- gstat(g, id = "ltom", form = ltom ~ 1, data=meuse.co) #co-variable organic matter (om).
v.cross <- variogram(g) #generate direct variograms and the cross-variogram.
g <- gstat(g, id = "ltpb", model = m.ltpb.f, fill.all=T) #add variogram models to gstat object. In this case, it has been used the variogram model. fitted to the target variable in previous chapter, for both target variable and the co-variable as starting points.
g <- fit.lmc(v.cross, g) #fit theoretical variograms to experimental ones (uses linear model of co-regionalisation).
k.c <- predict.gstat(g, meuse.grid) #predicts values for target variable in the prediction grid.
*Rossiter, D.G. 2012. Technical Note: Co-kriging with the gstat package of the R environment for statistical computing. University of Twente, Faculty of Geo-Information Science & Earth. Observation (ITC). Enschede (NL). Revision 2.3. 84p.
Best Answer
In
sp,SpatialPoints*,SpatialPixels*andSpatialGrid*(with*omitted or replaced byDataFrame) do support more than 2 spatial dimensions, as OP has done, butSpatialPolygons*andSpatialLines*do not. Withgstatyou can do 3-D block kriging with 3-D blocks (usingblock = c(10,10,10)), but you cannot do this for non-rectangular blocks, as OP wants. It is perfectly OK to substitute time for the third dimension, but you are constrained to the metric ST variogram.gives you more options for variogram models, but not for predicting block mean values (this is FYI, not an answer to the question).
The only answer to the question would be to do 3D conditional simulations, and aggregate point values over your arbitrary 3D (2D polygon + time extent) blocks. Tedious, but possible; also only along the 3D path, not along the path described in the ST vignette (
krigeSTdoes not do simulation - yet!).