Experimental variogram parameters¶

This tutorial is meant to describe the aim of the parameters for computing the experimental variogram (or covariance, ...)

Let us first recall that a Variogram is the calculation of a function, applied to a Data Base, and following a set of rules contained in a VarioParam. The VarioParam is a set of rules defined per direction; the characteristices of each Direction are stored in a DirParam.

In this script, we will consider a data set (in 2D and in 3D), define the parametrization of a single direction, and check all the samples which will be kept for comparison with the target sample iech.

In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import matplotlib.pyplot as plt
import numpy as np

import plotly.graph_objects as go
import gstlearn.plot3D as gop

import gstlearn.document as gdoc

gdoc.setNoScroll()
figsize = (8, 8)

2D case¶

We define a 2D data set with samples generated randomly within a square of size 1

In [2]:
ndat = 500
db = gl.Db.createFillRandom(ndat, ndim=2, nvar=1)

Visuallize the data set

In [3]:
fix, ax = plt.subplots(figsize=figsize)
ax.symbol(db, nameSize="z", flagCst=True)
ax.geometry(aspect=1)
ax.decoration(title="Sample locations")
plt.show()
No description has been provided for this image

Find the rank of a sample close to the center of the plot

In [4]:
center = [0.5, 0.5]
iech0 = db.getSampleClosestTo(center)

Look for the samples connected to the central point through the variogram criterium.

In [5]:
ilag0 = 2
dlag = 0.15
codir = [1.0, 0.0]
toldis = 0.5
tolang = 20.0
bench = 0.07
cylrad = 0.1
ranks = gl.variogramPerPoint(
    db, iech0, ilag0, dlag, codir, toldis, tolang, bench, cylrad
)
print("Number of samples selected =", len(ranks))

Number of samples selected = 21

Graphic representation of the samples selected for the calculation of the variogram:

  • for a given direction (defined by the vector codir)
  • for a given variogram lag
  • for a given choise of cylrad parameter
  • for a given choice of bench parameter
In [6]:
xtab0 = db.getSamplesOneCoordinate([iech0], idim=0)
ytab0 = db.getSamplesOneCoordinate([iech0], idim=1)

xtab = db.getSamplesOneCoordinate(ranks, idim=0)
ytab = db.getSamplesOneCoordinate(ranks, idim=1)

angles = gl.GH.rotationGetAngles(codir)
mini, center, maxi = gp.lagDefine(ilag0, dlag, toldis * 100)
In [7]:
fix, ax = plt.subplots(figsize=figsize)
ax.symbol(db, nameSize="z", flagCst=True)
ax.decoration(title="Sample locations")
ax.geometry(aspect=1)

plt.scatter(xtab0, ytab0, s=100, label="Selected pairs", c="black")
plt.scatter(xtab, ytab, s=50, label="Selected pairs", c="blue", marker="s")

gp.drawCircles(mini, maxi, True, xtab0, ytab0)

gp.drawCone(angles[0], tolang=tolang, x0=xtab0, y0=ytab0)

gp.drawCylrad(angles[0], cylrad, x0=xtab0, y0=ytab0)

gp.drawBench(angles[0], bench, x0=xtab0, y0=ytab0)

plt.show()
No description has been provided for this image

3D case¶

We define a 3D data set with samples generated randomly within a cube of size 1

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

ndat = 500
db = gl.Db.createFillRandom(ndat, ndim=ndim, nvar=1)

Visualize the data set

In [9]:
points = gop.PointDb(db, size=1, color="red")
fig = go.Figure(data=[points])
fig.update_layout(
    width=1000,  # largeur en pixels
    height=800,  # hauteur en pixels
)
f = fig.show()

Select the sample closest to the centrum of the cube.

In [10]:
center = [0.5, 0.5, 0.5]
iech0 = db.getSampleClosestTo(center)

Look for the samples connected to the central point through the variogram criterium.

In [11]:
ilag0 = 2
dlag = 0.2
codir = [0.0, 1.0, 1.0]
toldis = 0.5
tolang = 30.0
bench = 0.4
cylrad = 0.1
ranks = gl.variogramPerPoint(
    db, iech0, ilag0, dlag, codir, toldis, tolang, bench, cylrad
)
print("Number of samples selected =", len(ranks))

Number of samples selected = 8

Graphic representation of the samples selected for the calculation of the variogram:

  • for a given direction (defined by the vector codir)
  • for a given variogram lag
  • for a given choise of cylrad parameter
  • for a given choice of bench parameter
In [12]:
xtab0 = db.getSamplesOneCoordinate([iech0], idim=0)
ytab0 = db.getSamplesOneCoordinate([iech0], idim=1)
ztab0 = db.getSamplesOneCoordinate([iech0], idim=2)
apex = tuple(np.concatenate([xtab0, ytab0, ztab0]))
axis = tuple(codir)

xtab = db.getSamplesOneCoordinate(ranks, idim=0)
ytab = db.getSamplesOneCoordinate(ranks, idim=1)
ztab = db.getSamplesOneCoordinate(ranks, idim=2)

angles = gl.GH.rotationGetAngles(codir)
mini, center, maxi = gp.lagDefine(ilag0, dlag, toldis * 100)
In [13]:
opacity = 0.1
points = gop.PointDb(db, size=1, color="red")
target = gop.Scatter(xtab0, ytab0, ztab0, mode="markers", m_size=3, m_color="black")
neighs = gop.Scatter(xtab, ytab, ztab, mode="markers", m_size=2, m_color="blue")
cone1 = gop.Cone(
    apex, axis, tolang, orientation=+1, length=maxi, color="blue", opacity=opacity
)
cone2 = gop.Cone(
    apex, axis, tolang, orientation=-1, length=maxi, color="blue", opacity=opacity
)
sphere1 = gop.Sphere(mini, apex, color="red", opacity=0.6)
sphere2 = gop.Sphere(maxi, apex, color="red", opacity=opacity)
cylinder = gop.Cylinder(
    apex, axis, cylrad, length=2 * maxi, color="green", opacity=opacity
)
bench1 = gop.Bench(
    apex, bench, orientation=+1, length=2 * maxi, color="orange", opacity=opacity
)
bench2 = gop.Bench(
    apex, bench, orientation=-1, length=2 * maxi, color="orange", opacity=opacity
)

fig = go.Figure(
    data=[
        t
        for t in [
            points,
            target,
            neighs,
            sphere1,
            sphere2,
            cone1,
            cone2,
            cylinder,
            bench1,
            bench2,
        ]
        if t is not None
    ]
)

fig.update_layout(
    width=1000,  # largeur en pixels
    height=800,  # hauteur en pixels
    showlegend=False,
)
f = fig.show()