Theoretical variograms

We provide various theoretical variogram models from the literature, which can can be combined with ellipsoid distances to model geometric anisotropy and nested with scalars or matrix coefficients to express multivariate relations.

Models

In an intrinsic isotropic model, the variogram is only a function of the distance between any two points $\x_1,\x_2 \in \R^m$:

\[\gamma(\x_1,\x_2) = \gamma(||\x_1 - \x_2||) = \gamma(h)\]

Under the additional assumption of 2nd-order stationarity, the well-known covariance is directly related via $\gamma(h) = \cov(0) - \cov(h)$. This package implements a few commonly used as well as other more exotic variogram models. Most of these models share a set of default parameters (e.g. sill=1.0, range=1.0), which can be set with keyword arguments.

Functions are provided to query properties of variogram models programmatically:

isstationary(f)

isisotropic(f)

Corresponding covariance models are available for convenience. For example, the GaussianVariogram and the GaussianCovariance are two alternative representations of the same geostatistical behavior.

Gaussian

\[\gamma(h) = (s - n) \left[1 - \exp\left(-3\left(\frac{h}{r}\right)^2\right)\right] + n \cdot \1_{(0,\infty)}(h)\]

GaussianVariogram(; range, sill, nugget)
GaussianVariogram(ball; sill, nugget)

varioplot(GaussianVariogram())

Spherical

\[\gamma(h) = (s - n) \left[\left(\frac{3}{2}\left(\frac{h}{r}\right) + \frac{1}{2}\left(\frac{h}{r}\right)^3\right) \cdot \1_{(0,r)}(h) + \1_{[r,\infty)}(h)\right] + n \cdot \1_{(0,\infty)}(h)\]

SphericalVariogram(; range, sill, nugget)
SphericalVariogram(ball; sill, nugget)

varioplot(SphericalVariogram())

Exponential

\[\gamma(h) = (s - n) \left[1 - \exp\left(-3\left(\frac{h}{r}\right)\right)\right] + n \cdot \1_{(0,\infty)}(h)\]

ExponentialVariogram(; range, sill, nugget)
ExponentialVariogram(ball; sill, nugget)

varioplot(ExponentialVariogram())

Matern

\[\gamma(h) = (s - n) \left[1 - \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\sqrt{2\nu}\ 3\left(\frac{h}{r}\right)\right)^\nu K_\nu\left(\sqrt{2\nu}\ 3\left(\frac{h}{r}\right)\right)\right] + n \cdot \1_{(0,\infty)}(h)\]

MaternVariogram(; range, sill, nugget, order)
MaternVariogram(ball; sill, nugget, order)

varioplot(MaternVariogram())

Cubic

\[\gamma(h) = (s - n) \left[\left(7\left(\frac{h}{r}\right)^2 - \frac{35}{4}\left(\frac{h}{r}\right)^3 + \frac{7}{2}\left(\frac{h}{r}\right)^5 - \frac{3}{4}\left(\frac{h}{r}\right)^7\right) \cdot \1_{(0,r)}(h) + \1_{[r,\infty)}(h)\right] + n \cdot \1_{(0,\infty)}(h)\]

CubicVariogram(; range, sill, nugget)
CubicVariogram(ball; sill, nugget)

varioplot(CubicVariogram())

Pentaspherical

\[\gamma(h) = (s - n) \left[\left(\frac{15}{8}\left(\frac{h}{r}\right) - \frac{5}{4}\left(\frac{h}{r}\right)^3 + \frac{3}{8}\left(\frac{h}{r}\right)^5\right) \cdot \1_{(0,r)}(h) + \1_{[r,\infty)}(h)\right] + n \cdot \1_{(0,\infty)}(h)\]

PentasphericalVariogram(; range, sill, nugget)
PentasphericalVariogram(ball; sill, nugget)

varioplot(PentasphericalVariogram())

Sine hole

\[\gamma(h) = (s - n) \left[1 - \frac{\sin(\pi h / r)}{\pi h / r}\right] + n \cdot \1_{(0,\infty)}(h)\]

SineHoleVariogram(; range, sill, nugget)
SineHoleVariogram(ball; sill, nugget)

varioplot(SineHoleVariogram())

Circular

