Local anisotropies and nonstationary spatial processes for the GeoStats.jl framework
Author rmcaixeta
6 Stars
Updated Last
5 Months Ago
Started In
September 2020


Build Status Coverage

Local anisotropies ellipses extracted from reference data


This package deals with local anisotropies in geostatistics. It offers some solutions to extract them from a reference input. Then, it's possible to use the local anisotropies to create nonstationary spatial models, which can be used into estimation methods such as kriging. There are also some extra tools, like anisotropy interpolation and conversion to/between different rotation conventions. It is designed for 2-D and 3-D data and is developed to be used as an extension of GeoStats.jl when dealing with nonstationarity of second order moments. A list of current implementations:

  • Local anisotropies extraction methods:
    • Gradients
    • SVD / PCA

  • Local anisotropies operations:
    • Interpolation of directional data (via quaternions)
    • Import/export/convert to different conventions
    • Unconventional local variography using local anisotropies

  • Nonstationary spatial methods:
    • Moving windows
    • Kernel convolution
    • Spatial deformation

  • Solvers adapted to these nonstationary spatial methods:
    • Kriging

Warning: This package is still under (slow) development. Some of the implementations were not fully validated and may give inconsistent results. Debugging, contributions and suggestions are welcome.


First, it is necessary to install Julia. Installation instructions for Windows, Linux and macOS are available here.

To install the package: open the Julia REPL and then install the package with the following command. The GeoStats.jl, GeoStatsPlots and Plots.jl packages are installed together to run the usage example.

using Pkg; Pkg.add(["LocalAnisotropies","GeoStats","GeoStatsPlots","Plots"])


Introduction to local anisotropies extraction methods:

Lillah & Boisvert (2015) Inference of locally varying anisotropy fields from diverse data sources

Introduction to nonstationary spatial methods:

Sampson (2010) Constructions for Nonstationary Spatial Processes

Moving windows:

Haas (1990) Kriging and automated variogram modeling within a moving window
Stroet & Snepvangers (2005) Mapping curvilinear structures with local anisotropy kriging

Kernel convolution:

Higdon (1998) A process-convolution approach to modelling temperatures in the North Atlantic Ocean
Fouedjio et al. (2016) A generalized convolution model and estimation for non-stationary random functions

Spatial deformation:

Sampson & Guttorp (1992) Nonparametric estimation of nonstationary spatial covariance structure
Boisvert (2010) Geostatistics with locally varying anisotropy


The documentation of the main functions are available as docstrings Check below an usage example that illustrate the package applications.

Usage example

# load libraries for the example
using LocalAnisotropies
using GeoStats
using GeoStatsPlots
using Plots
using Random

# create a reference scenario for tests
D = georef((P=[25-abs(0.2*i^2-j) for i in -10:9, j in 1:20],))
S = sample(D, 80, replace=false)
S = georef(values(S), PointSet(centroid.(domain(S))))
G = CartesianGrid(20,20)

# create an estimation problem
P = EstimationProblem(S, G, :P)

# plot reference scenario and samples extracted for further estimations
splot = plot(G)
splot = plot!(S)

# get local anisotropies using gradients of a 8x8 radius window
rawlpars = localanisotropies(Gradients, D, :P, 8)
plot(D, alpha=0.6, colorbar=false)

# average 10 nearest local anisotropies
searcher = KNearestSearch(G, 10)
lpars = smoothpars(rawlpars, searcher)

# rescale magnitude of range2 to between 0.25 and 1.0
rescale_magnitude!(lpars, r2=(0.25,1.0))

# set possible reference axes to X and Y; ranges will be fixed in those directions
lparsx = reference_magnitude(lpars, :X)
lparsy = reference_magnitude(lpars, :Y)

plot(D, alpha=0.6, colorbar=false)

# plot(lpars, data) will only work for 2D spatial data
# for 3D or custom visualizations, it's possible to export it to VTK
to_vtk("ellipses", D, lpars)
# below the file "ellipses.vtu" loaded in Paraview using TensorGlyph (disable extract eigenvalues)

