interpolation - Preparing dataset to perform co-kriging in R gstat?

This is the first time I'm using co-kriging in gstat. My problem is that I'm not sure how to prepare the data frame to supply to co-kriging when the variable of interest and auxiliary variables are not measured at the same locations.

Data frame with the variable of interest (z_ar) is voi_gs

head(voi_gs)
             x        y     z_ar
8974  312216.6 530439.8 49.03470
8283  312084.6 530559.8 57.15355
5057  311772.6 530883.8 49.13453

1551  311976.6 531207.8 44.73679
10116 312168.6 530163.8 50.45549
6155  312528.6 530787.8 62.70750

Data frame with the auxiliary variable (z_pts) is aux_gs

head(aux_gs)
            x        y    z_pts
9564 312198.6 530313.8 75.91368
6584 311766.6 530745.8 38.87462
2138 311562.6 531153.8 58.62555

9110 312534.6 530421.8 68.35654
9525 312366.6 530325.8 54.26653
7497 311442.6 530649.8 38.95024

I combined them into one dataframe to supply to variogram() and krige() functions. Since none of the locations are same between voi_gs and aux_gs, I introduced NA values and combined them in the following way:

aux_gs$z_ar=NA
voi_gs$z_pts=NA
comb_gs=rbind(aux_gs,voi_gs)
head(comb_gs)


            x        y    z_pts z_ar
9564 312198.6 530313.8 75.91368   NA
6584 311766.6 530745.8 38.87462   NA
2138 311562.6 531153.8 58.62555   NA
9110 312534.6 530421.8 68.35654   NA
9525 312366.6 530325.8 54.26653   NA
7497 311442.6 530649.8 38.95024   NA

tail(comb_gs)
            x        y z_pts     z_ar

180  312468.6 531363.8    NA 70.54528
8633 312264.6 530511.8    NA 44.34631
7694 312492.6 530631.8    NA 57.30173
1079 312108.6 531255.8    NA 46.96482
2230 311124.6 531135.8    NA 40.36449
2201 312312.6 531147.8    NA 44.85896

and then, constructed cross variogram:

coordinates(comb_gs) = ~x+y
g = gstat(formula=z_pts~z_ar, data=comb_gs)

vg = variogram(g)

but variogram function does not accept NA values and gives me an error. I know I can't have NA values in the data, but I don't know how else I can create the data frame to construct the cross-variogram.

Answer

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.}