# Algorithm

Using semi-variograms instead of covariances is only justified under the unbiasedness conditions for the weights: S w i = 1 and S h j = 0 (see below), and with the Ordinary Kriging method (no trend, second order stationarity) (Deutsch & Journel 1992).

The direct use of variogram values (given the semi-variogram models g A , g B and cross-variogram model g AB), leads in case of m observations of predictand Ai and n observations of covariable Bj to the following system of CoKriging equations:

 (1)

If we call map A the collection of points with predictand values Ai, and map B the points with covariable values Bj and if A È B is the 'glued' combination of A and B it means that:

 GAA has m x m values of gA(h) from the lag vectors hij found in map A, h = || hij || GBA = G'AB containing m x n values of gAB(h) from the vectors hij found in map A È B. More precisely: the vectors defined by the set A x B (product set consisting of ordered pairs of points taken from A and B). GBB has n x n values of gB(h) from the lag vectors hij found in map B, h = || hij || 1m and 1'm are a column and a row vector both of order m; all elements in the vector are equal to 1 1n and 1'n are a column and a row vector both of order n; all elements in the vector are equal to 1 w is a column with m weights w i assigned to predictand A, and S w i = 1 h is a column with n weights h i assigned to covariable B, and S h j = 0 g pA is a column with m semi-variogram values g A (h), where h = || h pi ||, the distance from the estimated output pixel p to all visited (sampled) points in map A (predictand map). g pAB is a column with n cross-variogram values g AB (h), where h = || h pi ||, the distance from the estimated output pixel p to all visited (sampled) points in map B (covariable map). m1 is the Lagrange parameter introduced to formalize the unbiasedness and used to find the error variance of the prediction. m2 is the Lagrange parameter for the covariable.

The solution of the above system gives optimal values for w, h and m1 and m2.

These solutions (weights) do not depend on the sampled values in map A or B, but solely on the variogram models used and on the geometric distribution of the measurements, the so-called sampling scheme. They are used to get the prediction (Formula 2) as a linear combination of predictand measurements Ai and covariable measurements and Bj:

 (2)

The variance of the prediction error (s2 , Formula 3) is obtained from inner vector products and m1:

 s2 = S wi gA(hi) + S hj g AB(hj) + m1 (3)

This expression depends solely on the variogram models g A , g B and g AB and on the geometric layout of the sampled points (the sampling scheme) in both map A and map B.

Note: When the spherical distance option is used, the length of all lag vectors (h) are calculated over the sphere using the projection of the coordinate system that is used by the georeference of the output raster map.

References:

• Deutsch, C.V., and A.G. Journel. 1992. Geostatistical software library and user's guide. Oxford University Press, New York. 340 pp.
• Isaaks, E. H., and R. M. Srivastava. 1989. An introduction to applied Geostatistics. Oxford University Press, New York. 561 pp.
• Meer, F. D. van der. 1993. Introduction to geostatistics. ITC Lecture Notes. 72 pp.
• Pebesma, E.J. 1996. Mapping Groundwater Quality in the Netherlands. KNAG/Faculteit Ruimtelijke Wetenschappen, Utrecht University. 128 pp.
• Stein, A. and Corsten, L.C.A. 1991. Universal Kriging and CoKriging as a Regression Procedure. Biometrics June 1991. 13 pp.
• Stein, A. 1998. Spatial statistics for soils and the environment. ITC lecture notes. 47 pp.