Potential Model¶

This file is meant to demonstrate the use of gstlearn for Potential Model.

In [1]:
import numpy as np
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.0, 1.0, 80.0, 1.0, 60.0, 1.0, 40.0, 2.0, 50.0, 2.0]
dbiso1D = gl.Db.createFromSamples(
    5, gl.ELoadBy.SAMPLE, tabiso, ["x", "iso"], ["x1", "layer"]
)

The Gradient information

In [4]:
tabgrd = [0.0, 1.0]
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.0)
modelRef = gl.Model.createFromParam(gl.ECov.GAUSSIAN, range=range)
covp = gl.CovGradientAnalytic(modelRef.getCovAniso(0))
model = gl.Model()
model.setCov(covp);

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 [7]:
gl.OptDbg.setReference(-1)
gl.krigingPotential(
    dbiso1D,
    dbgrd1D,
    None,
    grid1D,
    model,
    flag_pot=True,
    flag_grad=True,
    flag_save_data=True,
    verbose=False,
)
gl.OptDbg.setReference(-1)

Graphic representation

In [8]:
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.0, 100.0)
No description has been provided for this image

Case in 2-D¶

Defining the Information¶

We define the space dimension

In [9]:
ndim = 2
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)

Defining the Iso-Potential information

In [10]:
tabiso = [
    7.0,
    6.0,
    1.0,
    5.0,
    6.0,
    1.0,
    6.0,
    5.0,
    1.0,
    3.0,
    6.0,
    2.0,
    7.0,
    7.0,
    2.0,
    8.0,
    3.0,
    2.0,
    8.0,
    1.0,
    3.0,
    7.0,
    9.0,
    3.0,
    10.0,
    5.0,
    3.0,
    3.0,
    1.0,
    3.0,
]
dbiso2D = gl.Db.createFromSamples(
    10, gl.ELoadBy.SAMPLE, tabiso, ["x", "y", "iso"], ["x1", "x2", "layer"]
)

Defining the Gradient information

In [11]:
tabgrd = [1.0, 6.0, 1.0, 0.0, 9.0, 2.0, -1.0, 1.0, 7.0, 8.0, 0.0, -1]
dbgrd2D = gl.Db.createFromSamples(
    3, gl.ELoadBy.SAMPLE, tabgrd, ["x", "y", "gx", "gy"], ["x1", "x2", "g1", "g2"]
)

Defining the Tangent information

In [12]:
tabtgt = [3.0, 7.0, 1.0, 0.0, 9.0, 7.0, 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 [13]:
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 [14]:
modelRef = gl.Model.createFromParam(gl.ECov.CUBIC, range=8.0)
covp = gl.CovGradientAnalytic(modelRef.getCovAniso(0))
model = gl.Model()
model.setCov(covp);

Launching the estimation¶

In [15]:
err = gl.krigingPotential(
    dbiso2D,
    dbgrd2D,
    dbtgt2D,
    grid2D,
    model,
    flag_pot=True,
    flag_save_data=True,
    verbose=False,
)

Visualization

In [16]:
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 = [0.     0.2837 1.6514]
Potential values at Gradient samples = [0.9826 1.2818 1.6521]
Potential values at Tangent samples = [0.3921 0.9813]
No description has been provided for this image

Case in 3-D¶

Defining the Information¶

We define the space dimension

In [17]:
ndim = 3
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)

Defining the Iso-Potential information

In [18]:
tabiso = [
    0.0,
    8.0,
    0.0,
    1.0,
    0.0,
    12.0,
    -8.0,
    1.0,
    5.0,
    0.0,
    0.0,
    1.0,
    -5.0,
    0.0,
    0.0,
    1.0,
    0.0,
    10.0,
    -5.0,
    1.0,
    0.0,
    5.0,
    -12.0,
    1.0,
    4.0,
    0.0,
    0.0,
    2.0,
    -4.0,
    0.0,
    0.0,
    2.0,
    0.0,
    0.0,
    -6.0,
    2.0,
]
dbiso3D = gl.Db.createFromSamples(
    9, gl.ELoadBy.SAMPLE, tabiso, ["x", "y", "z", "iso"], ["x1", "x2", "x3", "layer"]
)

Defining the Gradient information

In [19]:
tabgrd = [
    5.0,
    0.0,
    0.0,
    -1.0,
    0.0,
    0.0,
    -5.0,
    0.0,
    0.0,
    1.0,
    0.0,
    0.0,
    0.0,
    8.0,
    0.0,
    0.0,
    -1.0,
    0.0,
    0.0,
    5.0,
    0.0,
    0.0,
    -1.0,
    0.0,
    4.0,
    0.0,
    0.0,
    -1.0,
    0.0,
    0.0,
    -4.0,
    0.0,
    0.0,
    1.0,
    0.0,
    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 [20]:
tabtgt = [
    0.0,
    10.0,
    -5.0,
    0.0,
    -1.0,
    1.0,
    0.0,
    5.0,
    -11.0,
    0.0,
    1.0,
    0.0,
    0.0,
    0.0,
    -6.0,
    0.0,
    1.0,
    0.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 [21]:
nx = [61, 81, 61]
dx = [0.2, 0.2, 0.2]
x0 = [-6.0, 0.0, -12.0]
grid3D = gl.DbGrid.create(nx=nx, dx=dx, x0=x0)

Defining the Model (as IRF-1) and turn it into a Gradient Model

In [22]:
modelRef = gl.Model.createFromParam(gl.ECov.CUBIC, range=8.0)
modelRef.setDriftIRF(1)
covp = gl.CovGradientAnalytic(modelRef.getCovAniso(0))
model = gl.Model()
model.setCov(covp);

Launching the estimation¶

In [23]:
err = gl.krigingPotential(
    dbiso3D,
    dbgrd3D,
    dbtgt3D,
    grid3D,
    model,
    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 [24]:
levels = np.unique(np.round(dbiso3D.getColumn("Potential"), 4))

data = [
    gop.IsoSurfaceOnDbGrid(
        grid3D, "Potential", useSel=False, isomin=levels[0], isomax=levels[0]
    ),
    gop.IsoSurfaceOnDbGrid(
        grid3D, "Potential", useSel=False, isomin=levels[1], isomax=levels[1]
    ),
    gop.PointDb(dbiso3D, nameColor="Potential"),
]
fig = go.Figure(data=data)
f = fig.show()
In [ ]: