[SOLVED] Anisotropic variogram ranges and sill values

[ben_esp]

Hello all,

I am trying to get the ranges for the variograms in the following directions: 0, 45, 90 and 135 deg (anti-clockwise directions from East).

I am using the following code:

Code: Select all

rm(list=ls()) # clean global environment

library("RGeostats") # Loading 'RGeostats' library
library(sp) # Loading 'sp' library

SG = read.csv("export.csv") # Load the data
SG = SG[complete.cases(SG),] # Get rid of incomplete data

SG$export.log = log(SG$export) # calculate the log of export

datRG = db.create(SG) # Create the data base
datRG = db.locate(datRG, names=c("lon", "lat") ,"x")  # define coordinates
datRG = db.locate(datRG, names=c("export.log") ,"z")  # define working variables


v = vario.calc(datRG,lag=1,nlag=30,dirvect=c(0,45,90,135)) # experimental variogram with directions 0, 45, 90 and 135
plot(v,title="Directional variograms",npairdw=TRUE,npairpt=1, pos.legend = 1) # plot the experimental variograms
m <-  model.auto(v, struct=c("Nugget Effect",melem.name(c(3))),title="Directionnal variogram for export", pos.legend = 2) # fit a model
print(m) # print the model


This is the output I get when printing the model:

Model characteristics
=====================
Space dimension = 2
Number of variable(s) = 1
Number of basic structure(s) = 2
Number of drift function(s) = 1
Number of drift equation(s) = 1

Covariance Part
---------------
- Nugget Effect
Sill = 0.106
- Spherical
Sill = 0.188
Anisotropy = 1.000 0.631
Ranges = 2.385 1.504
Rotation =
[,1] [,2]
[1,] 0.999 0.053
[2,] -0.053 0.999
Angle = -3.022 degrees
Total Sill = 0.294

Drift Part
----------
Universality Condition

I am confused, because I have 4 model fits and I only see 2 ranges.

How can I get the range and sill for each direction?

When I do this:

Code: Select all

m[1]$range

I get:

8.628392

I don't know what this number is, it doesn't appear in the "print(m)" summary.

And when I do this:

Code: Select all

m[2]$range

I get:

2.384706

Which seems to correspond to the first range in the Spherical part. I can't reach the second value as m[2]$range[2] gives a 'NA'.

And when I do this:

Code: Select all

m[3]$range

I get:

Error in model@basics[[istr]] : subscript out of bounds

I might have misunderstood the documentations somehow.

Thank you in advance,

Benoit

[Didier Renard]

Hello Benoit

I think there is some misunderstanding on the Model in RGeostats. Your model is composed of two basic structures (i.e. a nugget effect and a spherical model).
Then you should not be surprised that the only available assessors are based on model[1] and model[2] (any other argument is erroneous returning "index out of bounds").
Then each basic structure can be isotropic or anisotropic.
Let us recall that the nugget effect is, by construction, isotropic... so as the answer of the assessor: model[1]$flag.aniso.
Conversely, the second basic structure (spherical) is anisotropic ... so as the answer of the assessor: model[2]$flag.aniso

Let us now concentrate on the second basic structure (i.e. spherical component).
You can retrieve a range (model[2]$range): the returned value is 2.385, whereas the range, as printed in the model, appears to be a vector: 2.385 and 1.504.
In fact, the printed result is a combination of the UNIQUE range value (2.385) and the vector of anisotropy coefficients: model[2]$aniso.coeffs = 1. and 0.630608.
Note that the value 1.504 = 0.630608 * 2.385

You understood that, as defined in geostatistical course, the second basic structure follows a geometrical anisotropy: the range follows an ellipse, with axes equal to 2.385 and 1.504; the longest axis is rotated at -3.022 degrees.

Now your question is: how can be extract the value of the apparent range in any direction defined in 2D.
Here is the small piece of program that answers to your question.
I have defined a small function (named e) which takes the model and the rank of the basic structure in input, as well as the target angle and returns the value of this apparent range:

Code: Select all

e <- function(model,is=1,angle=0)
{
   if (model$ndim != 2) stop("Only coded in 2-D")
   ndim  = model$ndim
   melem = model[is]   
   range = melem$range
   if (melem$flag.aniso)
   {
      rotmat = melem$aniso.rotmat
      coeffs = melem$aniso.coeffs
      vect   = c(cosd(angle) * range * coeffs[1],
               sind(angle) * range * coeffs[2])
      u = rotmat %*% vect
      range = sqrt(u[1] * u[1] + u[2] * u[2])
   }   
   range
}


As a demonstration, we can ask the value of the apparent range:
- e(model,2,-3.022) which obviously returns the valus of the longest axis, i.e. 2.383 (rounded value for 2.385)
- e(model,2,-3.022+90) returns 1.507163 (instead of 1.504)

I have even imagined a function "ep" calling the internal "e" function for any angle ranging from 0 to 360 degrees.

Code: Select all

ep <- function(model,is)
{
   angles = seq(0,360)
   number = length(angles)
   tab = numeric(number)
   for (i in 1:number)
      tab[i] = e(model,is,angles[i])
   print(range(tab))
   plot(angles,tab)
}


This function shows the shape of the apparent range as a function of the angle. The range of ranges is: 1.504 to 2.385.

Hope this will help

[ben_esp]

Hello Didier,

Thanks, this has helped a lot.

I have noticed that if a model has two directions (0 and 90), you don't get an angle. Is that because I only used two directions or is that related to my direction choices?

Also, when you have the long and short axis, are they the major and minor axes of the ellipse or the semi-major and semi-minor axes?

Thank you,

Benoit

[Didier Renard]

As you noticed, model.auto() cannot find a direction when only 2 directions are calculated (in a 2-D problem). This is an obvious result. Therefore model.auto() considers that the 2 directions in which the calculation of the experimental variograms have been performed, should correspond to the major and minor axes of the anisotropy ellipse.
The function model.auto() then concentrates on the anisotropy ranges.
When the experimental variogram is calculated in more than 2 directions, then model.auto() can determine the angle as well as the ranges in the minor and major axes of the anisotropy ellipse. These two ranges are not systematically ordered (shorter before longer or longer before shorter). The only rule is that the ranges are provided along the U, then V axes of the rotated ellipse. And the angle corresponds to the rotation from the X axis to the U axis.

[ben_esp]

Hello Didier,

Sorry I did not answer you earlier, I was in a week-long trip and just got back.

Thank you very much for your answers. They are valuable.

All the best,

Benoit