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 estimated along all directions in a given plane of reference are called varioplanes. Both variograms and varioplanes can be plotted directly with the following options:

GeoStatsFunctions.varioplot — Function

varioplot(γ; [options])

Plot the variogram or varioplane γ with given options.

Empirical variogram options:

vcolor - color of variogram
psize - size of points of variogram
ssize - size of segments of variogram
tshow - show text counts
tsize - size of text counts
hshow - show histogram
hcolor - color of histogram

Empirical varioplane options:

vscheme - color scheme of varioplane
rshow - show range of theoretical model
rmodel - theoretical model (e.g. SphericalVariogram)
rcolor - color of range curve

Theoretical variogram options:

maxlag - maximum lag for theoretical model

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.

(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 (default to 1/10 of maximum lag of data)
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.

References

Chilès, JP and Delfiner, P. 2012. Geostatistics: Modeling Spatial Uncertainty
Webster, R and Oliver, MA. 2007. Geostatistics for Environmental Scientists
Hoffimann, J and Zadrozny, B. 2019. Efficient variography with partition variograms

GeoStatsFunctions.DirectionalVariogram — Function

DirectionalVariogram(direction, data, var₁, var₂=var₁; dtol=1e-6, [parameters])

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

Optionally, forward parameters for the underlying EmpiricalVariogram.

GeoStatsFunctions.PlanarVariogram — Function

PlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6, [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.

Optionally, forward parameters for the underlying EmpiricalVariogram.

Base.values — Method

values(γ)

Returns the abscissa, the ordinate, and the bin counts of the empirical variogram γ.

GeoStatsFunctions.distance — Method

distance(γ)

Return the distance used to compute the empirical variogram γ.

GeoStatsFunctions.estimator — Method

estimator(γ)

Return the estimator used to compute the empirical variogram γ.

Base.merge — Method

merge(γα, γβ)

Merge the empirical variogram γα with the empirical variogram γβ assuming that both variograms have the same number of lags, distance and estimator.

Consider the following example image:

using GeoArtifacts

𝒟 = GeoArtifacts.image("Gaussian30x10")

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

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

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

Mke.plot(γ)

directional variograms along a specific direction:

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

Mke.plot(γₕ, vcolor = :maroon, hcolor = :maroon)
Mke.plot!(γᵥ)
Mke.current_figure()

or planar variograms over a specific plane:

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

Mke.plot(γₕ, vcolor = :maroon, hcolor = :maroon)
Mke.plot!(γᵥ)
Mke.current_figure()

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

GeoStatsFunctions.EmpiricalVarioplane — Type

EmpiricalVarioplane(data, var₁, var₂=var₁;
                    normal=spheredir(0,0),
                    nangs=50, ptol=0.5, dtol=0.5,
                    [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 for the plane partition, the tolerance dtol for the direction partition, the number of angles nangs in the plane, and forward the parameters to the underlying EmpiricalVariogram.

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

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

fig = Mke.Figure()
ax = Mke.PolarAxis(fig[1,1], title = "Varioplane")
Mke.plot!(ax, γ)
fig