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, Makie.jl and CairoMakie.jl packages are installed together to run the usage example.
using Pkg; Pkg.add(["LocalAnisotropies","GeoStats","Makie","CairoMakie"])Lillah & Boisvert (2015) Inference of locally varying anisotropy fields from diverse data sources
Sampson (2010) Constructions for Nonstationary Spatial Processes
Haas (1990) Kriging and automated variogram modeling within a moving window
Stroet & Snepvangers (2005) Mapping curvilinear structures with local anisotropy kriging
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
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.
# load libraries for the example
using LocalAnisotropies
using GeoStats
using Random
import CairoMakie as Mke
Random.seed!(1234)
Mke.set_theme!(size=(500, 500))
# 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, UniformSampling(80, replace=false))
S = georef(values(S), PointSet(centroid.(domain(S))))
G = CartesianGrid(20,20)
# plot reference scenario and samples extracted for further estimations
fig = Mke.Figure(size=(700, 350))
Mke.plot(fig[1,1],D.geometry,color=D.P)
Mke.plot(fig[1,2],G,showsegments=true,color=:white)
Mke.plot!(fig[1,2],S.geometry,color=S.P,pointsize=10)
Mke.current_figure()
# get local anisotropies using gradients of a 8x8 radius window
rawlpars = localanisotropies(Gradients, D, :P, 8)
Mke.plot(D.geometry, color=D.P, alpha=0.6)
Mke.plot!(rawlpars,D.geometry)
Mke.current_figure()
# 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)
Mke.plot(D.geometry, color=D.P, alpha=0.6)
Mke.plot!(lpars,D.geometry)
Mke.current_figure()
# for 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)
Mke.plot(expvario)
γx = ExponentialVariogram(sill=32., range=40.)
γy = GaussianVariogram(sill=32., range=8.)
Mke.plot!(γx)
Mke.current_figure()
# local inverse distance weighting using exponent 3
ID = LocalIDW(3, lparsx)
s1 = S |> LocalInterpolate(G, :P=>ID, maxneighbors=20)
Mke.plot(s1.geometry,color=s1.P)
# kriging using moving windows method
MW = LocalKriging(:MovingWindows, lparsx, γx)
s2 = S |> LocalInterpolate(G, :P=>MW, maxneighbors=20)
Mke.plot(s2.geometry,color=s2.P)
# kriging using kernel convolution method (smaller search; unstable with many local variations)
KC = LocalKriging(:KernelConvolution, lparsy, γy)
s3 = S |> LocalInterpolate(G, :P=>KC, maxneighbors=6)
Mke.plot(s3.geometry,color=s3.P)
# deform space using kernel variogram as dissimilarity input
Sd1, Dd1 = deformspace(S, G, lparsx, KernelVariogram, γx, anchors=1500)
γ1 = GaussianVariogram(sill=21.3, range=22.5)
s4 = Sd1 |> Interpolate(Dd1, :P=>Kriging(γ1))
# plot
fig4 = Mke.Figure(size=(700, 350))
Mke.plot(fig4[1,1],to_3d(s4).geometry,color=s4.P)
Mke.plot(fig4[1,2],G,color=s4.P)
Mke.current_figure()
# 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
s5 = Sd2 |> Interpolate(Dd2, :P=>Kriging(γ2))
# plot
fig5 = Mke.Figure(size=(700, 350))
Mke.plot(fig5[1,1],to_3d(s5).geometry,color=s5.P)
Mke.plot(fig5[1,2],G,color=s5.P)
Mke.current_figure()
# 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
s6 = Sd3 |> Interpolate(Dd3, :P=>Kriging(γ3))
# plot
fig6 = Mke.Figure(size=(700, 350))
Mke.plot(fig6[1,1],to_3d(s6).geometry,color=s6.P)
Mke.plot(fig6[1,2],G,color=s6.P)
Mke.current_figure()
γomni = GaussianVariogram(sill=32., range=11.)
OK = Kriging(γomni)
s0 = S |> InterpolateNeighbors(G, :P=>OK, maxneighbors=20)
Mke.plot(s0.geometry,color=s0.P)
# comparison of the different estimates
mse(a,b) = sum((a .- b) .^ 2)/length(b)
solvers = ["OK","LIDW","MW","KC","SD1","SD2","SD3"]
errors = [mse(getproperty(x,:P),getproperty(D,:P)) for x in [s0,s1,s2,s3,s4,s5,s6]]
Mke.barplot(1:7,errors,axis=(xticks=(1:7,solvers),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))
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)# normal score transformation
ns = Quantile()
S_ns, ref_S_ns = apply(ns, S)
Sd1_ns, ref_Sd1_ns = apply(ns, Sd1)
# normal scores unconventional variography along local X axis
expvario_ns = localvariography(S_ns, spars, :P, tol=2, maxlag=20, nlags=20, axis=:X)
γx_ns = SphericalVariogram(sill=1.0, range=15.)
# normal scores omni varios
expomni1_ns = EmpiricalVariogram(S_ns, :P, maxlag=20, nlags=20)
γomni1_ns = SphericalVariogram(sill=1.0, range=11.)
expomni2_ns = EmpiricalVariogram(Sd1_ns, :P, maxlag=20, nlags=20)
γomni2_ns = ExponentialVariogram(sill=1.0, range=50.)
figv = Mke.Figure(size=(700, 235))
Mke.plot(figv[1,1],expvario_ns)
Mke.plot!(figv[1,1],γx_ns,color=:red)
Mke.plot(figv[1,2],expomni1_ns)
Mke.plot!(figv[1,2],γomni1_ns,color=:red)
Mke.plot(figv[1,3],expomni2_ns)
Mke.plot!(figv[1,3],γomni2_ns,maxlag=20,color=:red)
Mke.current_figure()
# standard sequential gaussian simulation with isotropic variogram
sgs = SEQMethod(maxneighbors=25)
sims1 = rand(GaussianProcess(γomni1_ns), G, S_ns, 100, sgs)
med1 = quantile(sims1, 0.5)
sims1_bt = [revert(ns, x, ref_S_ns) for x in (sims1[1],med1)]
fig7 = Mke.Figure(size=(700, 350))
Mke.plot(fig7[1,1],G,color=sims1_bt[1].P,colormap=:jet)
Mke.plot(fig7[1,2],G,color=sims1_bt[2].P,colormap=:jet)
Mke.current_figure()
# sequential gaussian simulation with local anisotropies (MW method)
# not theoretically very correct, but can give good results
local_sgs = LocalSGS(lparsx,maxneighbors=25)
sims2 = rand(GaussianProcess(γx_ns), G, S_ns, 100, local_sgs)
med2 = quantile(sims2, 0.5)
sims2_bt = [revert(ns, x, ref_S_ns) for x in (sims2[1],med2)]
fig8 = Mke.Figure(size=(700, 350))
Mke.plot(fig8[1,1],G,color=sims2_bt[1].P,colormap=:jet)
Mke.plot(fig8[1,2],G,color=sims2_bt[2].P,colormap=:jet)
Mke.current_figure()
# sequential gaussian simulation after spatial deformation
# standard simulation after local anisotropies is "removed"
sims3 = rand(GaussianProcess(γomni2_ns), Dd1, Sd1_ns, 100, sgs)
med3 = quantile(sims3, 0.5)
sims3_bt = [revert(ns, x, ref_Sd1_ns) for x in (sims3[1],med3)]
fig9 = Mke.Figure(size=(700, 350))
Mke.plot(fig9[1,1],G,color=sims3_bt[1].P,colormap=:jet)
Mke.plot(fig9[1,2],G,color=sims3_bt[2].P,colormap=:jet)
Mke.current_figure()
# comparison of the different simulation methods
sim_solvers = ["SGS","SGS_MW","SGS_SD1"]
sim_errors = [mse(getproperty(x,:P),getproperty(D,:P)) for x in [sims1_bt[2],sims2_bt[2],sims3_bt[2]]]
fig10 = Mke.Figure(size=(700, 350))
Mke.barplot(fig10[1,1],1:3,sim_errors,axis=(xticks=(1:3,sim_solvers),ylabel="Mean squared error",xlabel="Simulation method (median)"))
Mke.current_figure()