# pass local anisotropies to samples
spars = nnpars(lpars, D, S)

# do an unconventional variography along local X axis (same can be done for Y)
expvario = localvariography(S, spars, :P, tol=2, maxlag=20, nlags=20, axis=:X)
γx = ExponentialVariogram(sill=32., range=40.)
γy = GaussianVariogram(sill=32., range=8.)

# kriging using moving windows method
MW = LocalKriging(:P => (variogram=γx, localaniso=lparsx, method=:MovingWindows))
s1 = solve(P, MW)

# kriging using kernel convolution method (smaller search; unstable with many local variations)
KC = LocalKriging(:P => (variogram=γy, localaniso=lparsy, method=:KernelConvolution, maxneighbors=6))
s2 = solve(P, KC)

# deform space using kernel variogram as dissimilarity input
Sd1, Dd1 = deformspace(S, G, lparsx, KernelVariogram, γx, anchors=1500)
Pd1 = EstimationProblem(Sd1, Dd1, :P)
γ1 = GaussianVariogram(sill=21.3, range=22.5)
s3 = solve(Pd1, Kriging(:P => (variogram=γ1,)))
plot(plot(to_3d(s3),[:P]), plot(georef(values(s3),G),[:P],colorbar=false))

# deform space based on a graph built with average anisotropic distances
# of the ten nearest data; do variography in multidimensional space
LDa = graph(S, G, lparsx, AnisoDistance, searcher)
Sd2, Dd2 = deformspace(LDa, GraphDistance, anchors=1500)
γ2 = GaussianVariogram(sill=22., range=22.)

# traditional kriging in the new multidimensional space
Pd2 = EstimationProblem(Sd2, Dd2, :P)
s4 = solve(Pd2, Kriging(:P => (variogram=γ2,)))
plot(plot(to_3d(s4),[:P]), plot(georef(values(s4),G),[:P],colorbar=false))

# deform space based on a graph built with kernel variogram of the ten
# nearest data; do variography in multidimensional space
LDv = graph(S, G, lparsy, KernelVariogram, γy, searcher)
Sd3, Dd3 = deformspace(LDv, GraphDistance, anchors=1500)
γ3 = NuggetEffect(1.0) + GaussianVariogram(sill=21., range=22.)

# traditional kriging in the new multidimensional space
Pd3 = EstimationProblem(Sd3, Dd3, :P)
s5 = solve(Pd3, Kriging(:P => (variogram=γ3,)))
plot(plot(to_3d(s5),[:P]), plot(georef(values(s5),G),[:P],colorbar=false))

γomni = GaussianVariogram(sill=32., range=11.)
OK = Kriging(:P => (variogram=γomni, maxneighbors=20))
s0 = solve(P, OK)

# comparison of the different estimates
mse(a,b) = sum((a .- b) .^ 2)/length(b)
solvers = ["OK","MW","KC","SD1","SD2","SD3"]
errors  = [mse(getproperty(x,:P),getproperty(D,:P)) for x in [s0,s1,s2,s3,s4,s5]]
bar(solvers,errors,legend=false,ylabel="Mean squared error",xlabel="Estimation method")

Some extra tools to work with local anisotropies:

# import external local anisotropies in GSLIB convention
dummy = georef((az=1:10, r1=1:10, r2=1:10), PointSet(rand(2,10)))
pars  = localanisotropies(dummy, [:az], [:r1,:r2], :GSLIB)

# interpolate local anisotropies into a coarser grid
G_ = CartesianGrid((10,10),(0.5,0.5),(2.0,2.0))
lpars_ = idwpars(lpars, searcher, G_, power=2.0)

# convert between different rotation conventions
angs1 = convertangles([30,30,30], :GSLIB, :Datamine)
angs2 = convertangles.(pars.rotation, :GSLIB)