Standard 2-D case study¶

Import packages¶

In [1]:
import numpy as np
import pandas as pd
import sys
import os
import matplotlib.pyplot as plt
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc

Global variables

In [2]:
verbose  = True
graphics = True
gl.OptCst.define(gl.ECst.NTCOL,6)
gdoc.setNoScroll()

Reading data¶

The data are stored in a CSV format in the file called Pollution.dat

In [3]:
filepath = gdoc.loadData("Pollution", "Pollution.dat")
mydb = gl.Db.createFromCSV(filepath,gl.CSVformat())
mydb.setLocators(["X","Y"],gl.ELoc.X)
mydb.setLocator("Zn",gl.ELoc.Z)
if verbose:
    dbfmt = gl.DbStringFormat()
    dbfmt.setFlags(flag_resume=True, flag_vars=True, flag_extend=True) 
    mydb.display(dbfmt)
Data Base Characteristics
=========================

Data Base Summary
-----------------
File is organized as a set of isolated points
Space dimension              = 2
Number of Columns            = 5
Total number of samples      = 102

Data Base Extension
-------------------
Coor #1 - Min =    109.850 - Max =    143.010 - Ext = 33.16
Coor #2 - Min =    483.660 - Max =    513.040 - Ext = 29.38

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = X - Locator = x1
Column = 2 - Name = Y - Locator = x2
Column = 3 - Name = Zn - Locator = z1
Column = 4 - Name = Pb - Locator = p1
 

Accessing to the variable names

In [4]:
print("List of all variable names =",mydb.getAllNames())
List of all variable names = ('rank', 'X', 'Y', 'Zn', 'Pb')

Extracting the vector containing the Zn variable in order to perform a selection

In [5]:
tabZn = mydb.getColumn('Zn') # equivalent to mydb["Zn"] or mydb[3]
selZn = tabZn < 20
mydb.addSelection(selZn,'sel')
mydb.setLocator('Pb',gl.ELoc.Z)
dbfmt = gl.DbStringFormat.createFromFlags(flag_stats=True)
if verbose:
    mydb.display(dbfmt)
Data Base Characteristics
=========================

Data Base Summary
-----------------
File is organized as a set of isolated points
Space dimension              = 2
Number of Columns            = 6
Total number of samples      = 102
Number of active samples     = 99

Data Base Statistics
--------------------
1 - Name rank - Locator NA
 Nb of data          =        102
 Nb of active values =         99
 Minimum value       =      1.000
 Maximum value       =    102.000
 Mean value          =     51.808
 Standard Deviation  =     29.114
 Variance            =    847.650
2 - Name X - Locator x1
 Nb of data          =        102
 Nb of active values =         99
 Minimum value       =    109.850
 Maximum value       =    143.010
 Mean value          =    119.924
 Standard Deviation  =      6.549
 Variance            =     42.890
3 - Name Y - Locator x2
 Nb of data          =        102
 Nb of active values =         99
 Minimum value       =    483.660
 Maximum value       =    513.040
 Mean value          =    498.622
 Standard Deviation  =      8.000
 Variance            =     64.001
4 - Name Zn - Locator NA
 Nb of data          =        102
 Nb of active values =         99
 Minimum value       =      1.090
 Maximum value       =     12.100
 Mean value          =      2.881
 Standard Deviation  =      1.667
 Variance            =      2.779
5 - Name Pb - Locator z1
 Nb of data          =        102
 Nb of active values =         99
 Minimum value       =      3.000
 Maximum value       =     12.700
 Mean value          =      5.658
 Standard Deviation  =      1.697
 Variance            =      2.881
6 - Name sel - Locator sel
 Nb of data          =        102
 Nb of active values =         99
 Minimum value       =      1.000
 Maximum value       =      1.000
 Mean value          =      1.000
 Standard Deviation  =      0.000
 Variance            =      0.000

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = X - Locator = x1
Column = 2 - Name = Y - Locator = x2
Column = 3 - Name = Zn - Locator = NA
Column = 4 - Name = Pb - Locator = z1
Column = 5 - Name = sel - Locator = sel
 

Display my Data (with samples represented by color and size)

In [6]:
if graphics:
    ax = mydb.plot(nameColor="Pb")
    ax.decoration(title="Data Set")

Variograms¶

We first define the geometry of the variogram calculations

In [7]:
myVarioParamOmni = gl.VarioParam()
mydir = gl.DirParam.create(npas=10, dpas=1.)
myVarioParamOmni.addDir(mydir)

We use the variogram definition in order to calculate the variogram cloud.

In [8]:
dbcloud = gl.db_variogram_cloud(mydb, myVarioParamOmni)

We recall that the Variogram cloud is calculated by filling an underlying grid where each cell is painted according to the number of pairs at the given distance and given variability. Representing the variogram cloud

In [9]:
if graphics:
    ax = dbcloud.plot("Cloud*")
    ax.decoration(title="Variogram Cloud")

Calculating the experimental omni-directional variogram

In [10]:
myVarioOmni = gl.Vario(myVarioParamOmni)
err = myVarioOmni.compute(mydb, gl.ECalcVario.VARIOGRAM)
if verbose:
    myVarioOmni.display()
Variogram characteristics
=========================
Number of variable(s)       = 1
Number of direction(s)      = 1
Space dimension             = 2
Variance-Covariance Matrix     2.881

Direction #1
------------
Number of lags              = 10
Direction coefficients      =      1.000     0.000
Direction angles (degrees)  =      0.000     0.000
Tolerance on direction      =     90.000 (degrees)
Calculation lag             =      1.000
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0     3.000     0.389     0.462
         1   123.000     1.081     1.495
         2   183.000     2.038     1.620
         3   205.000     3.006     2.526
         4   231.000     4.013     2.240
         5   229.000     5.036     2.524
         6   198.000     5.962     2.396
         7   187.000     7.000     2.708
         8   204.000     7.996     2.772
         9   184.000     8.990     2.868
 

The variogram is represented graphically for a quick check

In [11]:
if graphics:
    ax = myVarioOmni.plot()
    ax.decoration(title="Omni-directional Variogram for Pb")

Calculate a variogram in several directions

In [12]:
myvarioParam = gl.VarioParam()
mydirs = gl.DirParam.createMultiple(4, 10, 1.)
myvarioParam.addMultiDirs(mydirs)
myvario = gl.Vario(myvarioParam)
myvario.compute(mydb,gl.ECalcVario.VARIOGRAM)
if verbose:
    myvario.display()
Variogram characteristics
=========================
Number of variable(s)       = 1
Number of direction(s)      = 4
Space dimension             = 2
Variance-Covariance Matrix     2.881

Direction #1
------------
Number of lags              = 10
Direction coefficients      =      1.000     0.000
Direction angles (degrees)  =      0.000     0.000
Tolerance on direction      =     22.500 (degrees)
Calculation lag             =      1.000
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0     1.000     0.410     0.180
         1    29.000     1.094     1.634
         2    47.000     2.079     1.415
         3    53.000     3.003     2.824
         4    63.000     3.999     2.348
         5    66.000     5.035     2.319
         6    60.000     5.978     3.115
         7    52.000     7.045     2.746
         8    52.000     8.020     3.927
         9    37.000     8.980     2.554

Direction #2
------------
Number of lags              = 10
Direction coefficients      =      0.707     0.707
Direction angles (degrees)  =     45.000     0.000
Tolerance on direction      =     22.500 (degrees)
Calculation lag             =      1.000
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0     1.000     0.344     0.080
         1    31.000     1.051     1.113
         2    50.000     1.960     1.890
         3    62.000     2.999     2.443
         4    58.000     4.014     2.701
         5    51.000     5.016     2.702
         6    36.000     5.999     1.833
         7    37.000     7.015     2.130
         8    50.000     7.997     2.060
         9    53.000     8.995     2.381

Direction #3
------------
Number of lags              = 10
Direction coefficients      =      0.000     1.000
Direction angles (degrees)  =     90.000     0.000
Tolerance on direction      =     22.500 (degrees)
Calculation lag             =      1.000
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         1    32.000     1.149     1.631
         2    39.000     2.080     1.670
         3    39.000     2.979     2.511
         4    48.000     4.012     2.120
         5    51.000     5.029     3.055
         6    47.000     5.939     2.856
         7    49.000     6.965     2.386
         8    42.000     7.952     2.708
         9    41.000     9.018     2.320

Direction #4
------------
Number of lags              = 10
Direction coefficients      =     -0.707     0.707
Direction angles (degrees)  =    135.000     0.000
Tolerance on direction      =     22.500 (degrees)
Calculation lag             =      1.000
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0     1.000     0.411     1.125
         1    31.000     1.028     1.606
         2    47.000     2.044     1.496
         3    51.000     3.040     2.330
         4    62.000     4.028     1.791
         5    61.000     5.058     2.155
         6    55.000     5.939     1.587
         7    49.000     6.975     3.425
         8    60.000     8.004     2.408
         9    53.000     8.972     3.996
 
In [13]:
if graphics:
    ax = myvario.plot(idir=-1)
    ax.decoration(title="Multi-Directional Variogram of Pb")

Calculating the Variogram Map

In [14]:
myvmap = gl.db_vmap_compute(mydb,gl.ECalcVario.VARIOGRAM,[20,20])
if verbose:
    myvmap.display()
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 5
Total number of samples      = 1681

Grid characteristics:
---------------------
Origin :    -33.160   -29.380
Mesh   :      1.658     1.469
Number :         41        41

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = VMAP.Pb.Var - Locator = z1
Column = 4 - Name = VMAP.Pb.Nb - Locator = NA
 
In [15]:
if graphics:
    ax = myvmap.plot("*Var")
    ax.decoration(title="Variogram Map")

Model¶

Fitting a Model. We call the Automatic Fitting procedure providing the list of covariance functions to be tested.

In [16]:
mymodel = gl.Model.createFromDb(mydb)
err = mymodel.fit(myvario,[gl.ECov.EXPONENTIAL,gl.ECov.SPHERICAL])

Visualizing the resulting model, overlaid on the experimental variogram

In [17]:
if graphics:
    ax = gp.varmod(myvario,mymodel)
    ax.decoration(title="Model for Pb")

Model with equality constraints¶

We can impose some constraints on the parameters during the fit. For instance here, we impose an equality constraint on the range (range = 1).

In [18]:
myModelConstrained = gl.Model.createFromDb(mydb)
constr = gl.Constraints()
paramid = gl.CovParamId(0,0,gl.EConsElem.RANGE,0,0)
constr.addItem(gl.ConsItem(paramid,gl.EConsType.EQUAL,1.))
err = myModelConstrained.fit(myVarioOmni,[gl.ECov.SPHERICAL],constr)
if (err > 0): print("Error while fitting the model")
myModelConstrained
Out[18]:
Model characteristics
=====================
Space dimension              = 2
Number of variable(s)        = 1
Number of basic structure(s) = 1
Number of drift function(s)  = 0
Number of drift equation(s)  = 0

Covariance Part
---------------
Spherical
- Sill         =      2.101
- Range        =      1.000
Total Sill     =      2.101

We can impose inequality constraints by using EConsType.LOWER or EConsType.UPPER.

Adding a drift¶

In [19]:
mymodel.setDriftIRF()
if verbose:
    mymodel.display()
Model characteristics
=====================
Space dimension              = 2
Number of variable(s)        = 1
Number of basic structure(s) = 2
Number of drift function(s)  = 1
Number of drift equation(s)  = 1

Covariance Part
---------------
Exponential
- Sill         =      1.035
- Ranges       =      1.786     0.366
- Theo. Ranges =      0.596     0.122
- Angles       =     45.023     0.000
- Rotation Matrix
               [,  0]    [,  1]
     [  0,]     0.707    -0.707
     [  1,]     0.707     0.707
Spherical
- Sill         =      1.621
- Ranges       =      7.051     5.132
- Angles       =    136.897     0.000
- Rotation Matrix
               [,  0]    [,  1]
     [  0,]    -0.730    -0.683
     [  1,]     0.683    -0.730
Total Sill     =      2.656

Drift Part
----------
Universality_Condition
 

Defining the Neighborhood¶

We initiate a Neigborhood (Moving with a small number of samples for Demonstration)

In [20]:
myneigh = gl.NeighMoving.create(flag_xvalid=False,nmaxi=6,radius=10)
if verbose:
    myneigh.display()
Moving Neighborhood
===================
Minimum number of samples           = 1
Maximum number of samples           = 6
Maximum horizontal distance         = 10
 

Checking the Moving Neighborhood¶

We must first create a Grid which covers the area of interest

In [21]:
mygrid = gl.DbGrid.createCoveringDb(mydb,[],[0.5,0.5],[],[2,2])
if verbose:
    mygrid.display()
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 2
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
 

We can now test the neighborhood characteristics for each node of the previously defined grid.

In [22]:
err = gl.test_neigh(mydb,mygrid,mymodel,myneigh)
if (err > 0): print("Error while running test_neigh")
if verbose:
    mygrid.display()
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 7
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
Column = 2 - Name = Neigh.Pb.Number - Locator = NA
Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA
Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA
Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA
Column = 6 - Name = Neigh.Pb.NbCESect - Locator = z1
 

We can visualize some of the newly created variables, such as:

  • the number of points per neighborhood
In [23]:
if graphics:
    ax = mygrid.plot("Neigh*Number")
    ax.decoration(title="Number of Samples per Neighborhood")
  • the one giving the maximum distance per neighborhood
In [24]:
if graphics:
    ax=mygrid.plot("Neigh*MaxDist")
    ax.decoration(title="Maximum Distance per Neighborhood")

Cross-validation¶

We can now process the cross-validation step

In [25]:
err = gl.xvalid(mydb,mymodel,myneigh)
if (err > 0): print("Error while running xvalid")
if verbose:
    mydb.display()
Data Base Characteristics
=========================

Data Base Summary
-----------------
File is organized as a set of isolated points
Space dimension              = 2
Number of Columns            = 8
Total number of samples      = 102
Number of active samples     = 99

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = X - Locator = x1
Column = 2 - Name = Y - Locator = x2
Column = 3 - Name = Zn - Locator = NA
Column = 4 - Name = Pb - Locator = NA
Column = 5 - Name = sel - Locator = sel
Column = 6 - Name = Xvalid.Pb.esterr - Locator = z1
Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA
 
In [26]:
if graphics:
    gp.histogram(mydb,name="Xvalid.Pb.stderr",bins=50,range=(-4,3))

Estimation by Kriging¶

We now perform the Estimation by Ordinary Kriging. The Neighborhood is changed into a Unique Neighborhood.

In [27]:
mydb.setLocator("Pb",gl.ELoc.Z)
myneigh = gl.NeighUnique.create()
err = gl.kriging(mydb,mygrid,mymodel,myneigh)
if (err > 0): print("Error while running kriging")
if verbose:
    mygrid.display(dbfmt)
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 9
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Data Base Statistics
--------------------
1 - Name x1 - Locator x1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    107.850
 Maximum value       =    144.850
 Mean value          =    126.350
 Standard Deviation  =     10.824
 Variance            =    117.167
