Understanding Collocated CoKriging¶

This test is meant to explain the Collocated CoKriging principle and to answer to several questions posted a while ago.

In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import numpy as np

gdoc.setNoScroll()

The principle of Collocated CoKriging¶

The Collocated CoKriging stands in the multivariate case (nvar stands for the number of variables. The Kriging (or CoKriging) estimation consists in interpolating the variables at some target sites, starting from information known at data sites. Obviously the Model must be defined consistently.

Notes:

  • all the information may not be kwnown at all the data sites: we say that we are in an heterotopic case
  • some target sites may coincide with data sites: we will call them coincidence case the variables can be considered as stationary (we must specify their means and we may practice Simple Kriging) or intrinsic (we may practice Ordinary Kriging). Higher drift orders are possible but they will not be envisaged in this document.

The principle of Colocated CoKriging is to add to the Data information (of the Data file), some information available in the Target File: we must simply define the connection between the ofvariables in both files. This is done using the rankColCok parameter.

It consists of a vector of nvar ranks, such that rankColCok[ivar] = jvar means that:

  • for the variable 'ivar' (in the Data file and in the Model)
  • the corresponding information is available in the Target file as variable jvar
  • if no matching variable is available, simply set jvar = -1

This method clearly requires that the Target file already contains some variables defined.

Data¶

This tutorial is built as an demonstration task: i.e. we will consider very reduced number of data and target sites to better understand the functionality.

We consider:

  • a case with 3 variables
  • an input Data file with 3 samples, located randomly
  • a Model constituted of a single exponential structure
  • an output Target file with 3 sites located randomly for demonstration sake, the coordinates of the second target site are updated so as to match the ones of the third data site
  • we finally consider using either a Moving or a Unique neighborhood.
In [2]:
ndim = 2
ndat = 3
nvar = 3
nout = 3

# Global parameters
gl.law_set_random_seed(32131)
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)

We build the Models: a stationary Model (called modelSK) and an Intrinsic Model (called modelOK)

In [3]:
types = [gl.ECov.EXPONENTIAL]

modelSK = gl.Model.createFillRandom(ndim, nvar, types, 1.0, -1)
means = gl.VectorHelper.simulateGaussian(nvar)
modelSK.setMeans(means)

modelOK = gl.Model.createFillRandom(ndim, nvar, types, 1.0, 0)

We build the Neighborhoods: a Unique neighborhood (called neighU) and a large Moving neighborhood (called neighM)

In [4]:
neighU = gl.NeighUnique.create()

neighM = gl.NeighMoving.create(False, ndat, 5.0)

We build the input Data file (called data). We transform the information into an heterotopic one by setting some values to NA.

In [5]:
data = gl.Db.createFillRandom(ndat, ndim, nvar, seed=243415)
data.setLocVariable(gl.ELoc.Z, 1, 0, gl.TEST)

We build the output Data file (called targetRef). We copy the coordinates of the third data site (called matchDat) into the coordinates of the second target site (called matchOut).

In [6]:
target = gl.Db.createFillRandom(nout, ndim, nvar, seed=534243)

matchDat = 2
matchOut = 1
target.setCoordinate(matchOut, 0, data.getCoordinate(matchDat, 0))
target.setCoordinate(matchOut, 1, data.getCoordinate(matchDat, 1))

We finally parametrize the Collocation Option:

In [7]:
rank_colcok = [1, -1, 0]
krigoptColCok = gl.KrigOpt()
err = krigoptColCok.setColCok(rank_colcok)

krigopt = gl.KrigOpt()

Let us have an exhaustive printout of the information in both Data and Target files.

In [8]:
# Data file
data.getContentsAsTable().display()

        rank         x-1         x-2         z-1         z-2         z-3
       1.000       0.278       0.508       0.156       0.968       0.028
       2.000       0.181       0.380         N/A       0.069       1.427
       3.000       0.052       0.914       0.751       1.303       0.163
In [9]:
# Target file
target.getContentsAsTable().display()

        rank         x-1         x-2         z-1         z-2         z-3
       1.000       0.805       0.715      -0.953       0.938      -0.656
       2.000       0.052       0.914       2.280       1.416       1.511
       3.000       0.407       0.139      -0.621      -1.376      -0.402

We can check that the sample 'matchDat' of the Data file coincides with the sample 'matchOut' of the Target file.

