Fitting variograms

Overview

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.fitFunction
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")

varioplot(g)
Example block output

We can fit specific models to the empirical variogram:

γ = GeoStatsFunctions.fit(SineHoleVariogram, g)

varioplot(g)
varioplot!(γ, maxlag = 25u"m")
Mke.current_figure()
Example block output

or let the framework find the model with minimum error:

γ = GeoStatsFunctions.fit(Variogram, g)

varioplot(g)
varioplot!(γ, maxlag = 25u"m")
Mke.current_figure()
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)

varioplot(g)
varioplot!(γ, maxlag = 25u"m")
Mke.current_figure()
Example block output

Methods

Weighted least squares

GeoStatsFunctions.WeightedLeastSquaresType
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