2 - Name x2 - Locator x2
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    481.660
 Maximum value       =    515.160
 Mean value          =    498.410
 Standard Deviation  =      9.814
 Variance            =     96.313
3 - Name Neigh.Pb.Number - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      6.000
 Mean value          =      5.381
 Standard Deviation  =      1.538
 Variance            =      2.366
4 - Name Neigh.Pb.MaxDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.935
 Maximum value       =      9.999
 Mean value          =      5.598
 Standard Deviation  =      2.497
 Variance            =      6.233
5 - Name Neigh.Pb.MinDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.014
 Maximum value       =      9.978
 Mean value          =      3.484
 Standard Deviation  =      2.514
 Variance            =      6.320
6 - Name Neigh.Pb.NbNESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      1.000
 Mean value          =      1.000
 Standard Deviation  =      0.000
 Variance            =      0.000
7 - Name Neigh.Pb.NbCESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.000
 Maximum value       =      0.000
 Mean value          =      0.000
 Standard Deviation  =      0.000
 Variance            =      0.000
8 - Name Kriging.Pb.estim - Locator z1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.412
 Maximum value       =     11.456
 Mean value          =      6.123
 Standard Deviation  =      0.648
 Variance            =      0.419
9 - Name Kriging.Pb.stdev - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      0.248
 Maximum value       =      1.659
 Mean value          =      1.541
 Standard Deviation  =      0.159
 Variance            =      0.025

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
Column = 2 - Name = Neigh.Pb.Number - Locator = NA
Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA
Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA
Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA
Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA
Column = 7 - Name = Kriging.Pb.estim - Locator = z1
Column = 8 - Name = Kriging.Pb.stdev - Locator = NA
 

Visualizing the results

In [28]:
if graphics:
    ax = mygrid.plot("Kriging.Pb.estim")
    ax = mydb.plot("Pb")
    ax.decoration(title="Estimate of Pb")
In [29]:
if graphics:
    ax = mygrid.plot("Kriging.Pb.stdev")
    ax = mydb.plot("Pb")
    ax.decoration(title="St. Deviation of Pb")

Simulations¶

We must first transform the Data into Gaussian

In [30]:
myanamPb = gl.AnamHermite(30)
myanamPb.fitFromLocator(mydb)
if verbose:
    myanamPb.display()
Hermitian Anamorphosis
----------------------
Minimum absolute value for Y  = -2.7
Maximum absolute value for Y  = 2.6
Minimum absolute value for Z  = 3.0029
Maximum absolute value for Z  = 12.9777
Minimum practical value for Y = -2.7
Maximum practical value for Y = 2.6
Minimum practical value for Z = 3.0029
Maximum practical value for Z = 12.9777
Mean                          = 5.65758
Variance                      = 2.86296
Number of Hermite polynomials = 30
Normalized coefficients for Hermite polynomials (punctual variable)
               [,  0]    [,  1]    [,  2]    [,  3]    [,  4]    [,  5]    [,  6]
     [  0,]     5.658    -1.625     0.440    -0.069    -0.017     0.082    -0.061
     [  7,]     0.001     0.036    -0.044     0.004     0.047    -0.030    -0.029
     [ 14,]     0.037     0.007    -0.031     0.010     0.018    -0.019    -0.003
     [ 21,]     0.019    -0.010    -0.014     0.019     0.006    -0.023     0.004
     [ 28,]     0.022    -0.013
 

We can produce the Gaussian Anamorphosis graphically within its definition domain.

In [31]:
if graphics:
    gp.anam(myanamPb)

The next step consists in translating the target variable ('Pb') into its Gaussian transform

In [32]:
mydb.setLocator("Pb",gl.ELoc.Z)
err = myanamPb.rawToGaussianByLocator(mydb)
if verbose:
    mydb.display()
Data Base Characteristics
=========================

Data Base Summary
-----------------
File is organized as a set of isolated points
Space dimension              = 2
Number of Columns            = 9
Total number of samples      = 102
Number of active samples     = 99

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = X - Locator = x1
Column = 2 - Name = Y - Locator = x2
Column = 3 - Name = Zn - Locator = NA
Column = 4 - Name = Pb - Locator = NA
Column = 5 - Name = sel - Locator = sel
Column = 6 - Name = Xvalid.Pb.esterr - Locator = NA
Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA
Column = 8 - Name = Y.Pb - Locator = z1
 

We quickly calculate experimental (omni-directional) variograms of the gaussian variable (locater Z) using the already defined directions

In [33]:
myvarioParam = gl.VarioParam()
mydir = gl.DirParam(10,1.)
myvarioParam.addDir(mydir)
myVario = gl.Vario(myvarioParam)
err = myvario.compute(mydb,gl.ECalcVario.VARIOGRAM)

We fit the model by automatic fit (with the constraints that the total sill be equal to 1).

In [34]:
mymodelG = gl.Model.createFromDb(mydb)
err = mymodelG.fit(myvario,[gl.ECov.EXPONENTIAL])
if graphics:
    ax = gp.varmod(myvario,mymodelG)
    ax.decoration(title="Model for Gaussian Pb")

We perform a set of 10 conditional simulations using the Turning Bands Method.

In [35]:
err = gl.simtub(mydb,mygrid,mymodel,myneigh,10)
if verbose:
    mygrid.display()
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 19
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
Column = 2 - Name = Neigh.Pb.Number - Locator = NA
Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA
Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA
Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA
Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA
Column = 7 - Name = Kriging.Pb.estim - Locator = NA
Column = 8 - Name = Kriging.Pb.stdev - Locator = NA
Column = 9 - Name = Simu.Y.Pb.1 - Locator = z1
Column = 10 - Name = Simu.Y.Pb.2 - Locator = z2
Column = 11 - Name = Simu.Y.Pb.3 - Locator = z3
Column = 12 - Name = Simu.Y.Pb.4 - Locator = z4
Column = 13 - Name = Simu.Y.Pb.5 - Locator = z5
Column = 14 - Name = Simu.Y.Pb.6 - Locator = z6
Column = 15 - Name = Simu.Y.Pb.7 - Locator = z7
Column = 16 - Name = Simu.Y.Pb.8 - Locator = z8
Column = 17 - Name = Simu.Y.Pb.9 - Locator = z9
Column = 18 - Name = Simu.Y.Pb.10 - Locator = z10
 

Some statistics on the Conditional simulations in Gaussian scale

In [36]:
mygrid.deleteColumn("Stats.*")
mygrid.statisticsBySample(["Simu.Y.*"],[gl.EStatOption.MINI,
                                gl.EStatOption.MAXI,
                                gl.EStatOption.MEAN,
                                gl.EStatOption.STDV],
                                True)
if verbose:
    dbfmt = gl.DbStringFormat()
    dbfmt.setFlags(flag_resume=True, flag_stats=True) 
    mygrid.display(dbfmt)
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 23
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Data Base Statistics
--------------------
1 - Name x1 - Locator x1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    107.850
 Maximum value       =    144.850
 Mean value          =    126.350
 Standard Deviation  =     10.824
 Variance            =    117.167
2 - Name x2 - Locator x2
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    481.660
 Maximum value       =    515.160
 Mean value          =    498.410
 Standard Deviation  =      9.814
 Variance            =     96.313
3 - Name Neigh.Pb.Number - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      6.000
 Mean value          =      5.381
 Standard Deviation  =      1.538
 Variance            =      2.366
4 - Name Neigh.Pb.MaxDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.935
 Maximum value       =      9.999
 Mean value          =      5.598
 Standard Deviation  =      2.497
 Variance            =      6.233
