anis_GSLIBpar2A: Produce anisotropy scaling matrix from angle and anisotropy...

In gmGeostats: Geostatistics for Compositional Analysis

Produce anisotropy scaling matrix from angle and anisotropy ratios

Description

Produce anisotropy matrix (as the transposed of the Cholesky decomposition) from angle and anisotropy ratios

Usage

anis_GSLIBpar2A(ratios, angles, inv = FALSE)

anis2D_par2A(ratio, angle, inv = FALSE)

anis3D_par2A(ratios, angles, inv = FALSE)

Arguments

ratios	vector of two values between 0 and 1 giving the anisotropy ratios of medium/largest smallest/largest ranges
angles	as defined in gstat::vgm (and indeed GSLIB). For anis2D_par2A 'angle' is the direction of maximum range, i.e. largest spatial continuity, measured clockwise from North
inv	boolean or integer, see return for details
ratio	an anisotropy ratio (min/max range)
angle	direction of maximum range, i.e. largest spatial continuity, measured clockwise from North

Value

a 3x3 matrix of anisotropy.

If inv=TRUE (or 1) the output is a matrix A such that norm(h %*% A) allows to use isotropic variograms, being h = c(hx, hy, hz) the lag vectors.

If inv=FALSE (or 0) the output is a matrix A such that norm(h %*% solve(A)) allows to use isotropic variograms.

Other values are meaningless.

Functions

• anis2D_par2A: 2D case

• anis3D_par2A: 3D case

See Also

Other anisotropy: AnisotropyRangeMatrix(), AnisotropyScaling(), as.AnisotropyRangeMatrix(), as.AnisotropyScaling(), is.anisotropySpecification()

Examples

## ratio=0.5, azimuth 30?? (i.e. direction 60??)
A = anis2D_par2A(1, 30)
A
AAt = A %*% t(A)
 #  project the bisector 1:1 (i.e. 45??)
(k = c(1,1,0) %*% A)
atan2(k[2], k[1]) * 180/pi  # should be 15
sqrt(sum(k^2))
sqrt( c(1,1,0) %*% AAt %*% c(1,1,0) )
A = anis2D_par2A(0.5, 60)
rd = 60 * pi/180
A
A %*% t(A)
c(cos(rd), sin(rd),0) %*% A #  should be 1
c(-sin(rd), cos(rd),0) %*% A #  should be +/- sqrt(2)
c60 = cos(60*pi/180)
s60 = sin(60*pi/180)
c30 = cos(30*pi/180)
s30 = sin(30*pi/180)
#  in the new coordinates, 60cwN is (0,1,0)
R60p = anis3D_par2A(ratios=c(1,1), angles=c(60,0,0))
c(s60, c60, 0) %*% R60p
R6030 = anis3D_par2A(ratios=c(1,1), angles=c(60,30,0))
# the original X axis is positive on newX and newY, but negative on newZ
c(1,0,0) %*% R6030
# rotate first direction 60 degrees azimuth, then dip 30degrees upwards
c( c60*c30, -s60*c30, s30) %*% R6030
(Ranis = anis3D_par2A(ratios=c(0.5,0.25), angles=c(60,30,0)) )

gmGeostats documentation built on April 18, 2023, 5:08 p.m.