Field processes

Builtin

The following field processes are available upon loading GeoStats.jl.

GeoStatsProcesses.GaussianProcess — Type

GaussianProcess(function, [mean])

Gaussian process with given geostatistical function (e.g. variogram) and mean (default to zero mean).

Examples

# univariate processes
GaussianProcess(GaussianVariogram())
GaussianProcess(SphericalCovariance(), 0.5)

# multivariate processes
GaussianProcess(LinearTransiogram(), [0.5, 0.5])

source

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

# Gaussian process
proc = GaussianProcess(GaussianVariogram(range=30.0))

# unconditional simulation
real = rand(proc, grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].field)
viz(fig[1,2], real[2].geometry, color = real[2].field)
fig

The GaussianProcess can be simulated with different methods. Heuristics are used to select the most appropriate method for the data and domain at hand.

GeoStatsProcesses.LUSIM — Type

LUSIM(; [options])

The LU simulation method introduced by Alabert 1987.

The full covariance matrix is built to include all locations of the data, and samples from the multivariate Gaussian are drawn via a lower-upper (LU) matrix factorization.

Options

factorization - Factorization function from LinearAlgebra (default to cholesky)

References

Alabert 1987. The practice of fast conditional simulations through the LU decomposition of the covariance matrix
Oliver 2003. Gaussian cosimulation: modeling of the cross-covariance

Notes

The method is only adequate for domains with relatively small number of elements (e.g. 100x100 grids) where it is feasible to factorize the full covariance.

For larger domains (e.g. 3D grids), other methods are preferred such as SEQSIM and FFTSIM.

source

GeoStatsProcesses.FFTSIM — Type

FFTSIM(; [options])

The FFT simulation method introduced by Gutjahr 1997.

The covariance function is perturbed in the frequency domain after a fast Fourier transform. White noise is added to the phase of the spectrum, and a realization is produced by an inverse Fourier transform.

Options

minneighbors - Minimum number of neighbors (default to 1)
maxneighbors - Maximum number of neighbors (default to 26)
neighborhood - Search neighborhood (default to nothing)
distance - Distance used to find nearest neighbors (default to Euclidean())

References

Gutjahr 1997. General joint conditional simulations using a fast Fourier transform method
Gómez-Hernández, J. & Srivastava, R. 2021. One Step at a Time: The Origins of Sequential Simulation and Beyond

Notes

The method is limited to simulations on regular grids, and care must be taken to make sure that the correlation length is small enough compared to the grid size. As a general rule of thumb, avoid correlation lengths greater than 1/3 of the grid.

Visual artifacts can appear near the boundaries of the grid if the correlation length is large compared to the grid itself.

source

GeoStatsProcesses.SEQSIM — Type

SEQSIM(; [options])

The sequential simulation method introduced by Gomez-Hernandez 1993.

The method traverses all locations of the geospatial domain according to a path, approximates the conditional distribution within a neighborhood using a geostatistical model, and assigns a value to the center of the neighborhood by sampling from this distribution.

Options

path - Process path (default to LinearPath())
minneighbors - Minimum number of neighbors (default to 1)
maxneighbors - Maximum number of neighbors (default to 26)
neighborhood - Search neighborhood (default to :range)
distance - Distance used to find nearest neighbors (default to Euclidean())

For each location in the process path, a maximum number of neighbors maxneighbors is used to fit the conditional Gaussian distribution. The neighbors are searched according to a neighborhood.

The neighborhood can be a MetricBall, the symbol :range or nothing. The symbol :range is converted to MetricBall(range(f)) where f is the geostatistical function of the Gaussian process. If neighborhood is nothing, the nearest neighbors are used without additional constraints.

References

Gomez-Hernandez & Journel 1993. Joint Sequential Simulation of MultiGaussian Fields

Notes

This method is very sensitive to the neighbor search options. Care must be taken to make sure that enough neighbors are used in the underlying geostatistical model.

source

GeoStatsProcesses.LindgrenProcess — Type

LindgrenProcess(range=1.0, sill=1.0; init=NearestInit())

Lindgren process with given range (correlation length) and sill (total variance) as described in Lindgren 2011.

The process relies relies on a discretization of the Laplace-Beltrami operator on meshes and is adequate for highly curved domains (e.g. surfaces).

References

Lindgren et al. 2011. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach

source

# domain of interest
mesh = simplexify(Sphere((0, 0, 0), 1))

# Lindgren process
proc = LindgrenProcess()

# unconditional simulation
real = rand(proc, mesh, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].field)
viz(fig[1,2], real[2].geometry, color = real[2].field)
fig

External

The following field processes are available upon loading external packages.

ImageQuilting.jl

GeoStatsProcesses.QuiltingProcess — Type

QuiltingProcess(trainimg, tilesize; [options])

