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