GitHub

Variograms

Empirical variograms

Variograms are widely used in geostatistics due to their intimate connection with (co)variance and visual interpretability. The following video explains the concept in detail:

The Matheron's estimator of the empirical variogram is given by

\[\widehat{\gamma_M}(h) = \frac{1}{2|N(h)|} \sum_{(i,j) \in N(h)} (z_i - z_j)^2\]

where $N(h) = \left\{(i,j) \mid ||\p_i - \p_j|| = h\right\}$ is the set of pairs of locations at a distance $h$ and $|N(h)|$ is the cardinality of the set. Alternatively, the robust Cressie's estimator is given by

\[\widehat{\gamma_C}(h) = \frac{1}{2}\frac{\left\{\frac{1}{|N(h)|} \sum_{(i,j) \in N(h)} |z_i - z_j|^{1/2}\right\}^4}{0.457 + \frac{0.494}{|N(h)|} + \frac{0.045}{|N(h)|^2}}\]

Both estimators are available and can be used with general distance functions in order to for example:

  • Model anisotropy (e.g. ellipsoid distance)
  • Perform simulation on sphere (e.g. haversine distance)

Please see Distances.jl for a complete list of distance functions.

The high-performance estimation procedure implemented in the framework can consider all pairs of locations regardless of direction (ominidirectional) or a specified partition of the geospatial data (e.g. directional, planar).

(Omini)directional variograms

GeoStatsFunctions.EmpiricalVariogram — Type
EmpiricalVariogram(data, var₁, var₂=var₁; [options])

Computes the empirical (a.k.a. experimental) omnidirectional (cross-)variogram for variables var₁ and var₂ stored in geospatial data.

Options

  • nlags - number of lags (default to 20)
  • maxlag - maximum lag in length units (default to 1/2 of minimum side of bounding box)
  • distance - custom distance function (default to Euclidean distance)
  • estimator - variogram estimator (default to :matheron estimator)
  • algorithm - accumulation algorithm (default to :ball)

Available estimators:

  • :matheron - simple estimator based on squared differences
  • :cressie - robust estimator based on 4th power of differences

Available algorithms:

  • :full - loop over all pairs of points in the data
  • :ball - loop over all points inside maximum lag ball

All implemented algorithms produce the exact same result. The :ball algorithm is considerably faster when the maximum lag is much smaller than the bounding box of the domain of the data.

See also: DirectionalVariogram, PlanarVariogram, EmpiricalVariogramSurface.

References

source
GeoStatsFunctions.DirectionalVariogram — Function
DirectionalVariogram(direction, data, var₁, var₂=var₁; dtol=1e-6u"m", [options])

Computes the empirical (cross-)variogram for the variables var₁ and var₂ stored in geospatial data along a given direction with band tolerance dtol in length units.

Forwards options to the underlying EmpiricalVariogram.

source
GeoStatsFunctions.PlanarVariogram — Function
PlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6u"m", [options])

Computes the empirical (cross-)variogram for the variables var₁ and var₂ stored in geospatial data along a plane perpendicular to a normal direction with plane tolerance ntol in length units.

Forwards options to the underlying EmpiricalVariogram.

source

Consider the following example image:

using GeoStatsImages

img = geostatsimage("Gaussian30x10")

img |> viewer
Example block output

We can estimate ominidirectional variograms, which consider pairs of points along all directions:

γ = EmpiricalVariogram(img, :Z, maxlag = 50.)

funplot(γ)
Example block output

directional variograms along a specific direction:

γₕ = DirectionalVariogram((1.,0.), img, :Z, maxlag = 50.)
γᵥ = DirectionalVariogram((0.,1.), img, :Z, maxlag = 50.)

fig = funplot(γₕ, color = :maroon, histcolor = :maroon)
funplot!(fig, γᵥ)
Example block output

or planar variograms over a specific plane:

γᵥ = PlanarVariogram((1.,0.), img, :Z, maxlag = 50.)
γₕ = PlanarVariogram((0.,1.), img, :Z, maxlag = 50.)

fig = funplot(γₕ, color = :maroon, histcolor = :maroon)
funplot!(fig, γᵥ)
Example block output
Note

The directional and planar variograms coincide in this example because planes are equal to lines in 2-dimensional space. These concepts are most useful in 3-dimensional space where we may be interested in comparing the horizontal planar range to the vertical directional range.

