# Empirical variograms

An empirical variogram has the form:

$$$\hat{\gamma}(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. Empirical variograms can be estimated using general distance functions. These can be used in order to for example:

• Model anisotropy (e.g. ellipsoid distance)
• Perform geostatistical simulation on spherical coordinate systems (e.g. haversine distance)

Please see Distances.jl for a complete list of options.

Additionally, given two empirical variograms $\hat{\gamma}_\alpha$ and $\hat{\gamma}_\beta$, they can be merged as described in Hoffimann & Zadrozny 2019:

$$$\hat{\gamma}_{\alpha+\beta}(h) = \frac{|N_\alpha(h)| \cdot \hat{\gamma}_\alpha(h) + |N_\beta(h)| \cdot \hat{\gamma}_\beta(h)}{|N_\alpha(h)| + |N_\beta(h)|}$$$
Base.mergeMethod
merge(γα, γβ)

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

## Variograms

Consider the following image for illustration purposes:

using GeoStatsImages

𝒟 = geostatsimage("Gaussian30x10")

plot(𝒟)
Variography.EmpiricalVariogramType
EmpiricalVariogram(data, var₁, var₂=var₁; [optional parameters])

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

Parameters

• nlags - number of lags (default to 20)
• maxlag - maximum lag (default to half of maximum lag of data)
• distance - custom distance function (default to Euclidean distance)
• algo - accumulation algorithm (default to :ball)

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 data.

The function values can be used to retrieve the abscissa, ordinate and bin count of an empirical variogram:

julia> x, y, n = values(γ)

See also: DirectionalVariogram

References

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

plot(γ)
Variography.DirectionalVariogramFunction
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.

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

plot(γₕ, label="horizontal")
plot!(γᵥ, label="vertical")
Variography.PlanarVariogramFunction
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.

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

plot(γₕ, label="horizontal")
plot!(γᵥ, label="vertical")

## Varioplanes

Variography.EmpiricalVarioplaneType
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.

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

plot(γ)