Edit on GitHub

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

  • pool - Pool of worker processes (default to [myid()])
  • threads - Number of threads (default to cpucores())
  • verbose - Show progress and other information (default to true)

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])

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

GeoStatsBase.Ensemble — Type
Ensemble

An ensemble of geostatistical realizations from a geostatistical process.

# 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{2,Float64}
  variables: Z (Float64)
  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

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

# ------------
# main script
# ------------

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

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

# generate three realizations with three processes
real = rand(proc, grid, [:Z => Float64], 3, pool = workers())

Please consult The ultimate guide to distributed computing in Julia.

Builtin

The following processes are shipped with the framework.

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

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.

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

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.

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

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

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

using ImageQuilting
using GeoArtifacts

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

# quilting process
img  = GeoArtifacts.image("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  = GeoArtifacts.image("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  = GeoArtifacts.image("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 = GeoArtifacts.image("WalkerLakeTruth")

# training image with similar patterns
img = GeoArtifacts.image("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

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

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