Empirical surfaces

GeoStatsFunctions.EmpiricalVariogramSurface — Type
EmpiricalVariogramSurface(data, var₁, var₂=var₁;
                          normal=Vec(0,0,1), nangs=50,
                          ptol=0.5u"m", dtol=0.5u"m",
                          [options])

Given a normal direction, estimate the (cross-)variogram of variables var₁ and var₂ along all directions in the corresponding plane of variation.

Optionally, specify the tolerance ptol in length units for the plane partition, the tolerance dtol in length units for the direction partition, the number of angles nangs in the plane, and forward the options to the underlying EmpiricalVariogram.

source
γ = EmpiricalVariogramSurface(img, :Z, maxlag = 50.)

surfplot(γ)
Example block output

Theoretical variograms

We provide various theoretical variogram models from the literature, which can be combined with ellipsoid distances to model geometric anisotropy and combined 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 given points $\p_1,\p_2 \in \R^m$:

\[\gamma(\p_1,\p_2) = \gamma(||\p_1 - \p_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.

Gaussian

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

GeoStatsFunctions.GaussianVariogram — Type
GaussianVariogram(; range, sill, nugget)

A Gaussian variogram with range in length units, and sill and nugget contributions.

GaussianVariogram(; ranges, rotation, sill, nugget)

Alternatively, use multiple ranges and rotation matrix to construct an anisotropic model.

GaussianVariogram(ball; sill, nugget)

Alternatively, use a custom metric ball.

Examples

# isotropic model
GaussianVariogram(range=2.0m)

# anisotropic model
GaussianVariogram(ranges=(1.0m, 2.0m))
source
funplot(GaussianVariogram())
Example block output

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)\]

GeoStatsFunctions.SphericalVariogram — Type
SphericalVariogram(; range, sill, nugget)

A spherical variogram with range in length units, and sill and nugget contributions.

SphericalVariogram(; ranges, rotation, sill, nugget)

Alternatively, use multiple ranges and rotation matrix to construct an anisotropic model.

SphericalVariogram(ball; sill, nugget)

Alternatively, use a custom metric ball.

Examples

# isotropic model
SphericalVariogram(range=2.0m)

# anisotropic model
SphericalVariogram(ranges=(1.0m, 2.0m))
source
funplot(SphericalVariogram())
Example block output

Exponential

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

GeoStatsFunctions.ExponentialVariogram — Type
ExponentialVariogram(; range, sill, nugget)

An exponential variogram with range in length units, and sill and nugget contributions.

ExponentialVariogram(; ranges, rotation, sill, nugget)

Alternatively, use multiple ranges and rotation matrix to construct an anisotropic model.

ExponentialVariogram(ball; sill, nugget)

Alternatively, use a custom metric ball.

Examples

# isotropic model
ExponentialVariogram(range=2.0m)

# anisotropic model
ExponentialVariogram(ranges=(1.0m, 2.0m))
source
funplot(ExponentialVariogram())
Example block output

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)\]

GeoStatsFunctions.MaternVariogram — Type
MaternVariogram(; range, sill, nugget, order)

A Matérn variogram with range in length units, sill and nugget contributions, and order of Bessel function.

MaternVariogram(; ranges, rotation, sill, nugget, order)

Alternatively, use multiple ranges and rotation matrix to construct an anisotropic model.

MaternVariogram(ball; sill, nugget, order)

Alternatively, use a custom metric ball.

Examples

# isotropic model
MaternVariogram(range=2.0m)

# anisotropic model
MaternVariogram(ranges=(1.0m, 2.0m))
source
funplot(MaternVariogram())
Example block output

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)\]

GeoStatsFunctions.CubicVariogram — Type
CubicVariogram(; range, sill, nugget)

A cubic variogram with range in length units, and sill and nugget contributions.

CubicVariogram(; ranges, rotation, sill, nugget)

Alternatively, use multiple ranges and rotation matrix to construct an anisotropic model.

CubicVariogram(ball; sill, nugget)

Alternatively, use a custom metric ball.

Examples

# isotropic model
CubicVariogram(range=2.0m)

# anisotropic model
CubicVariogram(ranges=(1.0m, 2.0m))
source
funplot(CubicVariogram())
Example block output

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)\]

GeoStatsFunctions.PentaSphericalVariogram — Type
PentaSphericalVariogram(; range, sill, nugget)