\[\gamma(h) = (s - n) \left[\left(1 - \frac{2}{\pi} \cos^{-1}\left(\frac{h}{r}\right) + \frac{2h}{\pi r} \sqrt{1 - \frac{h^2}{r^2}} \right) \cdot \1_{(0,r)}(h) + \1_{[r,\infty)}(h)\right] + n \cdot \1_{(0,\infty)}(h)\]

CircularVariogram(; range, sill, nugget)
CircularVariogram(ball; sill, nugget)

varioplot(CircularVariogram())

Power

\[\gamma(h) = sh^a + n \cdot \1_{(0,\infty)}(h)\]

PowerVariogram(; scaling, exponent, nugget)

varioplot(PowerVariogram())

Nugget

\[\gamma(h) = n \cdot \1_{(0,\infty)}(h)\]

NuggetEffect(; nugget)

varioplot(NuggetEffect(1.0))

Anisotropy

Anisotropic models are easily obtained by defining an ellipsoid metric in place of the default Euclidean metric as shown in the following example. First, we create an ellipsoid that specifies the ranges and angles of rotation:

ellipsoid = MetricBall((3.0, 2.0, 1.0), RotZXZ(0.0, 0.0, 0.0))

MetricBall((3.0 m, 2.0 m, 1.0 m), Mahalanobis)

All rotations from Rotations.jl are supported as well as the following additional rotations from commercial geostatistical software:

MinesightAngles(θ₁, θ₂, θ₃)

DatamineAngles(θ₁, θ₂, θ₃)

VulcanAngles(θ₁, θ₂, θ₃)

GslibAngles(θ₁, θ₂, θ₃)

We pass the ellipsoid as the first argument to the variogram model instead of specifying a single range with a keyword argument:

GaussianVariogram(ellipsoid, sill = 2.0)

GaussianVariogram
├─ sill: 2.0
├─ nugget: 0.0
├─ ranges: (3.0 m, 2.0 m, 1.0 m)
└─ distance: Mahalanobis

To illustrate the concept, consider the following 2D data set:

points = rand(Point, 50, crs=Cartesian2D) |> Scale(100)
geotable = georef((Z=rand(50),), points)

viz(geotable.geometry, color = geotable.Z)

We interpolate the data with different ellipsoids by varying the angle of rotation from $0$ to $2\pi$ clockwise:

θs = range(0.0, step = π/4, stop = 2π - π/4)

# initialize figure
fig = Mke.Figure(size = (800, 1600))

# helper function to position subfigures
pos = i -> CartesianIndices((4, 2))[i].I

# Kriging with different angles
for (i, θ) in enumerate(θs)
  # ellipsoid rotated clockwise by angle θ
  e = MetricBall((20.,5.), Angle2d(θ))

  # anisotropic variogram model
  γ = GaussianVariogram(e)

  # domain of interpolation
  grid = CartesianGrid(100, 100)

  # perform interpolation
  sol = geotable |> Interpolate(grid, Kriging(γ))

  # visualize solution at subplot i
  viz(fig[pos(i)...], sol.geometry, color = sol.Z,
    axis = (title = "θ = $(rad2deg(θ))ᵒ",)
  )
end

fig

Nesting

A nested variogram model $\gamma(h) = c_1\cdot\gamma_1(h) + c_2\cdot\gamma_2(h) + \cdots + c_n\cdot\gamma_n(h)$ can be constructed from multiple variogram models, including matrix coefficients. The individual structures can be recovered in canonical form with the structures function:

structures(γ)

γ₁ = GaussianVariogram(nugget=1, sill=2)
γ₂ = SphericalVariogram(nugget=2, sill=3)

γ = 2γ₁ + 3γ₂

cₒ, c, g = structures(γ)

γ₁ = GaussianVariogram(nugget=1, sill=2)
γ₂ = ExponentialVariogram(nugget=2, sill=3)

# nested model with matrix coefficients
γ = [1.0 0.0; 0.0 1.0] * γ₁ + [2.0 0.5; 0.5 3.0] * γ₂

# structures in canonical form
cₒ, c, g = structures(γ)

cₒ # matrix nugget

2×2 Matrix{Float64}:
 5.0  1.0
 1.0  7.0

c # matrix coefficients

([1.0 0.0; 0.0 1.0], [2.0 0.5; 0.5 3.0])

g # normalized structures

(GaussianVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean), ExponentialVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean))