LocalAnisotropies.jl

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

LocalAnisotropies

Build Status Coverage

Local anisotropies ellipses extracted from reference data

Introduction

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.

Installation

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

References

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

Documentation

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
Random.seed!(1234)

# 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)
plot(plot(D),splot)

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

# 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,D)

# 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)
plot(expvario)
ฮณx = ExponentialVariogram(sill=32., range=40.)
ฮณy = GaussianVariogram(sill=32., range=8.)
plot!(ฮณx)

# kriging using moving windows method
MW = LocalKriging(:P => (variogram=ฮณx, localaniso=lparsx, method=:MovingWindows))
s1 = solve(P, MW)
plot(s1,[:P])

# 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)
plot(s2,[:P])

# 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)
plot(s0,[:P])

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