Problems¶

We now run several cases of Collocated CoKriging and display the results for comparison.

Simple CoKriging (no Collocation Option) in Moving Neighborhood¶

In [10]:
target.deleteColumns(["Kriging*"])
target.setLocators(["z*"], gl.ELoc.Z)
gl.kriging(data, target, modelSK, neighM, True, True, True)
target.getContentsAsTable(["Kriging*"]).display()

  *z-1.estim  *z-2.estim  *z-3.estim  *z-1.stdev  *z-2.stdev  *z-3.stdev  *.z-1.varz  *.z-2.varz  *.z-3.varz
       0.610       0.078       1.059       2.272       1.956       2.470       0.006       0.005       0.007
       0.751       1.303       0.163       0.000       0.000       0.000       5.167       3.831       6.110
       0.697       0.087       1.087       2.249       1.936       2.445       0.110       0.082       0.131

We can check that the variances are systematically zero for the second target site, as a data coincides. The resulting values coincide with the data values. At other target sites, there is nothing special.

Simple CoKriging (no Collocation Option) in Unique Neighborhood¶

In [11]:
target.deleteColumns(["Kriging*"])
target.setLocators(["z*"], gl.ELoc.Z)
gl.kriging(data, target, modelSK, neighU, True, True, True)
target.getContentsAsTable(["Kriging*"]).display()

  *z-1.estim  *z-2.estim  *z-3.estim  *z-1.stdev  *z-2.stdev  *z-3.stdev  *.z-1.varz  *.z-2.varz  *.z-3.varz
       0.610       0.078       1.059       2.272       1.956       2.470       0.006       0.005       0.007
       0.751       1.303       0.163       0.000       0.000       0.000       5.167       3.831       6.110
       0.697       0.087       1.087       2.249       1.936       2.445       0.110       0.082       0.131

We obtain the same results as before which prooves that the two sets of equations (whether a Moving or a Unique neighborhood is used) lead to the same results.

Ordinary CoKriging (no Collocation Option) in Moving Neighborhood¶

In [12]:
target.deleteColumns(["Kriging*"])
target.setLocators(["z*"], gl.ELoc.Z)
gl.kriging(data, target, modelOK, neighM, True, True, True)
target.getContentsAsTable(["Kriging*"]).display()

  *z-1.estim  *z-2.estim  *z-3.estim  *z-1.stdev  *z-2.stdev  *z-3.stdev  *.z-1.varz  *.z-2.varz  *.z-3.varz
       1.124       0.835       0.495       2.686       2.313       2.921       2.238       1.659       2.646
       0.751       1.303       0.163       0.000       0.000       0.000       5.167       3.831       6.110
       1.284       0.744       0.597       2.570       2.213       2.794       2.283       1.692       2.699

The results are comparable to those obtained using Simple CoKriging. Note the variances (when not zero) are slightly larger as expected.

Ordinary CoKriging (no Collocation Option) in Unique Neighborhood¶

In [13]:
target.deleteColumns(["Kriging*"])
target.setLocators(["z*"], gl.ELoc.Z)
gl.kriging(data, target, modelOK, neighU, True, True, True)
target.getContentsAsTable(["Kriging*"]).display()

  *z-1.estim  *z-2.estim  *z-3.estim  *z-1.stdev  *z-2.stdev  *z-3.stdev  *.z-1.varz  *.z-2.varz  *.z-3.varz
       1.124       0.835       0.495       2.686       2.313       2.921       2.238       1.659       2.646
       0.751       1.303       0.163       0.000       0.000       0.000       5.167       3.831       6.110
       1.284       0.744       0.597       2.570       2.213       2.794       2.283       1.692       2.699

Once more, the results are the same for Moving and Unique neighborhoods.

Simple Collocated CoKriging in Moving Neighborhood¶

In [14]:
target.deleteColumns(["Kriging*"])
target.setLocators(["z*"], gl.ELoc.Z)
gl.kriging(data, target, modelSK, neighM, True, True, True, krigopt=krigoptColCok)
target.getContentsAsTable(["Kriging*"]).display()

  *z-1.estim  *z-2.estim  *z-3.estim  *z-1.stdev  *z-2.stdev  *z-3.stdev  *.z-1.varz  *.z-2.varz  *.z-3.varz
       0.938      -1.428      -0.953       0.000       0.088       0.000       5.167       3.823       6.110
       0.751       1.303       0.163       0.000       0.000       0.000       5.167       3.831       6.110
      -1.376       0.693      -0.621       0.000       0.087       0.000       5.167       3.823       6.110

Here is the interesting point of the Collocated CoKriging option:

At any target, the estimation of variable #1 is equal to the value of the variable #2 of the Target File; the estimation of variable #3 is equal to the value of the variable #1 of the Target File. This corresponds to the rule given in the rankColCok vector.

As an example, let us consider the first target site:

  • the value of variable #1 is -0.953
  • the value of variable #2 is 0.938

In the same file:

  • the estimated value of the variable #1 is 0.938
  • the estimated value of the variable #3 is -0.953

Another interesting point is the case of Target site #2 which coincides with Data site #3. There are two possible options for the results of variable 1 (for example):

  • the estimated value corresponds to the variable #1 of Data site #3 (i.e. 0.751) as Kriging is an exact interpolator
  • the estimated value corresponds to the variable #3 of Target site #2 (i.e. 1.511) carried by the Targer file

At the same Target site #2 and for variable #2, the options are:

  • the estimated value corresponds to the variable #2 of Data site #3 (i.e. 1.303) as Kriging is an exact interpolator
  • the estimated value corresponds to the variable #1 of Target site #2 (i.e. 2.200) carried by the Targer file

The first option is used. No specific message is provided even if the two possible values are not equal.

Simple Collocated CoKriging in Unique Neighborhood¶

In [15]:
target.deleteColumns(["Kriging*"])
target.setLocators(["z*"], gl.ELoc.Z)
gl.kriging(data, target, modelSK, neighU, True, True, True, krigopt=krigoptColCok)
target.getContentsAsTable(["Kriging*"]).display()

  *z-1.estim  *z-2.estim  *z-3.estim  *z-1.stdev  *z-2.stdev  *z-3.stdev  *.z-1.varz  *.z-2.varz  *.z-3.varz
       0.938      -1.428      -0.953       0.000       0.088       0.000       5.167       3.823       6.110
       0.751       1.303       0.163       0.000       0.000       0.000       5.167       3.831       6.110
      -1.376       0.693      -0.621       0.000       0.088       0.000       5.167       3.823       6.110

Again, the results are the same for Moving and Unique Neighborhoods.

Ordinary Collocated CoKriging in Moving Neighborhood¶

In [16]:
target.deleteColumns(["Kriging*"])
target.setLocators(["z*"], gl.ELoc.Z)
gl.kriging(data, target, modelOK, neighM, True, True, True, krigopt=krigoptColCok)
target.getContentsAsTable(["Kriging*"]).display()

  *z-1.estim  *z-2.estim  *z-3.estim  *z-1.stdev  *z-2.stdev  *z-3.stdev  *.z-1.varz  *.z-2.varz  *.z-3.varz
       0.938       0.089      -0.953       0.000       0.107       0.000       5.167       3.826       6.110
       0.751       1.303       0.163       0.000       0.000       0.000       5.167       3.831       6.110
      -1.376       2.122      -0.621       0.000       0.105       0.000       5.167       3.826       6.110

The results for Ordinary CoKriging (even with the Collocated option) gives larger variances (when positive) for Ordinary than for Simple case. The difference is very light as, in both cases, there is some collocated information which reduces the possible uncertainty.

Ordinary Collocated CoKriging in Unique Neighborhood¶

In [17]:
target.deleteColumns(["Kriging*"])
target.setLocators(["z*"], gl.ELoc.Z)
gl.kriging(data, target, modelOK, neighU, True, True, True, krigopt=krigoptColCok)
target.getContentsAsTable(["Kriging*"]).display()

  *z-1.estim  *z-2.estim  *z-3.estim  *z-1.stdev  *z-2.stdev  *z-3.stdev  *.z-1.varz  *.z-2.varz  *.z-3.varz
       0.938       0.089      -0.953       0.000       0.107       0.000       5.167       3.826       6.110
       0.751       1.303       0.163       0.000       0.000       0.000       5.167       3.831       6.110
      -1.376       2.122      -0.621       0.000       0.105       0.000       5.167       3.825       6.110

Same as for Moving Neighborhood.