5 - Name Neigh.Pb.MinDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.014
 Maximum value       =      9.978
 Mean value          =      3.484
 Standard Deviation  =      2.514
 Variance            =      6.320
6 - Name Neigh.Pb.NbNESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      1.000
 Mean value          =      1.000
 Standard Deviation  =      0.000
 Variance            =      0.000
7 - Name Neigh.Pb.NbCESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.000
 Maximum value       =      0.000
 Mean value          =      0.000
 Standard Deviation  =      0.000
 Variance            =      0.000
8 - Name Kriging.Pb.estim - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.412
 Maximum value       =     11.456
 Mean value          =      6.123
 Standard Deviation  =      0.648
 Variance            =      0.419
9 - Name Kriging.Pb.stdev - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      0.248
 Maximum value       =      1.659
 Mean value          =      1.541
 Standard Deviation  =      0.159
 Variance            =      0.025
10 - Name Simu.Y.Pb.1 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -5.862
 Maximum value       =      7.142
 Mean value          =      0.484
 Standard Deviation  =      1.573
 Variance            =      2.474
11 - Name Simu.Y.Pb.2 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -4.793
 Maximum value       =      5.789
 Mean value          =      0.142
 Standard Deviation  =      1.568
 Variance            =      2.459
12 - Name Simu.Y.Pb.3 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -4.882
 Maximum value       =      6.844
 Mean value          =      0.277
 Standard Deviation  =      1.603
 Variance            =      2.568
13 - Name Simu.Y.Pb.4 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -5.957
 Maximum value       =      5.443
 Mean value          =      0.117
 Standard Deviation  =      1.614
 Variance            =      2.604
14 - Name Simu.Y.Pb.5 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -5.515
 Maximum value       =      4.505
 Mean value          =      0.097
 Standard Deviation  =      1.432
 Variance            =      2.052
15 - Name Simu.Y.Pb.6 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -6.198
 Maximum value       =      5.125
 Mean value          =      0.214
 Standard Deviation  =      1.530
 Variance            =      2.340
16 - Name Simu.Y.Pb.7 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -5.859
 Maximum value       =      5.812
 Mean value          =      0.079
 Standard Deviation  =      1.560
 Variance            =      2.434
17 - Name Simu.Y.Pb.8 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -4.978
 Maximum value       =      5.598
 Mean value          =      0.271
 Standard Deviation  =      1.589
 Variance            =      2.526
18 - Name Simu.Y.Pb.9 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -5.340
 Maximum value       =      6.919
 Mean value          =      0.683
 Standard Deviation  =      1.603
 Variance            =      2.569
19 - Name Simu.Y.Pb.10 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -5.433
 Maximum value       =      5.800
 Mean value          =      0.157
 Standard Deviation  =      1.617
 Variance            =      2.614
20 - Name Stats.MINI - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -6.198
 Maximum value       =      1.423
 Mean value          =     -2.101
 Standard Deviation  =      0.964
 Variance            =      0.929
21 - Name Stats.MAXI - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -0.843
 Maximum value       =      7.142
 Mean value          =      2.624
 Standard Deviation  =      1.021
 Variance            =      1.043
22 - Name Stats.MEAN - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =     -1.983
 Maximum value       =      2.417
 Mean value          =      0.252
 Standard Deviation  =      0.594
 Variance            =      0.353
23 - Name Stats.STDV - Locator z1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      0.220
 Maximum value       =      2.788
 Mean value          =      1.419
 Standard Deviation  =      0.362
 Variance            =      0.131

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
Column = 2 - Name = Neigh.Pb.Number - Locator = NA
Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA
Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA
Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA
Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA
Column = 7 - Name = Kriging.Pb.estim - Locator = NA
Column = 8 - Name = Kriging.Pb.stdev - Locator = NA
Column = 9 - Name = Simu.Y.Pb.1 - Locator = NA
Column = 10 - Name = Simu.Y.Pb.2 - Locator = NA
Column = 11 - Name = Simu.Y.Pb.3 - Locator = NA
Column = 12 - Name = Simu.Y.Pb.4 - Locator = NA
Column = 13 - Name = Simu.Y.Pb.5 - Locator = NA
Column = 14 - Name = Simu.Y.Pb.6 - Locator = NA
Column = 15 - Name = Simu.Y.Pb.7 - Locator = NA
Column = 16 - Name = Simu.Y.Pb.8 - Locator = NA
Column = 17 - Name = Simu.Y.Pb.9 - Locator = NA
Column = 18 - Name = Simu.Y.Pb.10 - Locator = NA
Column = 19 - Name = Stats.MINI - Locator = NA
Column = 20 - Name = Stats.MAXI - Locator = NA
Column = 21 - Name = Stats.MEAN - Locator = NA
Column = 22 - Name = Stats.STDV - Locator = z1
 

We visualize a conditional simulation in Gaussian scale

In [37]:
if graphics:
    ax = mygrid.plot("Simu.Y.Pb.1")
    ax = mydb.plot("Pb")
    ax.decoration(title="One Simulation of Pb in Gaussian Scale")

We turn the Gaussian conditional simulations into Raw scale (using the Anamorphosis back transform) and get rid of the Gaussian conditional simulations.

In [38]:
myanamPb.gaussianToRaw(mygrid,"Simu.Y.*")
mygrid.deleteColumn("Simu.Y.*")
mygrid.deleteColumn("Stats.*")
if verbose:
    mygrid.display()
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 19
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
Column = 2 - Name = Neigh.Pb.Number - Locator = NA
Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA
Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA
Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA
Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA
Column = 7 - Name = Kriging.Pb.estim - Locator = NA
Column = 8 - Name = Kriging.Pb.stdev - Locator = NA
Column = 9 - Name = Z.Simu.Y.Pb.1 - Locator = z1
Column = 10 - Name = Z.Simu.Y.Pb.2 - Locator = z2
Column = 11 - Name = Z.Simu.Y.Pb.3 - Locator = z3
Column = 12 - Name = Z.Simu.Y.Pb.4 - Locator = z4
Column = 13 - Name = Z.Simu.Y.Pb.5 - Locator = z5
Column = 14 - Name = Z.Simu.Y.Pb.6 - Locator = z6
Column = 15 - Name = Z.Simu.Y.Pb.7 - Locator = z7
Column = 16 - Name = Z.Simu.Y.Pb.8 - Locator = z8
Column = 17 - Name = Z.Simu.Y.Pb.9 - Locator = z9
Column = 18 - Name = Z.Simu.Y.Pb.10 - Locator = z10
 

We calculate some statistics on the Conditional Simulations in Raw scale.

In [39]:
mygrid.deleteColumn("Stats.*")
mygrid.statisticsBySample(["Z.Simu.*"],[gl.EStatOption.MINI, 
                                gl.EStatOption.MAXI,
                                gl.EStatOption.MEAN,
                                gl.EStatOption.STDV],
                  True)
if verbose:
    dbfmt = gl.DbStringFormat()
    dbfmt.setFlags(flag_resume=True, flag_stats=True) 
    mygrid.display(dbfmt)
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 23
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Data Base Statistics
--------------------
1 - Name x1 - Locator x1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    107.850
 Maximum value       =    144.850
 Mean value          =    126.350
 Standard Deviation  =     10.824
 Variance            =    117.167
2 - Name x2 - Locator x2
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    481.660
 Maximum value       =    515.160
 Mean value          =    498.410
 Standard Deviation  =      9.814
 Variance            =     96.313
3 - Name Neigh.Pb.Number - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      6.000
 Mean value          =      5.381
 Standard Deviation  =      1.538
 Variance            =      2.366
4 - Name Neigh.Pb.MaxDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.935
 Maximum value       =      9.999
 Mean value          =      5.598
 Standard Deviation  =      2.497
 Variance            =      6.233