A pentaspherical variogram with range in length units, and sill and nugget contributions.

PentaSphericalVariogram(; ranges, rotation, sill, nugget)

Alternatively, use multiple ranges and rotation matrix to construct an anisotropic model.

PentaSphericalVariogram(ball; sill, nugget)

Alternatively, use a custom metric ball.

Examples

# isotropic model
PentaSphericalVariogram(range=2.0m)

# anisotropic model
PentaSphericalVariogram(ranges=(1.0m, 2.0m))
source
funplot(PentaSphericalVariogram())
Example block output

Sine hole

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

GeoStatsFunctions.SineHoleVariogram — Type
SineHoleVariogram(; range, sill, nugget)

A sinehole variogram with range in length units, and sill and nugget contributions.

SineHoleVariogram(; ranges, rotation, sill, nugget)

Alternatively, use multiple ranges and rotation matrix to construct an anisotropic model.

SineHoleVariogram(ball; sill, nugget)

Alternatively, use a custom metric ball.

Examples

# isotropic model
SineHoleVariogram(range=2.0m)

# anisotropic model
SineHoleVariogram(ranges=(1.0m, 2.0m))
source
funplot(SineHoleVariogram())
Example block output

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)\]

GeoStatsFunctions.CircularVariogram — Type
CircularVariogram(; range, sill, nugget)

A circular variogram with range in length units, and sill and nugget contributions.

CircularVariogram(; ranges, rotation, sill, nugget)

Alternatively, use multiple ranges and rotation matrix to construct an anisotropic model.

CircularVariogram(ball; sill, nugget)

Alternatively, use a custom metric ball.

Examples

# isotropic model
CircularVariogram(range=2.0m)

# anisotropic model
CircularVariogram(ranges=(1.0m, 2.0m))
source
funplot(CircularVariogram())
Example block output

Power

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

funplot(PowerVariogram())
Example block output

Nugget

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

funplot(NuggetEffect())
Example block output

Anisotropic models

Anisotropic models can be constructed from a list of ranges and rotation matrix from Rotations.jl:

GaussianVariogram(ranges = (3.0, 2.0, 1.0), rotation = RotZXZ(0.0, 0.0, 0.0))
GaussianVariogram
├─ ranges: (3.0 m, 2.0 m, 1.0 m)
├─ rotation: [1.0 -0.0 0.0; 0.0 1.0 -0.0; 0.0 0.0 1.0]
├─ sill: 1.0
└─ nugget: 0.0

Rotation angles from commercial geostatistical software are also provided:

GeoStatsBase.MinesightAngles — Type
MinesightAngles(θ₁, θ₂, θ₃)

MineSight z-x'-y'' intrinsic rotation convention following the right-hand rule. All angles are in degrees and the sign convention is CW, CCW, CW positive.

The first rotation θ₁ is a horizontal rotation around the Z-axis, with positive being clockwise. The second rotation θ₂ is a rotation around the new X-axis, which moves the Y-axis into the desired position, with positive direction of rotation is up. The third rotation θ₃ is a rotation around the new Y-axis, which moves the X-axis into the desired position, with positive direction of rotation is up.

References

source
GeoStatsBase.DatamineAngles — Type
DatamineAngles(θ₁, θ₂, θ₃)

Datamine ZXY rotation convention following the left-hand rule. All angles are in degrees and the signal convention is CW, CW, CW positive. Y is the principal axis.

source
GeoStatsBase.VulcanAngles — Type
VulcanAngles(θ₁, θ₂, θ₃)

Vulcan ZYX rotation convention following the right-hand rule. All angles are in degrees and the signal convention is CW, CCW, CW positive. X is the principal axis.

source
GeoStatsBase.GslibAngles — Type
GslibAngles(θ₁, θ₂, θ₃)

GSLIB z-x'-y'' intrinsic rotation convention following the left-hand rule. All angles are in degrees and the sign convention is CW, CCW, CW positive. Y is the principal axis.

The first rotation θ₁ is a rotation around the Z-axis, this is also called the azimuth. The second rotation θ₂ is a rotation around the new X-axis, this is also called the dip. The third rotation θ₃ is a rotation around the new Y-axis, this is also called the tilt.

References

source

The effect of anisotropy is illustrated below for different rotation angles:

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

viz(geotable.geometry, color = geotable.Z)
Example block output
θ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

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

