Positive Definite Variogram Models#
Michael J. Pyrcz, Professor, The University of Texas at Austin
Twitter | GitHub | Website | GoogleScholar | Book | YouTube | Applied Geostats in Python e-book | LinkedIn
Chapter of e-book “Applied Geostatistics in Python: a Hands-on Guide with GeostatsPy”.
Cite this e-Book as:
Pyrcz, M.J., 2024, Applied Geostatistics in Python: a Hands-on Guide with GeostatsPy, https://geostatsguy.github.io/GeostatsPyDemos_Book.
The workflows in this book and more are available here:
Cite the GeostatsPyDemos GitHub Repository as:
Pyrcz, M.J., 2024, GeostatsPyDemos: GeostatsPy Python Package for Spatial Data Analytics and Geostatistics Demonstration Workflows Repository (0.0.1). Zenodo. https://zenodo.org/doi/10.5281/zenodo.12667035
By Michael J. Pyrcz
© Copyright 2024.
This chapter is a tutorial for / demonstration of Permissible Variogram Models available in GeostatsPy.
YouTube Lecture: check out my lectures on:
For your convenience here’s a summary of salient points.
Spatial Continuity#
Spatial Continuity is the correlation between values over distance.
No spatial continuity – no correlation between values over distance, random values at each location in space regardless of separation distance.
Homogenous phenomenon have perfect spatial continuity, since all values as the same (or very similar) they are correlated.
We need a statistic to quantify spatial continuity! A convenient method is the Semivariogram.
The Variogram#
Function of difference over distance.
The expected (average) squared difference between values separated by a lag distance vector (distance and direction), \(h\):
where \(z(\bf{u}_\alpha)\) and \(z(\bf{u}_\alpha + \bf{h})\) are the spatial sample values at tail and head locations of the lag vector respectively.
Calculated over a suite of lag distances to obtain a continuous function.
the \(\frac{1}{2}\) term converts a variogram into a semivariogram, but in practice the term variogram is used instead of semivariogram.
We prefer the semivariogram because it relates directly to the covariance function, \(C_x(\bf{h})\) and univariate variance, \(\sigma^2_x\):
Note the correlogram is related to the covariance function as:
The correlogram provides of function of the \(\bf{h}-\bf{h}\) scatter plot correlation vs. lag offset \(\bf{h}\).
Nested Variogram Models#
Spatial continuity can be described with nested spatial continuity models:
where \(\Gamma_x(\bf{h})\) is the nested variogram model resulting from the summation of \(nst\) nested variogram structures, \(\gamma_i(\bf{h})\).
if each \(\gamma_i(\bf{h})\) structure is a positive definite variogram model, then the summation, \(\Gamma_x(\bf{h})\), is also positive definite.
Each one of these variogram structures, \(\gamma_i(\bf{h})\), is based on a geometric anisotropy model parameterized by the orientation and range in the major and minor directions. In 2D this is simply an azimuth and ranges, \(azi\), \(a_{maj}\) and \(a_{min}\). Note, the range in the minor direction (orthogonal to the major direction.
The geometric anisotropy model assumes that the range in all off-diagonal directions is based on an ellipse with the major and minor axes aligned with and set to the major and minor for the variogram.
\begin{equation}
\bf{h}i = \sqrt{\left(\frac{r{maj}}{a_{maj_i}}\right)^2 + \left(\frac{r_{maj}}{a_{maj_i}}\right)^2}
\end{equation}
Therefore, if we know the major direction, range in major and minor directions, we may completely describe each nested component of the complete spatial continuity of the variable of interest, \(i = 1,\dots,nst\).
In this workflow we will observe the common permissible positive definite variogram models in GeostatsPy that we can combine to build our nested variogram models.
for all we assume a contribution of the sill (single structure) and the sill is 1.0 (standardized feature).
Load the Required Libraries#
The following code loads the required libraries.
import geostatspy.GSLIB as GSLIB # GSLIB utilities, visualization and wrapper
import geostatspy.geostats as geostats # GSLIB methods convert to Python
import geostatspy
print('GeostatsPy version: ' + str(geostatspy.__version__))
GeostatsPy version: 0.0.72
We will also need some standard packages. These should have been installed with Anaconda 3.
import os # set working directory, run executables
from tqdm import tqdm # suppress the status bar
from functools import partialmethod
tqdm.__init__ = partialmethod(tqdm.__init__, disable=True)
ignore_warnings = True # ignore warnings?
import numpy as np # ndarrays for gridded data
import pandas as pd # DataFrames for tabular data
import matplotlib.pyplot as plt # for plotting
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks
plt.rc('axes', axisbelow=True) # plot all grids below the plot elements
if ignore_warnings == True:
import warnings
warnings.filterwarnings('ignore')
cmap = plt.cm.inferno # color map
Define Functions#
This is a convenience function to add major and minor gridlines to improve plot interpretability.
def add_grid():
plt.gca().grid(True, which='major',linewidth = 1.0); plt.gca().grid(True, which='minor',linewidth = 0.2) # add y grids
plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)
plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks
List of Lag Distances#
We define a list of lag distances for our variogram models.
h = np.linspace(0,100,101) # lag distance vector
Nugget Effect Variogram Model#
No spatial correlation
Does not have a range, nor directionality, i.e., acts over all distances and directions.
Should be a small component of the overall variance
Very uncommon in siliciclastic sedimentary systems
More common for mineral grades in mining
May be measurement error
The equation:
where \(𝑐_1\) is the contribution, and \(h\) is the lag distance
gamma_nugget = np.ones(101) # nugget effect
gamma_nugget[0] = 0
plt.subplot(111)
plot = plt.plot(h[1:],gamma_nugget[1:],color='red',lw=4)
plt.plot([0,100],[1.0,1.0],color='black',ls='--'); plt.annotate('Sill',[5,1.02])
plt.scatter(0,0,color='red',edgecolor='black',s=50)
gca = plt.gca()
plt.xlim([0,100]); plt.ylim([0,1.2])
plt.xlabel(r'Lag Distance, $\bf{h}$'); plt.ylabel(r'Variogram, $\gamma$ ($\bf{h}$)')
plt.title('Nugget Effect Variogram Model'); add_grid()
plt.subplots_adjust(left=0.0, bottom=0.0, right=0.7, top=0.8, wspace=0.2, hspace=0.2); plt.show()
Spherical Variogram Model#
A very commonly observed variogram / spatial continuity form in many settings
Piecewise, beyond the range is equal to the sill
The equation:
where \(𝑐_1\) is the contribution, \(𝑎\) is the range and \(\bf{𝐡}\) is the lag distance
vrange = 80 # range parameter
gamma_sph = 1.0*(1.5*h/vrange) - 0.5*np.power(h/vrange,3) # spherical model
gamma_sph[h > vrange] = 1.0
plt.subplot(111)
plot = plt.plot(h,gamma_sph,color='red',lw=4)
plt.plot([0,100],[1.0,1.0],color='black',ls='--'); plt.annotate('Sill',[5,1.02])
plt.plot([0,vrange*2/3],[0,1.0],color='black',ls='--'); plt.annotate(r'$\frac{2}{3}$ Range',[vrange*2/3,1.02])
plt.plot([vrange,vrange],[0.0,1.0],color='black',ls='--'); plt.annotate('Range',[vrange-3.5,0.4],rotation=90.0)
gca = plt.gca()
plt.xlim([0,100]); plt.ylim([0,1.2])
plt.xlabel(r'Lag Distance, $\bf{h}$'); plt.ylabel(r'Variogram, $\gamma$ ($\bf{h}$)')
plt.title('Spherical Variogram Model'); add_grid()
plt.subplots_adjust(left=0.0, bottom=0.0, right=0.7, top=0.8, wspace=0.2, hspace=0.2); plt.show()
Exponential Variogram Model#
Also very commonly observed variogram / spatial continuity form
Less short-scale continuity than spherical, and reaches sill asymptotically, range is at 95% of the sill
The equation:
where \(𝑐_1\) is the contribution, \(𝑎\) is the range and \(\bf{h}\) is the lag distance.
vrange = 80 # range parameter
gamma_exp = 1.0*(1.0-np.exp(-3*(h/vrange))) # exponential model
plt.subplot(111)
plot = plt.plot(h,gamma_exp,color='red',lw=4)
plt.plot([0,100],[1.0,1.0],color='black',ls='--'); plt.annotate('Sill',[5,1.02])
plt.plot([0,vrange],[0.95,0.95],color='grey',ls='-'); plt.annotate('95% Sill',[5,0.91])
plt.plot([vrange,vrange],[0.0,1.0],color='black',ls='--'); plt.annotate('Range',[vrange-3.5,0.4],rotation=90.0)
plt.plot([0,vrange*1/3],[0,1.0],color='black',ls='--'); plt.annotate(r'$\frac{1}{3}$ Range',[vrange*1/3,1.02])
gca = plt.gca()
plt.xlim([0,100]); plt.ylim([0,1.2])
plt.xlabel(r'Lag Distance, $\bf{h}$'); plt.ylabel(r'Variogram, $\gamma$ ($\bf{h}$)')
plt.title('Exponential Variogram Model'); add_grid()
plt.subplots_adjust(left=0.0, bottom=0.0, right=0.7, top=0.8, wspace=0.2, hspace=0.2); plt.show()
Gaussian Variogram Model#
Less commonly observed variogram / spatial continuity form, e.g., for thickness and elevation
Much more short-scale continuity than spherical, and reaches sill asymptotically, range is at 95% of the sill
The equation:
where \(𝑐_1\) is the contribution, \(𝑎\) is the range and \(\bf{h}\) is the lag distance.
vrange = 80 # range parameter
gamma_gaus = 1.0*(1.0-np.exp(-3*np.power(h/vrange,2))) # Gaussian model
plt.subplot(111)
plot = plt.plot(h,gamma_gaus,color='red',lw=4)
plt.plot([0,100],[1.0,1.0],color='black',ls='--'); plt.annotate('Sill',[5,1.02])
plt.plot([0,vrange],[0.95,0.95],color='grey',ls='-'); plt.annotate('95% Sill',[5,0.91])
plt.plot([vrange,vrange],[0.0,1.0],color='black',ls='--'); plt.annotate('Range',[vrange-3.5,0.4],rotation=90.0)
gca = plt.gca()
plt.xlim([0,100]); plt.ylim([0,1.2])
plt.xlabel(r'Lag Distance, $\bf{h}$'); plt.ylabel(r'Variogram, $\gamma$ ($\bf{h}$)')
plt.title('Gaussian Variogram Model'); add_grid()
plt.subplots_adjust(left=0.0, bottom=0.0, right=0.7, top=0.8, wspace=0.2, hspace=0.2); plt.show()
Want to Work Together?#
I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.
Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I’d be happy to drop by and work with you!
Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!
I can be reached at mpyrcz@austin.utexas.edu.
I’m always happy to discuss,
Michael
Michael Pyrcz, Ph.D., P.Eng. Professor, Cockrell School of Engineering and The Jackson School of Geosciences, The University of Texas at Austin
More Resources Available at: Twitter | GitHub | Website | GoogleScholar | Geostatistics Book | YouTube | Applied Geostats in Python e-book | Applied Machine Learning in Python e-book | LinkedIn
Comments#
This was a basic demonstration of permissible variogram models to support 3D model construction. Much more can be done, I have other demonstrations for modeling workflows with GeostatsPy in the GitHub repository GeostatsPy_Demos.
I hope this is helpful,
Michael