5 - Name Neigh.Pb.MinDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.014
 Maximum value       =      9.978
 Mean value          =      3.484
 Standard Deviation  =      2.514
 Variance            =      6.320
6 - Name Neigh.Pb.NbNESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      1.000
 Mean value          =      1.000
 Standard Deviation  =      0.000
 Variance            =      0.000
7 - Name Neigh.Pb.NbCESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.000
 Maximum value       =      0.000
 Mean value          =      0.000
 Standard Deviation  =      0.000
 Variance            =      0.000
8 - Name Kriging.Pb.estim - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.412
 Maximum value       =     11.456
 Mean value          =      6.123
 Standard Deviation  =      0.648
 Variance            =      0.419
9 - Name Kriging.Pb.stdev - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      0.248
 Maximum value       =      1.659
 Mean value          =      1.541
 Standard Deviation  =      0.159
 Variance            =      0.025
10 - Name Z.Simu.Y.Pb.1 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.871
 Standard Deviation  =      2.959
 Variance            =      8.757
11 - Name Z.Simu.Y.Pb.2 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.270
 Standard Deviation  =      2.785
 Variance            =      7.759
12 - Name Z.Simu.Y.Pb.3 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.532
 Standard Deviation  =      2.904
 Variance            =      8.434
13 - Name Z.Simu.Y.Pb.4 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.264
 Standard Deviation  =      2.835
 Variance            =      8.036
14 - Name Z.Simu.Y.Pb.5 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.125
 Standard Deviation  =      2.548
 Variance            =      6.492
15 - Name Z.Simu.Y.Pb.6 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.397
 Standard Deviation  =      2.731
 Variance            =      7.456
16 - Name Z.Simu.Y.Pb.7 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.174
 Standard Deviation  =      2.709
 Variance            =      7.340
17 - Name Z.Simu.Y.Pb.8 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.504
 Standard Deviation  =      2.871
 Variance            =      8.245
18 - Name Z.Simu.Y.Pb.9 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      7.261
 Standard Deviation  =      3.079
 Variance            =      9.481
19 - Name Z.Simu.Y.Pb.10 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.340
 Standard Deviation  =      2.852
 Variance            =      8.133
20 - Name Stats.MINI - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =      8.083
 Mean value          =      3.412
 Standard Deviation  =      0.570
 Variance            =      0.324
21 - Name Stats.MAXI - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      4.282
 Maximum value       =     12.978
 Mean value          =     11.377
 Standard Deviation  =      2.109
 Variance            =      4.446
22 - Name Stats.MEAN - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.438
 Maximum value       =     11.597
 Mean value          =      6.474
 Standard Deviation  =      1.063
 Variance            =      1.131
23 - Name Stats.STDV - Locator z1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      0.349
 Maximum value       =      4.357
 Mean value          =      2.539
 Standard Deviation  =      0.738
 Variance            =      0.544

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
Column = 2 - Name = Neigh.Pb.Number - Locator = NA
Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA
Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA
Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA
Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA
Column = 7 - Name = Kriging.Pb.estim - Locator = NA
Column = 8 - Name = Kriging.Pb.stdev - Locator = NA
Column = 9 - Name = Z.Simu.Y.Pb.1 - Locator = NA
Column = 10 - Name = Z.Simu.Y.Pb.2 - Locator = NA
Column = 11 - Name = Z.Simu.Y.Pb.3 - Locator = NA
Column = 12 - Name = Z.Simu.Y.Pb.4 - Locator = NA
Column = 13 - Name = Z.Simu.Y.Pb.5 - Locator = NA
Column = 14 - Name = Z.Simu.Y.Pb.6 - Locator = NA
Column = 15 - Name = Z.Simu.Y.Pb.7 - Locator = NA
Column = 16 - Name = Z.Simu.Y.Pb.8 - Locator = NA
Column = 17 - Name = Z.Simu.Y.Pb.9 - Locator = NA
Column = 18 - Name = Z.Simu.Y.Pb.10 - Locator = NA
Column = 19 - Name = Stats.MINI - Locator = NA
Column = 20 - Name = Stats.MAXI - Locator = NA
Column = 21 - Name = Stats.MEAN - Locator = NA
Column = 22 - Name = Stats.STDV - Locator = z1
 

We visualize a Conditional Simulation in Raw Scale

In [40]:
if graphics:
    ax = mygrid.plot("Z.Simu.Y.Pb.1")
    ax = mydb.plot("Pb")
    ax.decoration(title="One simulation of Pb in Raw Scale")

Let us now average the conditional simulations in order to have a comparison with the estimation by kriging.

In [41]:
mygrid.deleteColumn("Stats.*")
mygrid.statisticsBySample(["Z.Simu.*"],[gl.EStatOption.MEAN],True)
if verbose:
    dbfmt = gl.DbStringFormat()
    dbfmt.setFlags(flag_resume=True, flag_stats=True) 
    mygrid.display(dbfmt)
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 20
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Data Base Statistics
--------------------
1 - Name x1 - Locator x1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    107.850
 Maximum value       =    144.850
 Mean value          =    126.350
 Standard Deviation  =     10.824
 Variance            =    117.167
2 - Name x2 - Locator x2
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    481.660
 Maximum value       =    515.160
 Mean value          =    498.410
 Standard Deviation  =      9.814
 Variance            =     96.313
3 - Name Neigh.Pb.Number - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      6.000
 Mean value          =      5.381
 Standard Deviation  =      1.538
 Variance            =      2.366
4 - Name Neigh.Pb.MaxDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.935
 Maximum value       =      9.999
 Mean value          =      5.598
 Standard Deviation  =      2.497
 Variance            =      6.233
5 - Name Neigh.Pb.MinDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.014
 Maximum value       =      9.978
 Mean value          =      3.484
 Standard Deviation  =      2.514
 Variance            =      6.320
6 - Name Neigh.Pb.NbNESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      1.000
 Mean value          =      1.000
 Standard Deviation  =      0.000
 Variance            =      0.000
7 - Name Neigh.Pb.NbCESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.000
 Maximum value       =      0.000
 Mean value          =      0.000
 Standard Deviation  =      0.000
 Variance            =      0.000
8 - Name Kriging.Pb.estim - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.412
 Maximum value       =     11.456
 Mean value          =      6.123
 Standard Deviation  =      0.648
 Variance            =      0.419
9 - Name Kriging.Pb.stdev - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      0.248
 Maximum value       =      1.659
 Mean value          =      1.541
 Standard Deviation  =      0.159
 Variance            =      0.025
10 - Name Z.Simu.Y.Pb.1 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.871
 Standard Deviation  =      2.959
 Variance            =      8.757
11 - Name Z.Simu.Y.Pb.2 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.270
 Standard Deviation  =      2.785
 Variance            =      7.759
12 - Name Z.Simu.Y.Pb.3 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.532
 Standard Deviation  =      2.904
 Variance            =      8.434
13 - Name Z.Simu.Y.Pb.4 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.264
 Standard Deviation  =      2.835
 Variance            =      8.036
14 - Name Z.Simu.Y.Pb.5 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.125
 Standard Deviation  =      2.548
 Variance            =      6.492
15 - Name Z.Simu.Y.Pb.6 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.397
 Standard Deviation  =      2.731
 Variance            =      7.456
16 - Name Z.Simu.Y.Pb.7 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.174
 Standard Deviation  =      2.709
 Variance            =      7.340
17 - Name Z.Simu.Y.Pb.8 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.504
 Standard Deviation  =      2.871
 Variance            =      8.245
