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.)
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.,100.)
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., 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 [11]:
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 [12]:
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 [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.)
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.,  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 [19]:
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 [20]:
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 [21]:
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 [22]:
modelRef = gl.Model.createFromParam(gl.ECov.CUBIC, range=8.)
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 [ ]: