Potential Model¶
This file is meant to demonstrate the use of gstlearn for Potential Model.
In [1]:
import numpy as np
import pandas as pd
import sys
import os
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.plot3D as gop
import gstlearn.document as gdoc
import matplotlib.pyplot as plt
import plotly.graph_objects as go
gdoc.setNoScroll()
Case in 1-D¶
Defining the Information¶
Setting some global variables
In [2]:
ndim = 1
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)
Defining the data points used as controls
The set of iso-potential data points. Each datum is characterized by its coordinate and an indicator of the set to which it belongs.
In [3]:
tabiso = [30., 1.,
80., 1.,
60., 1.,
40., 2.,
50., 2.]
dbiso1D = gl.Db.createFromSamples(5, gl.ELoadBy.SAMPLE, tabiso, ["x", "iso"], ["x1", "layer"])
The Gradient information
In [4]:
tabgrd = [0., 1.]
dbgrd1D= gl.Db.createFromSamples(1, gl.ELoadBy.SAMPLE, tabgrd, ["x", "g"], ["x1", "g1"])
The 1-D grid covering the field of estimation
In [5]:
nx = [100]
grid1D = gl.DbGrid.create(nx)
The Model is composed of a Gaussian Covariance with a scale 20. The model must be turned into a Gradient model.
In [6]:
range = gl.scale2range(gl.ECov.GAUSSIAN, 20.)
model = gl.Model.createFromParam(gl.ECov.GAUSSIAN, range)
model.switchToGradient()
The Neighborhood (although parametrized) should be Unique
In [7]:
neighU = gl.NeighUnique.create()
Launching the Estimation¶
We launch the calculation of the Potential on the nodes of the grid. The use of OptDbg statement is maintained for spying the process when estimating a target node in particular (-1 stands for no control).
In [8]:
gl.OptDbg.setReference(-1)
gl.potential_kriging(dbiso1D, dbgrd1D, None, grid1D, model, neighU,
flag_pot=True, flag_grad=True, flag_save_data=True, verbose=False)
gl.OptDbg.setReference(-1)
Graphic representation
In [9]:
gp.grid1D(grid1D,name="Potential")
gp.geometry(dims=[10,6], aspect='auto')
gp.symbol(dbiso1D, nameCoorY="Potential", nameColor="iso", s=100)
gp.gradient(dbiso1D, nameCoorY="Potential")
hl = plt.hlines(dbiso1D.getColumn("Potential"),0.,100.)
Case in 2-D¶
Defining the Information¶
We define the space dimension
In [10]:
ndim = 2
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)
Defining the Iso-Potential information
In [11]:
tabiso = [ 7., 6., 1.,
5., 6., 1.,
6., 5., 1.,
3., 6., 2.,
7., 7., 2.,
8., 3., 2.,
8., 1., 3.,
7., 9., 3.,
10., 5., 3.,
3., 1., 3.]
dbiso2D = gl.Db.createFromSamples(10, gl.ELoadBy.SAMPLE, tabiso,["x","y","iso"],["x1","x2","layer"])
Defining the Gradient information
In [12]:
tabgrd = [ 1., 6., 1., 0.,
9., 2., -1., 1.,
7., 8., 0., -1 ]
dbgrd2D = gl.Db.createFromSamples(3, gl.ELoadBy.SAMPLE, tabgrd,["x","y","gx","gy"],["x1","x2","g1","g2"])
Defining the Tangent information
In [13]:
tabtgt = [ 3., 7., 1., 0.,
9., 7., 0.5, -0.5 ]
dbtgt2D = gl.Db.createFromSamples(2, gl.ELoadBy.SAMPLE, tabtgt, ["x","y","tx","ty"],["x1","x2","tangent1","tangent2"])
Defining the output grid
In [14]:
nx = [101,101]
dx = [0.1, 0.1]
grid2D = gl.DbGrid.create(nx, dx)
Defining the Model and turn it into a Gradient Model
In [15]:
model = gl.Model.createFromParam(gl.ECov.CUBIC, 8.)
model.switchToGradient()
model.display()
Model characteristics ===================== (Specific for Handling Gradient) 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 --------------- Cubic - Sill = 1.000 - Range = 8.000 Total Sill = 1.000 Known Mean(s) 0.000
Defining the Neighborhood
In [16]:
neighU = gl.NeighUnique.create()
Launching the estimation¶
In [17]:
err = gl.potential_kriging(dbiso2D, dbgrd2D, dbtgt2D, grid2D, model, neighU,
flag_pot = True, flag_save_data = True, verbose = False)
Visualization
In [18]:
showResults = True
if showResults:
gp.plot(grid2D)
gp.plot(dbiso2D, nameColor="iso", s=100)
if showResults:
levels = np.unique(np.round(dbiso2D.getColumn("Potential"),4))
print(f"Potential value at Intercept families = {levels}")
gp.isoline(grid2D, name="Potential", levels=levels, colors='red')
gp.gradient(dbgrd2D, color='black', scale=20)
if showResults:
levels = np.unique(np.round(dbgrd2D.getColumn("Potential"),4))
print(f"Potential values at Gradient samples = {levels}")
gp.isoline(grid2D, name="Potential", levels=levels, colors="black")
gp.tangent(dbtgt2D, color='blue', scale=20)
if showResults:
levels = np.unique(np.round(dbtgt2D.getColumn("Potential"),4))
print(f"Potential values at Tangent samples = {levels}")
gp.isoline(grid2D, name="Potential", levels=levels, colors="blue")
plt.show()
Potential value at Intercept families = [-1.6514 -0.2837 0. ] Potential values at Gradient samples = [-1.6521 -1.2818 -0.9826] Potential values at Tangent samples = [-0.9813 -0.3921]
Case in 3-D¶
Defining the Information¶
We define the space dimension
In [19]:
ndim = 3
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)
Defining the Iso-Potential information
In [20]:
tabiso = [ 0., 8., 0., 1.,
0., 12., -8., 1.,
5., 0., 0., 1.,
-5., 0., 0., 1.,
0., 10., -5., 1.,
0., 5.,-12., 1.,
4., 0., 0., 2.,
-4., 0., 0., 2.,
0., 0., -6., 2.]
dbiso3D = gl.Db.createFromSamples(9, gl.ELoadBy.SAMPLE, tabiso,["x","y","z","iso"],["x1","x2","x3","layer"])
Defining the Gradient information
In [21]:
tabgrd = [ 5., 0., 0., -1., 0., 0.,
-5., 0., 0., 1., 0., 0.,
0., 8., 0., 0., -1., 0.,
0., 5., 0., 0., -1., 0.,
4., 0., 0., -1., 0., 0.,
-4., 0., 0., 1., 0., 0.]
dbgrd3D = gl.Db.createFromSamples(6, gl.ELoadBy.SAMPLE, tabgrd,["x","y","z","gx","gy","gz"],
["x1","x2","x3","g1","g2","g3"])
Defining the Tangent information
In [22]:
tabtgt = [ 0., 10., -5., 0., -1., 1.,
0., 5., -11., 0., 1., 0.,
0., 0., -6., 0., 1., 0.]
dbtgt3D = gl.Db.createFromSamples(3, gl.ELoadBy.SAMPLE, tabtgt, ["x","y","z","tx","ty","tz"],
["x1","x2","x3","tangent1","tangent2","tangent3"])
Defining the output grid
In [23]:
nx = [61,81,61]
dx = [0.2, 0.2, 0.2]
x0 = [-6., 0., -12.]
grid3D = gl.DbGrid.create(nx=nx, dx=dx, x0=x0)
Defining the Model (as IRF-1) and turn it into a Gradient Model
In [24]:
model = gl.Model.createFromEnvironment(1, ndim)
model.addCovFromParam(gl.ECov.CUBIC, range=10.)
model.setDriftIRF(1)
model.switchToGradient()
model.display()
Model characteristics ===================== (Specific for Handling Gradient) Space dimension = 3 Number of variable(s) = 1 Number of basic structure(s) = 1 Number of drift function(s) = 4 Number of drift equation(s) = 4 Covariance Part --------------- Cubic - Sill = 1.000 - Range = 10.000 Total Sill = 1.000 Drift Part ---------- Universality_Condition Drift:x1 Drift:x2 Drift:x3
Defining the Neighborhood
In [25]:
neighU = gl.NeighUnique.create()
Launching the estimation¶
In [26]:
err = gl.potential_kriging(dbiso3D, dbgrd3D, dbtgt3D, grid3D, model, neighU,
flag_pot = True, flag_save_data = True, verbose = False)
dbiso3D.display()
Data Base Characteristics ========================= Data Base Summary ----------------- File is organized as a set of isolated points Space dimension = 3 Number of Columns = 6 Total number of samples = 9 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x - Locator = x1 Column = 2 - Name = y - Locator = x2 Column = 3 - Name = z - Locator = x3 Column = 4 - Name = iso - Locator = layer Column = 5 - Name = Potential - Locator = z1
3-D Visualization
In [27]:
levels = np.unique(np.round(dbiso3D.getColumn("Potential"),4))
data = [gop.IsoSurfaceOnDbGrid(grid3D, "Potential", useSel=False,
isomin=levels[0], isomax=levels[0]),
gop.PointDb(dbiso3D, nameColor="Potential")
]
fig = go.Figure(data=data)
f = fig.show()