18 - Name Z.Simu.Y.Pb.9 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      7.261
 Standard Deviation  =      3.079
 Variance            =      9.481
19 - Name Z.Simu.Y.Pb.10 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      6.340
 Standard Deviation  =      2.852
 Variance            =      8.133
20 - Name Stats.MEAN - Locator z1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.438
 Maximum value       =     11.597
 Mean value          =      6.474
 Standard Deviation  =      1.063
 Variance            =      1.131

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
Column = 2 - Name = Neigh.Pb.Number - Locator = NA
Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA
Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA
Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA
Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA
Column = 7 - Name = Kriging.Pb.estim - Locator = NA
Column = 8 - Name = Kriging.Pb.stdev - Locator = NA
Column = 9 - Name = Z.Simu.Y.Pb.1 - Locator = NA
Column = 10 - Name = Z.Simu.Y.Pb.2 - Locator = NA
Column = 11 - Name = Z.Simu.Y.Pb.3 - Locator = NA
Column = 12 - Name = Z.Simu.Y.Pb.4 - Locator = NA
Column = 13 - Name = Z.Simu.Y.Pb.5 - Locator = NA
Column = 14 - Name = Z.Simu.Y.Pb.6 - Locator = NA
Column = 15 - Name = Z.Simu.Y.Pb.7 - Locator = NA
Column = 16 - Name = Z.Simu.Y.Pb.8 - Locator = NA
Column = 17 - Name = Z.Simu.Y.Pb.9 - Locator = NA
Column = 18 - Name = Z.Simu.Y.Pb.10 - Locator = NA
Column = 19 - Name = Stats.MEAN - Locator = z1
 

Displaying the average of the Conditional Simulations

In [42]:
if graphics:
    ax = mygrid.plot("Stats*MEAN")
    ax = mydb.plot("Pb")
    ax.decoration(title="Mean of Pb simulations")

Multivariate case¶

The Gaussian transform of the Pb variable has already been calculated. It suffices to perform the Gaussian transform of the Zn variable

In [43]:
mydb.setLocator("Zn",gl.ELoc.Z)
myanamZn = gl.AnamHermite(30)
myanamZn.fitFromLocator(mydb)
if verbose:
    myanamZn.display()
Hermitian Anamorphosis
----------------------
Minimum absolute value for Y  = -2.5
Maximum absolute value for Y  = 2.6
Minimum absolute value for Z  = 1.1469
Maximum absolute value for Z  = 12.1276
Minimum practical value for Y = -2.5
Maximum practical value for Y = 2.6
Minimum practical value for Z = 1.1469
Maximum practical value for Z = 12.1276
Mean                          = 2.88061
Variance                      = 2.76263
Number of Hermite polynomials = 30
Normalized coefficients for Hermite polynomials (punctual variable)
               [,  0]    [,  1]    [,  2]    [,  3]    [,  4]    [,  5]    [,  6]
     [  0,]     2.881    -1.277     0.877    -0.447    -0.095     0.294    -0.121
     [  7,]    -0.087     0.134    -0.029    -0.087     0.069     0.034    -0.065
     [ 14,]     0.005     0.044    -0.026    -0.020     0.034     0.001    -0.033
     [ 21,]     0.010     0.027    -0.016    -0.019     0.016     0.012    -0.014
     [ 28,]    -0.005     0.011
 
In [44]:
if graphics:
    gp.anam(myanamZn)

We convert the raw data into its Gaussian equivalent

In [45]:
mydb.setLocator("Zn",gl.ELoc.Z)
err = myanamZn.rawToGaussianByLocator(mydb)
if verbose:
    mydb.display()
Data Base Characteristics
=========================

Data Base Summary
-----------------
File is organized as a set of isolated points
Space dimension              = 2
Number of Columns            = 10
Total number of samples      = 102
Number of active samples     = 99

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = X - Locator = x1
Column = 2 - Name = Y - Locator = x2
Column = 3 - Name = Zn - Locator = NA
Column = 4 - Name = Pb - Locator = NA
Column = 5 - Name = sel - Locator = sel
Column = 6 - Name = Xvalid.Pb.esterr - Locator = NA
Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA
Column = 8 - Name = Y.Pb - Locator = NA
Column = 9 - Name = Y.Zn - Locator = z1
 

We now perform the multivariate variogram calculation

In [46]:
mydb.setLocators(["Y.Pb","Y.Zn"],gl.ELoc.Z)
myvario = gl.Vario(myvarioParam)
err = myvario.compute(mydb,gl.ECalcVario.VARIOGRAM)
mymodelM = gl.Model.createFromDb(mydb)
err = mymodelM.fit(myvario,[gl.ECov.EXPONENTIAL])
if graphics:
    axs = gp.varmod(myvario,mymodelM)
    gp.decoration(axs,title="Multivariate Model")
    gp.geometry(axs,dims=[5,5])

We perform 10 bivariate conditional simulations (deleting the previous monovariate simulation outcomes first for better legibility)

In [47]:
mygrid.deleteColumn("Z.Simu*")
err = gl.simtub(mydb,mygrid,mymodelM,myneigh,nbsimu=10)
if verbose:
    mygrid.display()
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 30
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
Column = 2 - Name = Neigh.Pb.Number - Locator = NA
Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA
Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA
Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA
Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA
Column = 7 - Name = Kriging.Pb.estim - Locator = NA
Column = 8 - Name = Kriging.Pb.stdev - Locator = NA
Column = 9 - Name = Stats.MEAN - Locator = NA
Column = 10 - Name = Simu.Y.Pb.1 - Locator = z1
Column = 11 - Name = Simu.Y.Pb.2 - Locator = z2
Column = 12 - Name = Simu.Y.Pb.3 - Locator = z3
Column = 13 - Name = Simu.Y.Pb.4 - Locator = z4
Column = 14 - Name = Simu.Y.Pb.5 - Locator = z5
Column = 15 - Name = Simu.Y.Pb.6 - Locator = z6
Column = 16 - Name = Simu.Y.Pb.7 - Locator = z7
Column = 17 - Name = Simu.Y.Pb.8 - Locator = z8
Column = 18 - Name = Simu.Y.Pb.9 - Locator = z9
Column = 19 - Name = Simu.Y.Pb.10 - Locator = z10
Column = 20 - Name = Simu.Y.Zn.1 - Locator = z11
Column = 21 - Name = Simu.Y.Zn.2 - Locator = z12
Column = 22 - Name = Simu.Y.Zn.3 - Locator = z13
Column = 23 - Name = Simu.Y.Zn.4 - Locator = z14
Column = 24 - Name = Simu.Y.Zn.5 - Locator = z15
Column = 25 - Name = Simu.Y.Zn.6 - Locator = z16
Column = 26 - Name = Simu.Y.Zn.7 - Locator = z17
Column = 27 - Name = Simu.Y.Zn.8 - Locator = z18
Column = 28 - Name = Simu.Y.Zn.9 - Locator = z19
Column = 29 - Name = Simu.Y.Zn.10 - Locator = z20
 

We back-transform each set of simulation outcomes using its own Gaussian Anamorphosis function. Finally we delete the Gaussian variables and ask for the statistics on the simulated variables in the Raw Scale.

In [48]:
err = myanamZn.gaussianToRaw(mygrid,"Simu.Y.Zn*")
err = myanamPb.gaussianToRaw(mygrid,"Simu.Y.Pb*")
mygrid.deleteColumn("Simu.Y*")
mygrid.deleteColumn("Stats.*")
mygrid.statisticsBySample(["Z.Simu.*"],[gl.EStatOption.MINI,
                                gl.EStatOption.MAXI,
                                gl.EStatOption.MEAN,
                                gl.EStatOption.STDV],
                    True)
