GitHub

Empirical variograms

Variograms are widely used in geostatistics due to their intimate connection with (cross-)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 ||\x_i - \x_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).

Variograms can plotted with the following options:

GeoStatsFunctions.varioplot — Function
varioplot(γ; [options])

Plot the variogram γ with given options.

Common variogram options:

  • color - color of variogram
  • size - size of variogram
  • maxlag - maximum lag of variogram

Empirical variogram options:

  • pointsize - size of points of variogram
  • showtext - show text counts
  • textsize - size of text counts
  • showhist - show histogram
  • histcolor - color of histogram

Notes

  • This function will only work in the presence of a Makie.jl backend via package extensions in Julia v1.9 or later versions of the language.
source

(Omini)directional

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

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

Parameters

  • 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, EmpiricalVarioplane.

References

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

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.

Optionally, forward parameters for the underlying EmpiricalVariogram.

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

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.

Optionally, forward parameters for 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.)

varioplot(γ)
Example block output

directional variograms along a specific direction:

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

varioplot(γₕ, color = :maroon, histcolor = :maroon)
varioplot!(γᵥ)
Mke.current_figure()
Example block output

or planar variograms over a specific plane:

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

varioplot(γₕ, color = :maroon, histcolor = :maroon)
varioplot!(γᵥ)
Mke.current_figure()
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.

Varioplanes

Variograms estimated along all directions in a given plane of reference are called varioplanes.

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

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 parameters to the underlying EmpiricalVariogram.

source

The varioplane is plotted on a polar axis for all lags and angles:

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

planeplot(γ)
Example block output