# anisotropic variogram with different rotation angles
for (i, θ) in enumerate(θs)
  # anisotropic variogram model
  γ = GaussianVariogram(ranges = (20.0, 5.0), rotation = Angle2d(θ))

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

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

fig
Example block output

Composite models

A composite 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:

GeoStatsFunctions.structures — Function
structures(γ)

Return the individual structures of a (possibly composite) variogram as a tuple. The structures are the total nugget, and the coefficients (or contributions) for for the remaining non-trivial structures after normalization (i.e. sill=1, nugget=0).

Examples

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

structures(2γ₁ + 3γ₂)
source
structures(cov)

Return the individual structures of a (possibly composite) covariance as a tuple. The structures are the total nugget, and the coefficients (or contributions) for for the remaining non-trivial structures after normalization (i.e. sill=1, nugget=0).

Examples

cov₁ = GaussianCovariance(nugget=1, sill=2)
cov₂ = SphericalCovariance(nugget=2, sill=3)

structures(2cov₁ + 3cov₂)
source
γ₁ = GaussianVariogram(nugget=1, sill=2)
γ₂ = ExponentialVariogram(nugget=2, sill=3)

# composite 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 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 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(range: 1.0 m, sill: 1.0, nugget: 0.0), ExponentialVariogram(range: 1.0 m, sill: 1.0, nugget: 0.0))

Fitting models

Fitting theoretical variograms to empirical observations is an important modeling step to ensure valid mathematical models of geospatial continuity. Given an empirical variogram, the following function can be used to fit models:

GeoStatsFunctions.fit — Function
fit(F, f, algo=WeightedLeastSquares(); kwargs...)

Fit theoretical geostatistical function of type F to empirical function f using algorithm algo.

Optionally fix theoretical parameters like range, sill and nugget in the kwargs.

Examples

julia> fit(SphericalVariogram, g)
julia> fit(ExponentialVariogram, g)
julia> fit(ExponentialVariogram, g, sill=1.0)
julia> fit(ExponentialVariogram, g, maxsill=1.0)
julia> fit(GaussianVariogram, g, WeightedLeastSquares())
source
fit(Fs, f, algo=WeightedLeastSquares(); kwargs...)

Fit theoretical geostatistical functions of types Fs to empirical function f using algorithm algo and return the one with minimum error.

Examples

julia> fit([SphericalVariogram, ExponentialVariogram], g)
source
fit(F, f, weightfun; kwargs...)

Convenience method that forwards the weighting function weightfun to the WeightedLeastSquares algorithm.

Examples

fit(SphericalVariogram, g, h -> exp(-h))
fit(Variogram, g, h -> exp(-h/100))
source
fit(Variogram, g, algo=WeightedLeastSquares(); kwargs...)

Fit all stationary Variogram models to empirical variogram g using algorithm algo and return the one with minimum error.

Examples

julia> fit(Variogram, g)
julia> fit(Variogram, g, WeightedLeastSquares())
source

Example

# sinusoidal data
𝒟 = georef((Z=[sin(i/2) + sin(j/2) for i in 1:50, j in 1:50],))

# empirical variogram
g = EmpiricalVariogram(𝒟, :Z, maxlag = 25u"m")

funplot(g)
Example block output

We can fit specific models to the empirical variogram:

γ = GeoStatsFunctions.fit(SineHoleVariogram, g)

fig = funplot(g)
funplot!(fig, γ, maxlag = 25u"m")
Example block output

or let the framework find the model with minimum error:

γ = GeoStatsFunctions.fit(Variogram, g)

fig = funplot(g)
funplot!(fig, γ, maxlag = 25u"m")
Example block output

which should be a SineHoleVariogram given that the synthetic data of this example is sinusoidal.

Optionally, we can specify a weighting function to give different weights to the lags:

γ = GeoStatsFunctions.fit(SineHoleVariogram, g, h -> 1 / h^2)

fig = funplot(g)
funplot!(fig, γ, maxlag = 25u"m")
Example block output

Methods

Weighted least squares

GeoStatsFunctions.WeightedLeastSquares — Type
WeightedLeastSquares()
WeightedLeastSquares(w)

Fit theoretical variogram using weighted least squares with weighting function w (e.g. h -> 1/h). If no weighting function is provided, bin counts of empirical variogram are normalized and used as weights.

source