Standard 2-D case studyĀ¶
Import packagesĀ¶
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
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:
ax = mydb.plot(nameColor="Pb")
ax.decoration(title="Data Set")
VariogramsĀ¶
We first define the geometry of the variogram calculations
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.
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:
ax = dbcloud.plot("Cloud*")
ax.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:
ax = myVarioOmni.plot()
ax.decoration(title="Omni-directional Variogram for Pb")
Calculate a variogram in several directions
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 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:
ax = myvario.plot(idir=-1)
ax.decoration(title="Multi-Directional Variogram of Pb")
Calculating the Variogram Map
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
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.
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:
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).
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
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:
ax = mygrid.plot(nameRaster="Neigh*Number")
ax.decoration(title="Number of Samples per Neighborhood")
- the one giving the maximum distance per neighborhood
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
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:
ax = mygrid.plot(nameRaster="Kriging.Pb.estim")
ax = mydb.plot(nameColor="Pb")
ax.decoration(title="Estimate of Pb")
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
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.anam(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.)
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,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
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
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.
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.
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
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.
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 Summar