Image quilting process with training image trainimg and tile size tilesize as described in Hoffimann et al. 2017.

Options

overlap - Overlap size (default to (1/6, 1/6, ..., 1/6))
path - Process path (:raster (default), :dilation, or :random)
inactive - Vector of inactive voxels (i.e. CartesianIndex) in the grid
soft - A pair (data,dataTI) of geospatial data objects (default to nothing)
tol - Initial relaxation tolerance in (0,1] (default to 0.1)

References

Hoffimann et al 2017. Stochastic simulation by image quilting of process-based geological models
Hoffimann et al 2015. Geostatistical modeling of evolving landscapes by means of image quilting

source

using ImageQuilting
using GeoStatsImages

# domain of interest
grid = CartesianGrid(200, 200)

# quilting process
img  = geostatsimage("Strebelle")
proc = QuiltingProcess(img, (62, 62))

# unconditional simulation
real = rand(proc, grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].facies)
viz(fig[1,2], real[2].geometry, color = real[2].facies)
fig

# domain of interest
grid = CartesianGrid(200, 200)

# quilting process
img  = geostatsimage("StoneWall")
proc = QuiltingProcess(img, (13, 13))

# unconditional simulation
real = rand(proc, grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].Z)
viz(fig[1,2], real[2].geometry, color = real[2].Z)
fig

# domain of interest
grid = domain(img)

# pre-existing observations
img  = geostatsimage("Strebelle")
data = img |> Sample(20, replace=false)

# quilting process
proc = QuiltingProcess(img, (30, 30))

# conditional simulation
real = rand(proc, grid, 2, data=data)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].facies)
viz(fig[1,2], real[2].geometry, color = real[2].facies)
fig

Voxels marked with the special symbol NaN are treated as inactive. The algorithm will skip tiles that only contain inactive voxels to save computation and will generate realizations that are consistent with the mask. This is particularly useful with complex 3D models that have large inactive portions.

# domain of interest
grid = domain(img)

# skip circle at the center
nx, ny = size(grid)
r = 100; circle = []
for i in 1:nx, j in 1:ny
  if (i-nx÷2)^2 + (j-ny÷2)^2 < r^2
    push!(circle, CartesianIndex(i, j))
  end
end

# quilting process
proc = QuiltingProcess(img, (62, 62), inactive = circle)

# unconditional simulation
real = rand(proc, grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].facies)
viz(fig[1,2], real[2].geometry, color = real[2].facies)
fig

It is possible to incorporate auxiliary variables to guide the selection of patterns from the training image.

using ImageFiltering

# image assumed as ground truth (unknown)
truth = geostatsimage("WalkerLakeTruth")

# training image with similar patterns
img = geostatsimage("WalkerLake")

# forward model (blur filter)
function forward(data)
  img = asarray(data, :Z)
  krn = KernelFactors.IIRGaussian([10,10])
  fwd = imfilter(img, krn)
  georef((; fwd=vec(fwd)), domain(data))
end

# apply forward model to both images
data   = forward(truth)
dataTI = forward(img)

proc = QuiltingProcess(img, (27, 27), soft=(data, dataTI))

real = rand(proc, domain(truth), 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].Z)
viz(fig[1,2], real[2].geometry, color = real[2].Z)
fig

TuringPatterns.jl

GeoStatsProcesses.TuringProcess — Type

TuringProcess(; [options])

Turing process as described in Turing 1952.

Options

params - basic parameters (default to nothing)
blur - blur algorithm (default to nothing)
edge - edge condition (default to nothing)
iter - number of iterations (default to 100)

References

Turing 1952. The chemical basis of morphogenesis

source

using TuringPatterns

# domain of interest
grid = CartesianGrid(200, 200)

# unconditional simulation
real = rand(TuringProcess(), grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].field)
viz(fig[1,2], real[2].geometry, color = real[2].field)
fig

StratiGraphics.jl

GeoStatsProcesses.StrataProcess — Type

StrataProcess(environment; [options])

Strata process with given geological environment as described in Hoffimann 2018.

Options

state - Initial geological state
stack - Stacking scheme (:erosional or :depositional)
nepochs - Number of epochs (default to 10)

References

Hoffimann 2018. Morphodynamic analysis and statistical synthesis of geormorphic data

source

using StratiGraphics

# domain of interest
grid = CartesianGrid(50, 50, 20)

# stratigraphic environment
p = SmoothingProcess()
T = [0.5 0.5; 0.5 0.5]
Δ = ExponentialDuration(1.0)
ℰ = Environment([p, p], T, Δ)

# strata simulation
real = rand(StrataProcess(ℰ), grid, 2)

fig = Mke.Figure(size = (800, 400))
viz(fig[1,1], real[1].geometry, color = real[1].field)
viz(fig[1,2], real[2].geometry, color = real[2].field)
fig