Standard 2-D case study¶
Import packages¶
import numpy as np
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
Global variables
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
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
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
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)
if graphics:
gp.plot(mydb, nameColor="Pb")
gp.decoration(title="Data Set")
Variograms¶
We first define the geometry of the variogram calculations
myVarioParamOmni = gl.VarioParam()
mydir = gl.DirParam.create(nlag=10, dlag=1.0)
myVarioParamOmni.addDir(mydir)
We use the variogram definition in order to calculate the variogram cloud.
dbcloud = gl.db_vcloud(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
if graphics:
gp.plot(dbcloud, "Cloud*")
gp.decoration(title="Variogram Cloud")
Calculating the experimental omni-directional variogram
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
Variable(s) = [Pb]
Variance-Covariance Matrix 2.881
Direction #1
------------
Number of lags = 10
Direction coefficients = 1.000 0.000
Direction angles (degrees) = 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
if graphics:
gp.plot(myVarioOmni)
gp.decoration(title="Omni-directional Variogram for Pb")
Calculate a variogram in several directions
myvarioParam = gl.VarioParam()
mydirs = gl.DirParam.createMultiple(4, 10, 1.0)
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
Variable(s) = [Pb]
Variance-Covariance Matrix 2.881
Direction #1
------------
Number of lags = 10
Direction coefficients = 1.000 0.000
Direction angles (degrees) = 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
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
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
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
if graphics:
gp.plot(myvario, idir=-1)
gp.decoration(title="Multi-Directional Variogram of Pb")
Calculating the Variogram Map
myvmap = gl.db_vmap(mydb, calculType=gl.ECalcVario.VARIOGRAM, nxx=[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
if graphics:
gp.plot(myvmap, "*Var")
gp.decoration(title="Variogram Map")
Model¶
Fitting a Model. We call the Automatic Fitting procedure providing the list of covariance functions to be tested.
mymodel = gl.Model.createFromDb(mydb)
err = mymodel.fit(myvario, [gl.ECov.EXPONENTIAL, gl.ECov.SPHERICAL])
Visualizing the resulting model, overlaid on the experimental variogram
if graphics:
gp.varmod(myvario, mymodel)
gp.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).
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.0))
err = myModelConstrained.fit(myVarioOmni, [gl.ECov.SPHERICAL], constr)
if err > 0:
print("Error while fitting the model")
myModelConstrained
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 Known Mean(s) 0.000
We can impose inequality constraints by using EConsType.LOWER or EConsType.UPPER.
Adding a drift¶
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.039
- Ranges = 2.103 0.386
- Theo. Ranges = 0.702 0.129
- Angles = 44.853 0.000
- Rotation Matrix
[, 0] [, 1]
[ 0,] -0.709 0.705
[ 1,] -0.705 -0.709
Spherical
- Sill = 1.606
- Ranges = 6.909 5.009
- Angles = 136.493 0.000
- Rotation Matrix
[, 0] [, 1]
[ 0,] -0.725 -0.688
[ 1,] 0.688 -0.725
Total Sill = 2.645
Drift Part
----------
Universality_Condition
Defining the Neighborhood¶
We initiate a Neigborhood (Moving with a small number of samples for Demonstration)
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
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.
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
if graphics:
gp.plot(mygrid, "Neigh*Number")
gp.decoration(title="Number of Samples per Neighborhood")
- the one giving the maximum distance per neighborhood
if graphics:
gp.plot(mygrid, "Neigh*MaxDist")
gp.decoration(title="Maximum Distance per Neighborhood")
Cross-validation¶
We can now process the cross-validation step
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
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.
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.408 Maximum value = 11.493 Mean value = 6.114 Standard Deviation = 0.641 Variance = 0.410 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.234 Maximum value = 1.655 Mean value = 1.540 Standard Deviation = 0.160 Variance = 0.026 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
if graphics:
gp.plot(mygrid, "Kriging.Pb.estim")
gp.plot(mydb, nameColor="Pb")
gp.decoration(title="Estimate of Pb")
if graphics:
gp.plot(mygrid, "Kriging.Pb.stdev")
gp.plot(mydb, nameColor="Pb")
gp.decoration(title="St. Deviation of Pb")
Simulations¶
We must first transform the Data into Gaussian
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.
if graphics:
gp.plot(myanamPb)
The next step consists in translating the target variable ('Pb') into its Gaussian transform
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
myvarioParam = gl.VarioParam()
mydir = gl.DirParam(10, 1.0)
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).
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.
err = gl.simtub(mydb, mygrid, mymodel, 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 = 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.S1 - Locator = z1 Column = 10 - Name = Simu.Y.Pb.S2 - Locator = z2 Column = 11 - Name = Simu.Y.Pb.S3 - Locator = z3 Column = 12 - Name = Simu.Y.Pb.S4 - Locator = z4 Column = 13 - Name = Simu.Y.Pb.S5 - Locator = z5 Column = 14 - Name = Simu.Y.Pb.S6 - Locator = z6 Column = 15 - Name = Simu.Y.Pb.S7 - Locator = z7 Column = 16 - Name = Simu.Y.Pb.S8 - Locator = z8 Column = 17 - Name = Simu.Y.Pb.S9 - Locator = z9 Column = 18 - Name = Simu.Y.Pb.S10 - Locator = z10
Some statistics on the Conditional simulations in Gaussian scale
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.408 Maximum value = 11.493 Mean value = 6.114 Standard Deviation = 0.641 Variance = 0.410 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.234 Maximum value = 1.655 Mean value = 1.540 Standard Deviation = 0.160 Variance = 0.026 10 - Name Simu.Y.Pb.S1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.146 Maximum value = 5.909 Mean value = 0.295 Standard Deviation = 1.646 Variance = 2.708 11 - Name Simu.Y.Pb.S2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.775 Maximum value = 5.453 Mean value = -0.094 Standard Deviation = 1.546 Variance = 2.392 12 - Name Simu.Y.Pb.S3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.496 Maximum value = 6.225 Mean value = 0.119 Standard Deviation = 1.584 Variance = 2.509 13 - Name Simu.Y.Pb.S4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -4.897 Maximum value = 6.261 Mean value = 0.433 Standard Deviation = 1.587 Variance = 2.518 14 - Name Simu.Y.Pb.S5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -4.759 Maximum value = 6.695 Mean value = 0.420 Standard Deviation = 1.588 Variance = 2.521 15 - Name Simu.Y.Pb.S6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.215 Maximum value = 7.209 Mean value = 0.304 Standard Deviation = 1.603 Variance = 2.569 16 - Name Simu.Y.Pb.S7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -4.939 Maximum value = 6.154 Mean value = 0.370 Standard Deviation = 1.535 Variance = 2.358 17 - Name Simu.Y.Pb.S8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.058 Maximum value = 6.995 Mean value = 0.323 Standard Deviation = 1.568 Variance = 2.459 18 - Name Simu.Y.Pb.S9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -4.184 Maximum value = 5.380 Mean value = 0.399 Standard Deviation = 1.507 Variance = 2.271 19 - Name Simu.Y.Pb.S10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.266 Maximum value = 6.048 Mean value = 0.011 Standard Deviation = 1.591 Variance = 2.533 20 - Name Stats.MINI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.775 Maximum value = 0.847 Mean value = -2.099 Standard Deviation = 0.936 Variance = 0.876 21 - Name Stats.MAXI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -1.201 Maximum value = 7.209 Mean value = 2.681 Standard Deviation = 1.024 Variance = 1.048 22 - Name Stats.MEAN - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -2.125 Maximum value = 2.241 Mean value = 0.258 Standard Deviation = 0.577 Variance = 0.333 23 - Name Stats.STDV - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.267 Maximum value = 2.773 Mean value = 1.431 Standard Deviation = 0.365 Variance = 0.134 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.S1 - Locator = NA Column = 10 - Name = Simu.Y.Pb.S2 - Locator = NA Column = 11 - Name = Simu.Y.Pb.S3 - Locator = NA Column = 12 - Name = Simu.Y.Pb.S4 - Locator = NA Column = 13 - Name = Simu.Y.Pb.S5 - Locator = NA Column = 14 - Name = Simu.Y.Pb.S6 - Locator = NA Column = 15 - Name = Simu.Y.Pb.S7 - Locator = NA Column = 16 - Name = Simu.Y.Pb.S8 - Locator = NA Column = 17 - Name = Simu.Y.Pb.S9 - Locator = NA Column = 18 - Name = Simu.Y.Pb.S10 - 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
if graphics:
gp.plot(mygrid, "Simu.Y.Pb.S1")
gp.plot(mydb, nameColor="Pb")
gp.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.
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.S1 - Locator = z1 Column = 10 - Name = Z.Simu.Y.Pb.S2 - Locator = z2 Column = 11 - Name = Z.Simu.Y.Pb.S3 - Locator = z3 Column = 12 - Name = Z.Simu.Y.Pb.S4 - Locator = z4 Column = 13 - Name = Z.Simu.Y.Pb.S5 - Locator = z5 Column = 14 - Name = Z.Simu.Y.Pb.S6 - Locator = z6 Column = 15 - Name = Z.Simu.Y.Pb.S7 - Locator = z7 Column = 16 - Name = Z.Simu.Y.Pb.S8 - Locator = z8 Column = 17 - Name = Z.Simu.Y.Pb.S9 - Locator = z9 Column = 18 - Name = Z.Simu.Y.Pb.S10 - Locator = z10
We calculate some statistics on the Conditional Simulations in Raw scale.
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.408 Maximum value = 11.493 Mean value = 6.114 Standard Deviation = 0.641 Variance = 0.410 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.234 Maximum value = 1.655 Mean value = 1.540 Standard Deviation = 0.160 Variance = 0.026 10 - Name Z.Simu.Y.Pb.S1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.593 Standard Deviation = 2.989 Variance = 8.931 11 - Name Z.Simu.Y.Pb.S2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.889 Standard Deviation = 2.596 Variance = 6.741 12 - Name Z.Simu.Y.Pb.S3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.252 Standard Deviation = 2.771 Variance = 7.676 13 - Name Z.Simu.Y.Pb.S4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.777 Standard Deviation = 2.993 Variance = 8.957 14 - Name Z.Simu.Y.Pb.S5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.718 Standard Deviation = 2.901 Variance = 8.415 15 - Name Z.Simu.Y.Pb.S6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.600 Standard Deviation = 2.911 Variance = 8.474 16 - Name Z.Simu.Y.Pb.S7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.665 Standard Deviation = 2.893 Variance = 8.371 17 - Name Z.Simu.Y.Pb.S8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.545 Standard Deviation = 2.842 Variance = 8.079 18 - Name Z.Simu.Y.Pb.S9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.700 Standard Deviation = 2.870 Variance = 8.239 19 - Name Z.Simu.Y.Pb.S10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.088 Standard Deviation = 2.727 Variance = 7.439 20 - Name Stats.MINI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 6.887 Mean value = 3.400 Standard Deviation = 0.536 Variance = 0.287 21 - Name Stats.MAXI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.942 Maximum value = 12.978 Mean value = 11.503 Standard Deviation = 2.070 Variance = 4.283 22 - Name Stats.MEAN - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.340 Maximum value = 10.577 Mean value = 6.483 Standard Deviation = 1.017 Variance = 1.035 23 - Name Stats.STDV - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.323 Maximum value = 4.549 Mean value = 2.577 Standard Deviation = 0.734 Variance = 0.538 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.S1 - Locator = NA Column = 10 - Name = Z.Simu.Y.Pb.S2 - Locator = NA Column = 11 - Name = Z.Simu.Y.Pb.S3 - Locator = NA Column = 12 - Name = Z.Simu.Y.Pb.S4 - Locator = NA Column = 13 - Name = Z.Simu.Y.Pb.S5 - Locator = NA Column = 14 - Name = Z.Simu.Y.Pb.S6 - Locator = NA Column = 15 - Name = Z.Simu.Y.Pb.S7 - Locator = NA Column = 16 - Name = Z.Simu.Y.Pb.S8 - Locator = NA Column = 17 - Name = Z.Simu.Y.Pb.S9 - Locator = NA Column = 18 - Name = Z.Simu.Y.Pb.S10 - 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
if graphics:
gp.plot(mygrid, "Z.Simu.Y.Pb.S1")
gp.plot(mydb, nameColor="Pb")
gp.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.
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.408 Maximum value = 11.493 Mean value = 6.114 Standard Deviation = 0.641 Variance = 0.410 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.234 Maximum value = 1.655 Mean value = 1.540 Standard Deviation = 0.160 Variance = 0.026 10 - Name Z.Simu.Y.Pb.S1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.593 Standard Deviation = 2.989 Variance = 8.931 11 - Name Z.Simu.Y.Pb.S2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.889 Standard Deviation = 2.596 Variance = 6.741 12 - Name Z.Simu.Y.Pb.S3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.252 Standard Deviation = 2.771 Variance = 7.676 13 - Name Z.Simu.Y.Pb.S4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.777 Standard Deviation = 2.993 Variance = 8.957 14 - Name Z.Simu.Y.Pb.S5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.718 Standard Deviation = 2.901 Variance = 8.415 15 - Name Z.Simu.Y.Pb.S6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.600 Standard Deviation = 2.911 Variance = 8.474 16 - Name Z.Simu.Y.Pb.S7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.665 Standard Deviation = 2.893 Variance = 8.371 17 - Name Z.Simu.Y.Pb.S8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.545 Standard Deviation = 2.842 Variance = 8.079 18 - Name Z.Simu.Y.Pb.S9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.700 Standard Deviation = 2.870 Variance = 8.239 19 - Name Z.Simu.Y.Pb.S10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.088 Standard Deviation = 2.727 Variance = 7.439 20 - Name Stats.MEAN - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.340 Maximum value = 10.577 Mean value = 6.483 Standard Deviation = 1.017 Variance = 1.035 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.S1 - Locator = NA Column = 10 - Name = Z.Simu.Y.Pb.S2 - Locator = NA Column = 11 - Name = Z.Simu.Y.Pb.S3 - Locator = NA Column = 12 - Name = Z.Simu.Y.Pb.S4 - Locator = NA Column = 13 - Name = Z.Simu.Y.Pb.S5 - Locator = NA Column = 14 - Name = Z.Simu.Y.Pb.S6 - Locator = NA Column = 15 - Name = Z.Simu.Y.Pb.S7 - Locator = NA Column = 16 - Name = Z.Simu.Y.Pb.S8 - Locator = NA Column = 17 - Name = Z.Simu.Y.Pb.S9 - Locator = NA Column = 18 - Name = Z.Simu.Y.Pb.S10 - Locator = NA Column = 19 - Name = Stats.MEAN - Locator = z1
Displaying the average of the Conditional Simulations
if graphics:
gp.plot(mygrid, "Stats*MEAN")
gp.plot(mydb, nameColor="Pb")
gp.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
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
if graphics:
gp.plot(myanamZn)
We convert the raw data into its Gaussian equivalent
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
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(title="Multivariate Model")
We perform 10 bivariate conditional simulations (deleting the previous monovariate simulation outcomes first for better legibility)
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.S1 - Locator = z1 Column = 11 - Name = Simu.Y.Pb.S2 - Locator = z2 Column = 12 - Name = Simu.Y.Pb.S3 - Locator = z3 Column = 13 - Name = Simu.Y.Pb.S4 - Locator = z4 Column = 14 - Name = Simu.Y.Pb.S5 - Locator = z5 Column = 15 - Name = Simu.Y.Pb.S6 - Locator = z6 Column = 16 - Name = Simu.Y.Pb.S7 - Locator = z7 Column = 17 - Name = Simu.Y.Pb.S8 - Locator = z8 Column = 18 - Name = Simu.Y.Pb.S9 - Locator = z9 Column = 19 - Name = Simu.Y.Pb.S10 - Locator = z10 Column = 20 - Name = Simu.Y.Zn.S1 - Locator = z11 Column = 21 - Name = Simu.Y.Zn.S2 - Locator = z12 Column = 22 - Name = Simu.Y.Zn.S3 - Locator = z13 Column = 23 - Name = Simu.Y.Zn.S4 - Locator = z14 Column = 24 - Name = Simu.Y.Zn.S5 - Locator = z15 Column = 25 - Name = Simu.Y.Zn.S6 - Locator = z16 Column = 26 - Name = Simu.Y.Zn.S7 - Locator = z17 Column = 27 - Name = Simu.Y.Zn.S8 - Locator = z18 Column = 28 - Name = Simu.Y.Zn.S9 - Locator = z19 Column = 29 - Name = Simu.Y.Zn.S10 - 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.
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.408 Maximum value = 11.493 Mean value = 6.114 Standard Deviation = 0.641 Variance = 0.410 9 - Name Kriging.Pb.stdev - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.234 Maximum value = 1.655 Mean value = 1.540 Standard Deviation = 0.160 Variance = 0.026 10 - Name Z.Simu.Y.Zn.S1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.775 Standard Deviation = 1.509 Variance = 2.276 11 - Name Z.Simu.Y.Zn.S2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.931 Standard Deviation = 1.737 Variance = 3.016 12 - Name Z.Simu.Y.Zn.S3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.984 Standard Deviation = 1.780 Variance = 3.169 13 - Name Z.Simu.Y.Zn.S4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.830 Standard Deviation = 1.567 Variance = 2.455 14 - Name Z.Simu.Y.Zn.S5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.845 Standard Deviation = 1.598 Variance = 2.555 15 - Name Z.Simu.Y.Zn.S6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.883 Standard Deviation = 1.752 Variance = 3.071 16 - Name Z.Simu.Y.Zn.S7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.991 Standard Deviation = 1.769 Variance = 3.130 17 - Name Z.Simu.Y.Zn.S8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.871 Standard Deviation = 1.660 Variance = 2.756 18 - Name Z.Simu.Y.Zn.S9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.925 Standard Deviation = 1.686 Variance = 2.843 19 - Name Z.Simu.Y.Zn.S10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 12.128 Mean value = 2.829 Standard Deviation = 1.511 Variance = 2.283 20 - Name Z.Simu.Y.Pb.S1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.724 Standard Deviation = 1.762 Variance = 3.106 21 - Name Z.Simu.Y.Pb.S2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.793 Standard Deviation = 1.811 Variance = 3.281 22 - Name Z.Simu.Y.Pb.S3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.833 Standard Deviation = 1.721 Variance = 2.961 23 - Name Z.Simu.Y.Pb.S4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.842 Standard Deviation = 1.732 Variance = 2.999 24 - Name Z.Simu.Y.Pb.S5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.730 Standard Deviation = 1.784 Variance = 3.183 25 - Name Z.Simu.Y.Pb.S6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.535 Standard Deviation = 1.602 Variance = 2.566 26 - Name Z.Simu.Y.Pb.S7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.788 Standard Deviation = 1.783 Variance = 3.179 27 - Name Z.Simu.Y.Pb.S8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.848 Standard Deviation = 1.941 Variance = 3.767 28 - Name Z.Simu.Y.Pb.S9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.887 Standard Deviation = 1.814 Variance = 3.291 29 - Name Z.Simu.Y.Pb.S10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 5.840 Standard Deviation = 1.791 Variance = 3.208 30 - Name Stats.MINI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 1.147 Maximum value = 4.178 Mean value = 1.835 Standard Deviation = 0.286 Variance = 0.082 31 - Name Stats.MAXI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.445 Maximum value = 12.978 Mean value = 9.193 Standard Deviation = 2.062 Variance = 4.252 32 - Name Stats.MEAN - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 2.376 Maximum value = 9.400 Mean value = 4.334 Standard Deviation = 0.595 Variance = 0.354 33 - Name Stats.STDV - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.579 Maximum value = 3.887 Mean value = 2.116 Standard Deviation = 0.477 Variance = 0.227 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.S1 - Locator = NA Column = 10 - Name = Z.Simu.Y.Zn.S2 - Locator = NA Column = 11 - Name = Z.Simu.Y.Zn.S3 - Locator = NA Column = 12 - Name = Z.Simu.Y.Zn.S4 - Locator = NA Column = 13 - Name = Z.Simu.Y.Zn.S5 - Locator = NA Column = 14 - Name = Z.Simu.Y.Zn.S6 - Locator = NA Column = 15 - Name = Z.Simu.Y.Zn.S7 - Locator = NA Column = 16 - Name = Z.Simu.Y.Zn.S8 - Locator = NA Column = 17 - Name = Z.Simu.Y.Zn.S9 - Locator = NA Column = 18 - Name = Z.Simu.Y.Zn.S10 - Locator = NA Column = 19 - Name = Z.Simu.Y.Pb.S1 - Locator = NA Column = 20 - Name = Z.Simu.Y.Pb.S2 - Locator = NA Column = 21 - Name = Z.Simu.Y.Pb.S3 - Locator = NA Column = 22 - Name = Z.Simu.Y.Pb.S4 - Locator = NA Column = 23 - Name = Z.Simu.Y.Pb.S5 - Locator = NA Column = 24 - Name = Z.Simu.Y.Pb.S6 - Locator = NA Column = 25 - Name = Z.Simu.Y.Pb.S7 - Locator = NA Column = 26 - Name = Z.Simu.Y.Pb.S8 - Locator = NA Column = 27 - Name = Z.Simu.Y.Pb.S9 - Locator = NA Column = 28 - Name = Z.Simu.Y.Pb.S10 - 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
limits = gl.Limits([np.nan, 4.0, 6.0, 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.
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.S1 - Locator = z1 Column = 11 - Name = Indicator.Pb.Class.S2 - Locator = z2 Column = 12 - Name = Indicator.Pb.Class.S3 - Locator = z3
We calculate the variogram of the Indicators for future use
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
Variable(s) = [Indicator.Pb.Class.S1 Indicator.Pb.Class.S2 Indicator.Pb.Class.S3]
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
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
ax = gp.varmod(myvarioInd)
Then we build a categorical variable which gives the index of the class to which each sample belongs
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.S1 - Locator = NA Column = 11 - Name = Indicator.Pb.Class.S2 - Locator = NA Column = 12 - Name = Indicator.Pb.Class.S3 - Locator = NA Column = 13 - Name = Category.Pb - Locator = z1