Fields

Base.rand — Method

rand([rng], process::FieldProcess, domain, data, [nreals], [method]; [parameters])

Generate one or nreals realizations of the field process with method over the domain with data and optional paramaters. Optionally, specify the random number generator rng and global parameters.

The data can be a geotable, a pair, or an iterable of pairs of the form var => T, where var is a symbol or string with the variable name and T is the corresponding data type.

Parameters

workers - Worker processes (default to workers())
threads - Number of threads (default to cpucores())
verbose - Show progress and other information (default to true)
async - Return to master process immediately (default to false)

Examples

julia> rand(process, domain, geotable, 3)
julia> rand(process, domain, :z => Float64)
julia> rand(process, domain, "z" => Float64)
julia> rand(process, domain, [:a => Float64, :b => Float64])
julia> rand(process, domain, ["a" => Float64, "b" => Float64])

source

Realizations are stored in an Ensemble as illustrated in the following example:

GeoStatsProcesses.Ensemble — Type

Ensemble

An ensemble of geostatistical realizations from a geostatistical process.

source

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

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

# unconditional simulation
real = rand(proc, grid, [:Z => Float64], 100)

2D Ensemble
  domain:    100×100 CartesianGrid
  variables: Z
  N° reals:  100

We can visualize the first two realizations:

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

the mean and variance:

m, v = mean(real), var(real)

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

or the 25th and 75th percentiles:

q25 = quantile(real, 0.25)
q75 = quantile(real, 0.75)

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

The cummulative distribution function is obtained likewise:

p25 = cdf(real, 0.25)
p75 = cdf(real, 0.75)

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

All field processes can generate realizations in parallel using multiple Julia processes. Doing so requires using the Distributed standard library, like in the following example:

using Distributed

# request additional processes
addprocs(3)

# load code on every single process
@everywhere using GeoStats

# setup simulation
grid = CartesianGrid(100, 100)
proc = GaussianProcess(GaussianVariogram(range=30.0))

# perform simulation on all worker processes
real = rand(proc, grid, [:Z => Float64], 3, workers = workers())

Please consult The ultimate guide to distributed computing in Julia.

Builtin

The following processes are shipped with the framework.

GeoStatsProcesses.GaussianProcess — Type

GaussianProcess(variogram=GaussianVariogram(), mean=0.0)

Gaussian process with given theoretical variogram and global mean.

source

GeoStatsProcesses.LUMethod — Type

LUMethod(; [paramaters])

The LU Gaussian process method introduced by Alabert 1987. The full covariance matrix is built to include all locations of the process domain, and samples from the multivariate Gaussian are drawn via LU factorization.

Parameters

factorization - Factorization method (default to cholesky)
correlation - Correlation coefficient between two covariates (default to 0)
init - Data initialization method (default to NearestInit())

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 SEQMethod and FFTMethod.

source

GeoStatsProcesses.FFTMethod — Type

FFTMethod(; [paramaters])

The FFT Gaussian process 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.

Parameters

minneighbors - Minimum number of neighbors (default to 1)
maxneighbors - Maximum number of neighbors (default to 10)
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 Cartesian 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.
The method is extremely fast, and can be used to generate large 3D realizations.

source

GeoStatsProcesses.SEQMethod — Type

SEQMethod(; [paramaters])

The sequential process method introduced by Gomez-Hernandez 1993. It traverses all locations of the geospatial domain according to a path, approximates the conditional Gaussian distribution within a neighborhood using simple Kriging, and assigns a value to the center of the neighborhood by sampling from this distribution.

Parameters

path - Process path (default to LinearPath())
minneighbors - Minimum number of neighbors (default to 1)
maxneighbors - Maximum number of neighbors (default to 36)
neighborhood - Search neighborhood (default to :range)
distance - Distance used to find nearest neighbors (default to Euclidean())
init - Data initialization method (default to NearestInit())

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(γ)) where γ is the variogram 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 various parameters. Care must be taken to make sure that enough neighbors are used in the underlying Kriging model.

source

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

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

# unconditional simulation
real = rand(proc, grid, [:Z => Float64], 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

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. Optionally, specify the data initialization method init.

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, [:Z => Float64], 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

External

The following processes are available in external packages.

ImageQuilting.jl

GeoStatsProcesses.QuiltingProcess — Type

QuiltingProcess(trainimg, tilesize; [paramaters])

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

Parameters

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)
init - Data initialization method (default to NearestInit())

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, [:facies => Int], 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, [:Z => Int], 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, data, 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

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, [:facies => Float64], 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), [:Z => Float64], 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(; [paramaters])

Turing process as described in Turing 1952.

Parameters

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, [:z => Float64], 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

StratiGraphics.jl

GeoStatsProcesses.StrataProcess — Type

StrataProcess(environment; [paramaters])

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

Parameters

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, [:z => Float64], 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