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_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
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(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(nameRaster="*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 Known Mean(s) 0.000
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.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)
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(nameRaster="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(nameRaster="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.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
In [28]:
if graphics:
ax = mygrid.plot(nameRaster="Kriging.Pb.estim")
ax = mydb.plot(nameColor="Pb")
ax.decoration(title="Estimate of Pb")
In [29]:
if graphics:
ax = mygrid.plot(nameRaster="Kriging.Pb.stdev")
ax = mydb.plot(nameColor="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.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.1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.657 Maximum value = 5.385 Mean value = 0.096 Standard Deviation = 1.558 Variance = 2.428 11 - Name Simu.Y.Pb.2 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.567 Maximum value = 5.479 Mean value = -0.016 Standard Deviation = 1.485 Variance = 2.204 12 - Name Simu.Y.Pb.3 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -6.331 Maximum value = 5.947 Mean value = 0.202 Standard Deviation = 1.538 Variance = 2.366 13 - Name Simu.Y.Pb.4 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -4.987 Maximum value = 5.198 Mean value = 0.197 Standard Deviation = 1.509 Variance = 2.277 14 - Name Simu.Y.Pb.5 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.212 Maximum value = 6.757 Mean value = 0.253 Standard Deviation = 1.560 Variance = 2.432 15 - Name Simu.Y.Pb.6 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.511 Maximum value = 5.955 Mean value = 0.357 Standard Deviation = 1.530 Variance = 2.342 16 - Name Simu.Y.Pb.7 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.226 Maximum value = 5.411 Mean value = 0.206 Standard Deviation = 1.471 Variance = 2.163 17 - Name Simu.Y.Pb.8 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -4.746 Maximum value = 5.969 Mean value = 0.315 Standard Deviation = 1.501 Variance = 2.252 18 - Name Simu.Y.Pb.9 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.383 Maximum value = 6.563 Mean value = 0.181 Standard Deviation = 1.603 Variance = 2.570 19 - Name Simu.Y.Pb.10 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -5.238 Maximum value = 5.652 Mean value = 0.164 Standard Deviation = 1.542 Variance = 2.379 20 - Name Stats.MINI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -6.331 Maximum value = 1.009 Mean value = -2.095 Standard Deviation = 0.937 Variance = 0.879 21 - Name Stats.MAXI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -0.578 Maximum value = 6.757 Mean value = 2.507 Standard Deviation = 1.001 Variance = 1.003 22 - Name Stats.MEAN - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = -2.536 Maximum value = 2.182 Mean value = 0.195 Standard Deviation = 0.569 Variance = 0.324 23 - Name Stats.STDV - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.184 Maximum value = 2.840 Mean value = 1.379 Standard Deviation = 0.355 Variance = 0.126 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(nameRaster="Simu.Y.Pb.1")
ax = mydb.plot(nameColor="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.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.1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.202 Standard Deviation = 2.727 Variance = 7.434 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 = 5.950 Standard Deviation = 2.481 Variance = 6.155 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.360 Standard Deviation = 2.727 Variance = 7.438 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.333 Standard Deviation = 2.704 Variance = 7.314 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.456 Standard Deviation = 2.822 Variance = 7.963 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.621 Standard Deviation = 2.872 Variance = 8.249 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.332 Standard Deviation = 2.671 Variance = 7.137 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.542 Standard Deviation = 2.812 Variance = 7.907 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 = 6.333 Standard Deviation = 2.814 Variance = 7.916 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.306 Standard Deviation = 2.732 Variance = 7.466 20 - Name Stats.MINI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 7.224 Mean value = 3.400 Standard Deviation = 0.534 Variance = 0.285 21 - Name Stats.MAXI - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 4.561 Maximum value = 12.978 Mean value = 11.169 Standard Deviation = 2.196 Variance = 4.824 22 - Name Stats.MEAN - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.527 Maximum value = 10.444 Mean value = 6.343 Standard Deviation = 0.998 Variance = 0.997 23 - Name Stats.STDV - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 0.291 Maximum value = 4.395 Mean value = 2.445 Standard Deviation = 0.742 Variance = 0.551 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(nameRaster="Z.Simu.Y.Pb.1")
ax = mydb.plot(nameColor="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.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.1 - Locator NA Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.003 Maximum value = 12.978 Mean value = 6.202 Standard Deviation = 2.727 Variance = 7.434 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 = 5.950 Standard Deviation = 2.481 Variance = 6.155 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.360 Standard Deviation = 2.727 Variance = 7.438 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.333 Standard Deviation = 2.704 Variance = 7.314 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.456 Standard Deviation = 2.822 Variance = 7.963 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.621 Standard Deviation = 2.872 Variance = 8.249 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.332 Standard Deviation = 2.671 Variance = 7.137 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.542 Standard Deviation = 2.812 Variance = 7.907 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 = 6.333 Standard Deviation = 2.814 Variance = 7.916 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.306 Standard Deviation = 2.732 Variance = 7.466 20 - Name Stats.MEAN - Locator z1 Nb of data = 5100 Nb of active values = 5100 Minimum value = 3.527 Maximum value = 10.444 Mean value = 6.343 Standard Deviation = 0.998 Variance = 0.997 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(nameRaster="Stats*MEAN")
ax = mydb.plot(nameColor="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.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.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 [ ]: