This file is meant to demonstrate the use of gstlearn for Potential Model.
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
return false;
}
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 matplotlib.pyplot as plt
import plotly.graph_objects as go
Setting some global variables
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.
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
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
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.
range = gl.CovAniso.scale2range(gl.ECov.GAUSSIAN, 20.)
model = gl.Model.createFromParam(gl.ECov.GAUSSIAN, range)
model.switchToGradient()
The Neighborhood (although parametrized) should be Unique
neighU = gl.NeighUnique.create()
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).
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
ax = gp.grid1D(grid1D,name="Potential")
ax.geometry(dims=[10,6], aspect='auto')
ax = gp.point(dbiso1D, nameCoorY="Potential", nameColor="iso", size=100,
flagGradient=True)
hl = plt.hlines(dbiso1D.getColumn("Potential"),0.,100.)
We define the space dimension
ndim = 2
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)
Defining the Iso-Potential information
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
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
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
nx = [101,101]
dx = [0.1, 0.1]
grid2D = gl.DbGrid.create(nx, dx)
Defining the Model and turn it into a Gradient Model
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
Defining the Neighborhood
neighU = gl.NeighUnique.create()
err = gl.potential_kriging(dbiso2D, dbgrd2D, dbtgt2D, grid2D, model, neighU,
flag_pot = True, flag_save_data = True, verbose = False)
Visualization
grid2D
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 4 Total number of samples = 10201 Grid characteristics: --------------------- Origin : 0.000 0.000 Mesh : 0.100 0.100 Number : 101 101 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = x2 - Locator = x2 Column = 3 - Name = Potential - Locator = z1
gp.setDefaultGeographic(dims=[8,8], aspect='auto')
ax = grid2D.plot()
ax = dbiso2D.plot(nameColor="iso", size=100)
levels = np.unique(np.round(dbiso2D.getColumn("Potential"),4))
ax = grid2D.plot(nameContour="Potential", levels=levels, colors='red')
ax = dbgrd2D.plot(flagGradient=True, colorGradient='black', scaleGradient=5)
levels = np.unique(np.round(dbgrd2D.getColumn("Potential"),4))
ax = grid2D.plot(nameContour="Potential", levels=levels, colors="black")
ax = dbtgt2D.plot(flagTangent=True, colorTangent='blue', scaleTangent=10)
levels = np.unique(np.round(dbtgt2D.getColumn("Potential"),4))
ax = grid2D.plot(nameContour="Potential", levels=levels, colors="blue")
We define the space dimension
ndim = 3
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)
Defining the Iso-Potential information
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
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
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
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
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
neighU = gl.NeighUnique.create()
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
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()