if verbose:
    dbfmt = gl.DbStringFormat()
    dbfmt.setFlags(flag_resume=True, flag_stats=True) 
    mygrid.display(dbfmt)
Data Base Grid Characteristics
==============================

Data Base Summary
-----------------
File is organized as a regular grid
Space dimension              = 2
Number of Columns            = 33
Total number of samples      = 5100

Grid characteristics:
---------------------
Origin :    107.850   481.660
Mesh   :      0.500     0.500
Number :         75        68

Data Base Statistics
--------------------
1 - Name x1 - Locator x1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    107.850
 Maximum value       =    144.850
 Mean value          =    126.350
 Standard Deviation  =     10.824
 Variance            =    117.167
2 - Name x2 - Locator x2
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =    481.660
 Maximum value       =    515.160
 Mean value          =    498.410
 Standard Deviation  =      9.814
 Variance            =     96.313
3 - Name Neigh.Pb.Number - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      6.000
 Mean value          =      5.381
 Standard Deviation  =      1.538
 Variance            =      2.366
4 - Name Neigh.Pb.MaxDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.935
 Maximum value       =      9.999
 Mean value          =      5.598
 Standard Deviation  =      2.497
 Variance            =      6.233
5 - Name Neigh.Pb.MinDist - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.014
 Maximum value       =      9.978
 Mean value          =      3.484
 Standard Deviation  =      2.514
 Variance            =      6.320
6 - Name Neigh.Pb.NbNESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      1.000
 Maximum value       =      1.000
 Mean value          =      1.000
 Standard Deviation  =      0.000
 Variance            =      0.000
7 - Name Neigh.Pb.NbCESect - Locator NA
 Nb of data          =       5100
 Nb of active values =       4596
 Minimum value       =      0.000
 Maximum value       =      0.000
 Mean value          =      0.000
 Standard Deviation  =      0.000
 Variance            =      0.000
8 - Name Kriging.Pb.estim - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.412
 Maximum value       =     11.456
 Mean value          =      6.123
 Standard Deviation  =      0.648
 Variance            =      0.419
9 - Name Kriging.Pb.stdev - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      0.248
 Maximum value       =      1.659
 Mean value          =      1.541
 Standard Deviation  =      0.159
 Variance            =      0.025
10 - Name Z.Simu.Y.Zn.1 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.805
 Standard Deviation  =      1.580
 Variance            =      2.498
11 - Name Z.Simu.Y.Zn.2 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.733
 Standard Deviation  =      1.408
 Variance            =      1.983
12 - Name Z.Simu.Y.Zn.3 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.781
 Standard Deviation  =      1.566
 Variance            =      2.453
13 - Name Z.Simu.Y.Zn.4 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.907
 Standard Deviation  =      1.630
 Variance            =      2.658
14 - Name Z.Simu.Y.Zn.5 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.900
 Standard Deviation  =      1.551
 Variance            =      2.405
15 - Name Z.Simu.Y.Zn.6 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.919
 Standard Deviation  =      1.629
 Variance            =      2.654
16 - Name Z.Simu.Y.Zn.7 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.840
 Standard Deviation  =      1.500
 Variance            =      2.251
17 - Name Z.Simu.Y.Zn.8 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.851
 Standard Deviation  =      1.607
 Variance            =      2.581
18 - Name Z.Simu.Y.Zn.9 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.760
 Standard Deviation  =      1.395
 Variance            =      1.947
19 - Name Z.Simu.Y.Zn.10 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =     12.128
 Mean value          =      2.934
 Standard Deviation  =      1.701
 Variance            =      2.895
20 - Name Z.Simu.Y.Pb.1 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.745
 Standard Deviation  =      1.770
 Variance            =      3.133
21 - Name Z.Simu.Y.Pb.2 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.728
 Standard Deviation  =      1.769
 Variance            =      3.130
22 - Name Z.Simu.Y.Pb.3 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.822
 Standard Deviation  =      1.811
 Variance            =      3.281
23 - Name Z.Simu.Y.Pb.4 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.634
 Standard Deviation  =      1.715
 Variance            =      2.940
24 - Name Z.Simu.Y.Pb.5 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.639
 Standard Deviation  =      1.642
 Variance            =      2.695
25 - Name Z.Simu.Y.Pb.6 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.742
 Standard Deviation  =      1.714
 Variance            =      2.939
26 - Name Z.Simu.Y.Pb.7 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.742
 Standard Deviation  =      1.853
 Variance            =      3.435
27 - Name Z.Simu.Y.Pb.8 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.699
 Standard Deviation  =      1.751
 Variance            =      3.065
28 - Name Z.Simu.Y.Pb.9 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.726
 Standard Deviation  =      1.784
 Variance            =      3.184
29 - Name Z.Simu.Y.Pb.10 - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.003
 Maximum value       =     12.978
 Mean value          =      5.784
 Standard Deviation  =      1.837
 Variance            =      3.374
30 - Name Stats.MINI - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      1.147
 Maximum value       =      4.252
 Mean value          =      1.839
 Standard Deviation  =      0.287
 Variance            =      0.082
31 - Name Stats.MAXI - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      3.300
 Maximum value       =     12.978
 Mean value          =      9.189
 Standard Deviation  =      2.051
 Variance            =      4.206
32 - Name Stats.MEAN - Locator NA
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      2.386
 Maximum value       =      9.627
 Mean value          =      4.285
 Standard Deviation  =      0.491
 Variance            =      0.241
33 - Name Stats.STDV - Locator z1
 Nb of data          =       5100
 Nb of active values =       5100
 Minimum value       =      0.498
 Maximum value       =      4.139
 Mean value          =      2.093
 Standard Deviation  =      0.485
 Variance            =      0.235

Variables
---------
Column = 0 - Name = x1 - Locator = x1
Column = 1 - Name = x2 - Locator = x2
Column = 2 - Name = Neigh.Pb.Number - Locator = NA
Column = 3 - Name = Neigh.Pb.MaxDist - Locator = NA
Column = 4 - Name = Neigh.Pb.MinDist - Locator = NA
Column = 5 - Name = Neigh.Pb.NbNESect - Locator = NA
Column = 6 - Name = Neigh.Pb.NbCESect - Locator = NA
Column = 7 - Name = Kriging.Pb.estim - Locator = NA
Column = 8 - Name = Kriging.Pb.stdev - Locator = NA
Column = 9 - Name = Z.Simu.Y.Zn.1 - Locator = NA
Column = 10 - Name = Z.Simu.Y.Zn.2 - Locator = NA
Column = 11 - Name = Z.Simu.Y.Zn.3 - Locator = NA
Column = 12 - Name = Z.Simu.Y.Zn.4 - Locator = NA
Column = 13 - Name = Z.Simu.Y.Zn.5 - Locator = NA
Column = 14 - Name = Z.Simu.Y.Zn.6 - Locator = NA
Column = 15 - Name = Z.Simu.Y.Zn.7 - Locator = NA
Column = 16 - Name = Z.Simu.Y.Zn.8 - Locator = NA
Column = 17 - Name = Z.Simu.Y.Zn.9 - Locator = NA
Column = 18 - Name = Z.Simu.Y.Zn.10 - Locator = NA
Column = 19 - Name = Z.Simu.Y.Pb.1 - Locator = NA
Column = 20 - Name = Z.Simu.Y.Pb.2 - Locator = NA
Column = 21 - Name = Z.Simu.Y.Pb.3 - Locator = NA
Column = 22 - Name = Z.Simu.Y.Pb.4 - Locator = NA
Column = 23 - Name = Z.Simu.Y.Pb.5 - Locator = NA
Column = 24 - Name = Z.Simu.Y.Pb.6 - Locator = NA
Column = 25 - Name = Z.Simu.Y.Pb.7 - Locator = NA
Column = 26 - Name = Z.Simu.Y.Pb.8 - Locator = NA
Column = 27 - Name = Z.Simu.Y.Pb.9 - Locator = NA
Column = 28 - Name = Z.Simu.Y.Pb.10 - Locator = NA
Column = 29 - Name = Stats.MINI - Locator = NA
Column = 30 - Name = Stats.MAXI - Locator = NA
Column = 31 - Name = Stats.MEAN - Locator = NA
Column = 32 - Name = Stats.STDV - Locator = z1
 

Categorical Variable¶

We compare the initial variable 'Pb' with a set of disjoint intervals. The 'Pb' values varying from 3 to 12.7, we consider three classes:

  • values below 4
  • values between 4 and 6
  • values above 6

We first build the indicators for each class

In [49]:
limits = gl.Limits([np.nan, 4., 6., np.nan])
if verbose:
    limits.display()
Bound( 1 ) : ] -Inf ; 4 [
Bound( 2 ) : [ 4 ; 6 [
Bound( 3 ) : [ 6 ;  +Inf [
 

We apply the set of limits previously defined in order to transform the input variable into Indicators of the different classes.

In [50]:
err = limits.toIndicator(mydb,"Pb")
if verbose:
    mydb.display()
Data Base Characteristics
=========================

Data Base Summary
-----------------
File is organized as a set of isolated points
Space dimension              = 2
Number of Columns            = 13
Total number of samples      = 102
Number of active samples     = 99

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = X - Locator = x1
Column = 2 - Name = Y - Locator = x2
Column = 3 - Name = Zn - Locator = NA
Column = 4 - Name = Pb - Locator = NA
Column = 5 - Name = sel - Locator = sel
Column = 6 - Name = Xvalid.Pb.esterr - Locator = NA
Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA
Column = 8 - Name = Y.Pb - Locator = NA
Column = 9 - Name = Y.Zn - Locator = NA
Column = 10 - Name = Indicator.Pb.Class.1 - Locator = z1
Column = 11 - Name = Indicator.Pb.Class.2 - Locator = z2
Column = 12 - Name = Indicator.Pb.Class.3 - Locator = z3
 

We calculate the variogram of the Indicators for future use

In [51]:
myvarioindParam = gl.VarioParam()
myvarioindParam.addDir(mydir)
myvarioInd = gl.Vario(myvarioindParam)
err = myvarioInd.compute(mydb,gl.ECalcVario.VARIOGRAM)
if verbose:
    myvarioInd.display()
Variogram characteristics
=========================
Number of variable(s)       = 3
Number of direction(s)      = 1
Space dimension             = 2
Variance-Covariance Matrix
               [,  0]    [,  1]    [,  2]
     [  0,]     0.107    -0.062    -0.044
     [  1,]    -0.062     0.250    -0.187
     [  2,]    -0.044    -0.187     0.231

Direction #1
------------
Number of lags              = 10
Direction coefficients      =      1.000     0.000
Direction angles (degrees)  =      0.000     0.000
Tolerance on direction      =     90.000 (degrees)
Calculation lag             =      1.000
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0     3.000     0.389     0.000
         1   123.000     1.081     0.081
         2   183.000     2.038     0.126
         3   205.000     3.006     0.156
         4   231.000     4.013     0.132
         5   229.000     5.036     0.159
         6   198.000     5.962     0.152
         7   187.000     7.000     0.107
         8   204.000     7.996     0.096
         9   184.000     8.990     0.068

For variables 2 and 1
      Rank    Npairs  Distance     Value
         0     3.000     0.389     0.000
         1   123.000     1.081    -0.065
         2   183.000     2.038    -0.077
         3   205.000     3.006    -0.085
         4   231.000     4.013    -0.093
         5   229.000     5.036    -0.085
         6   198.000     5.962    -0.061
         7   187.000     7.000    -0.045
         8   204.000     7.996    -0.042
         9   184.000     8.990    -0.038

For variable 2
      Rank    Npairs  Distance     Value
         0     3.000     0.389     0.167
         1   123.000     1.081     0.199
         2   183.000     2.038     0.221
         3   205.000     3.006     0.251
         4   231.000     4.013     0.292
         5   229.000     5.036     0.258
         6   198.000     5.962     0.237
         7   187.000     7.000     0.254
         8   204.000     7.996     0.228
         9   184.000     8.990     0.234

For variables 3 and 1
      Rank    Npairs  Distance     Value
         0     3.000     0.389     0.000
         1   123.000     1.081    -0.016
         2   183.000     2.038    -0.049
         3   205.000     3.006    -0.071
         4   231.000     4.013    -0.039
         5   229.000     5.036    -0.074
         6   198.000     5.962    -0.091
         7   187.000     7.000    -0.061
         8   204.000     7.996    -0.054
         9   184.000     8.990    -0.030

For variables 3 and 2
      Rank    Npairs  Distance     Value
         0     3.000     0.389    -0.167
         1   123.000     1.081    -0.134
         2   183.000     2.038    -0.145
         3   205.000     3.006    -0.166
         4   231.000     4.013    -0.199
         5   229.000     5.036    -0.172
         6   198.000     5.962    -0.177
         7   187.000     7.000    -0.209
         8   204.000     7.996    -0.186
         9   184.000     8.990    -0.196

For variable 3
      Rank    Npairs  Distance     Value
         0     3.000     0.389     0.167
         1   123.000     1.081     0.150
         2   183.000     2.038     0.194
         3   205.000     3.006     0.237
         4   231.000     4.013     0.238
         5   229.000     5.036     0.247
         6   198.000     5.962     0.268
         7   187.000     7.000     0.270
         8   204.000     7.996     0.240
         9   184.000     8.990     0.226
 
In [52]:
ax = gp.varmod(myvarioInd)
gp.geometry(ax,dims=[4,4])

Then we build a categorical variable which gives the index of the class to which each sample belongs

In [53]:
err = limits.toCategory(mydb,"Pb")
if verbose:
    dbfmt = gl.DbStringFormat()
    dbfmt.setFlags(flag_stats=True)
    dbfmt.setNames(["Category*"])
    dbfmt.setMode(2)
    mydb.display(dbfmt)
Data Base Characteristics
=========================

Data Base Summary
-----------------
File is organized as a set of isolated points
Space dimension              = 2
Number of Columns            = 14
Total number of samples      = 102
Number of active samples     = 99

Data Base Statistics
--------------------
14 - Name Category.Pb - Locator z1
 Nb of data          =        102
 Nb of active values =         99
 Class         1 =         12 (    12.121%)
 Class         2 =         51 (    51.515%)
 Class         3 =         36 (    36.364%)

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = X - Locator = x1
Column = 2 - Name = Y - Locator = x2
Column = 3 - Name = Zn - Locator = NA
Column = 4 - Name = Pb - Locator = NA
Column = 5 - Name = sel - Locator = sel
Column = 6 - Name = Xvalid.Pb.esterr - Locator = NA
Column = 7 - Name = Xvalid.Pb.stderr - Locator = NA
Column = 8 - Name = Y.Pb - Locator = NA
Column = 9 - Name = Y.Zn - Locator = NA
Column = 10 - Name = Indicator.Pb.Class.1 - Locator = NA
Column = 11 - Name = Indicator.Pb.Class.2 - Locator = NA
Column = 12 - Name = Indicator.Pb.Class.3 - Locator = NA
Column = 13 - Name = Category.Pb - Locator = z1
 
In [ ]: