Geostatistics - Spatial Correlation and Variograms

Milkias Z. SEMEREAB

GEOL0097-Geostatistics, University of Liège

GitHub Link

October 31, 2022

---

• Task 1: Structural Analysis

Before we dive into variogram and spatial modelling, we can use the EDA tools that we know already to visualize and interpret the distribution and spatial continuity of the variables (heavy metals concentration) in the X, Y and Z directions.

# distribution of Zn values through x-direction?
plot(data$logZn ~ data$x, xlab="Easting", ylab="Log-Zn concentration")
abline(lm(data$logZn ~ data$x), col="red")

Note:
- abline() = add Straight line to a an existing Plot
- lm() = fits Linear Models (e.g. linear regression)

#  distribution of Zn values through y-direction?
plot(data$logZn ~ data$y, xlab="Northing", ylab="Log-Zn concentration")
abline(lm(data$logZn ~ data$y), col="red")

# distribution of Zn values through z-direction?
plot(data$logZn ~ data$elev, xlab="Elevation", ylab="Log-Zn concentration")
abline(lm(data$logZn ~ data$elev), col="red")

Questions:

• Discuss about the distribution/spatial continuity of Zn values through x, y and z directions?

• Task 2: Introducing RGeostats

• RGeostats package:

For the variography analysis we will use the RGeostats package.RGeostats is the Geostatistical R Package developed by the Geostatistics Lab of the Geosciences Research Center of MINES ParisTech. RGeostats offers most of the famous geostatistical techniques and it is designed to tackle problems with several variables defined in a space of any dimensions. It is particularly adapted for students or researchers who want to test geostatistical procedures using R scripts. Link to RGeostats website.

To instal RGeostats package: - Install the package Rcpp (an interface between R and C++ langauges that is required when using RGeostats). It is in CRAN archieve (R’s central packages repository) and you can easily install it just by writing its name in the package window in Rstudio (like we did for corrplot and plot3D).

• On the other hand, RGeostats is not yet published to CRAN and we need to download and install it manually. Download RGeostats here.

• Once you download the zip folder, Install RGeostats manually: Click Packages Tab → Install → Select Package Archive File (.zip, .tar.gz) in the Install from slot → Find the corresponding zipped file on your local PC, and click Open → Click Install

• You can check if it is successfully installed by going to the packages window and see if Rgeostats is in the list of installed packages.

• Load/import Rgeostats package by checking the box or using the library() function.

# Load Rgeostats package
library(RGeostats)

• Create standard database

First we need to turn our data (which is in data.frame class) into a database (DB class) using the function db.create(). Once database is created, you can Check its contents using the function print() or db.print().

We will focus on the log-transformed Zn concentrations. The same analysis can be repeated for other heavy metal concentrations (Cd, Cu, and Pb).

# create a standard DB and store it in a variable called db.data
db.data = db.create(x1=data[,1], x2=data[,2], z1=data[,17],flag.grid=FALSE)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

Where:

• x1=data[,1] is Easting (m)
• x2=data[,2] is Northing (m)
• z1=data[,17] is LogZn concentration values (ppm)

Remember: Help window and help function are alaways at your disposal.

# use help() function and ? help operator to access the documentation pages for R functions
help("db.create")
?db.create

#check the content of the database you just created
class(db.data)

[1] "db"
attr(,"package")
[1] "RGeostats"

db.print(db.data)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

Data Base Characteristics
=========================

Data Base Summary
-----------------
File is organized as a set of isolated points
Space dimension              = 2
Number of Columns            = 4
Maximum Number of UIDs       = 4
Total number of samples      = 155

Variables
---------
Column = 0 - Name = rank - Locator = NA
Column = 1 - Name = x1 - Locator = x1
Column = 2 - Name = x2 - Locator = x2
Column = 3 - Name = z1 - Locator = z1

• The first column in the database is a ranking column which contains the value that is the rank of the particular data point where 1 is assigned to the first data point and 155 is assigned to the highest data point.

• Second and third columns are related to the spatial information of the samples (easting and northing)

• The third corresponds to the attribute being analyzed (logZn).

• Graphic representation of the data (proportional representation) using the function db.plot() or plot() :

# using db.plot() function
db.plot(db.data,xlab="Easting(m)", ylab = "Northing(m)", title="LogZn concentration")

# using the usual plot() function
plot(db.data, xlab="Easting(m)", ylab = "Northing(m)", title="LogZn concentration")

Both functions produced the same plot as we did last time but using a different graphic representation.

• Task 3: Omni-Directional Variography Analysis

• Compute and plot omni-directional experimental variograms of various number of lags using functions of RGeostats package

# Define the Diagonal Geographic Domain
D = sqrt((max(data[,1])-min(data[,1]))^2 + (max(data[,2])-min(data[,2]))^2) 

# Display value of D
print(D)

[1] 4789.868

Note: lag distance should not exceed half of longest dimension in the field.

# Define number of lags. It is good practise to start with 10 lags
nlag = 10

# Compute the omni-directional experimental variogram with 10 lags using  vario.calc() function and store it in a variable called Vario_E1. 
Vario_E1 = vario.calc(db.data,nlag=nlag,lag=D/(2*nlag))

# display the variogram characteristics 
print(Vario_E1)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

Variogram characteristics
=========================
Number of variable(s)       = 1
Number of direction(s)      = 1
Space dimension             = 2
Variance-Covariance Matrix     0.098

Direction #1
------------
Number of lags              = 10
Direction coefficients      =      1.000     0.000
Direction angles (degrees)  =      0.000     0.000
Tolerance on direction      =     90.000 (degrees)
Calculation lag             =    239.493
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0    79.000    88.901     0.033
         1   862.000   249.508     0.052
         2  1156.000   480.023     0.091
         3  1293.000   718.387     0.113
         4  1255.000   955.540     0.126
         5  1083.000  1194.132     0.122
         6   993.000  1433.762     0.113
         7   930.000  1673.616     0.104
         8   822.000  1908.365     0.096
         9   711.000  2154.216     0.098

# Display the omni-directional experimental variogram with 10 lags using  vario.plot() function
vario.plot(Vario_E1,npairdw = 2, npairpt=1, varline =TRUE, title="Omni-directional Experimental Variogram - 10 Lags", col="blue", lty=2, xlab="Distance(m)", ylab="Variogram")

• Change the number of lags and observe what happens to the variogram

# change nlag to 20
nlag=20

# Compute the omni-directional experimental variogram with 20 lags. you can call it Vario_E2
Vario_E2=vario.calc(db.data,nlag=nlag,lag=D/(2*nlag))

# display the variogram characteristics
print(Vario_E2)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

Variogram characteristics
=========================
Number of variable(s)       = 1
Number of direction(s)      = 1
Space dimension             = 2
Variance-Covariance Matrix     0.098

Direction #1
------------
Number of lags              = 20
Direction coefficients      =      1.000     0.000
Direction angles (degrees)  =      0.000     0.000
Tolerance on direction      =     90.000 (degrees)
Calculation lag             =    119.747
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0     6.000    52.302     0.015
         1   244.000   132.785     0.032
         2   445.000   242.974     0.055
         3   516.000   360.886     0.075
         4   604.000   479.915     0.089
         5   598.000   599.009     0.102
         6   665.000   719.596     0.112
         7   635.000   838.913     0.125
         8   631.000   956.650     0.121
         9   578.000  1075.818     0.131
        10   556.000  1195.374     0.123
        11   505.000  1316.436     0.121
        12   512.000  1437.127     0.107
        13   460.000  1555.978     0.107
        14   471.000  1674.099     0.106
        15   467.000  1796.698     0.099
        16   415.000  1917.253     0.096
        17   356.000  2038.413     0.093
        18   334.000  2154.617     0.102
        19   322.000  2271.511     0.106

# Plot the omni-directional experimental variogram with 20 lags
vario.plot(Vario_E2,npairdw = 2,npairpt=1,varline =TRUE,title="Omni-directional Experimental Variogram - 20 Lags",col="blue",lty=2,xlab="Distance(m)", ylab="Variogram")

Questions: Which variogram should you choose? 10 lags or 20 lags?

There are no set rules for determining what lag size and number of lags should be used. It is the craft of the researcher, their knowledge of the phenomenon they are analyzing, and the reason(s) for modeling a variogram that help to determine the appropriate lag size. Consider what happens when we take lag size to its extremes:

• A lag size = 0 will produce a variogram cloud that perfectly displays all possible pairings, but makes interpretation of the variogram structure difficult.

• At a lag size = infinity (or a distance at least as large as the maximum distance between any two samples), we get one value represented by a single point that represents the average distance and average variogram value for all sample pairings.
  

• Selecting an appropriate lag size between these extremes allows for the creation of a manageable semivariogram to aid in interpretation.

In most cases, the goal is to have as many pairs of points as possible represented in any one variogram point. More pairs per variogram point, however, means a wider bin, and a wider lag distance typically results in less structure for the first points in a variogram.

There are two rules of thumb for selecting a lag size/lag number:

1. Have at least 30-50 pairs minimum for any one variogram point. Smaller bins or lag size means less pairs and probably better structure, but too small a bin or lag size typically introduces more noise into the variogram.

2. The lag size times the number of lags should be about half of the largest distance among all points (half diagonal distance). This condition taken care of by RGeostats for us. Remeber why we defined D.

3. A variogram with as small as possible lag size is preferred so that the spatial variability in the short-range will not be smoothed.

4. At least 3 variogram plots are needed before reaching sill (i.e. before losing spatial correlation)

• Fit automatically a variogram model on the omni-directional experimental variogram

The experimental (or observed) variograms, which represent your sample data, are then fit to each of the 3 types of variogram models within the course. This is so that the data, which were sampled at discreet units, can be modeled as a continuous function, and the value for any unknown point at any distance can be interpolated/estimated. Once the best fit has been made, you can proceed with the interpolation process by kriging methods (simple or ordinary kriging).

# Fit automatically an isotropic model on Vario_E1 using model.auto() and store the model as Vario_M1   
Vario_M1 = model.auto(Vario_E1, struct = melem.name(c(1,3)),draw=FALSE)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

# Display model characteristics
print(Vario_M1)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

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.011
Spherical
- Sill         =      0.105
- Range        =    873.347
Total Sill     =      0.116

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

# plot  experimental variogram Vario_E1 (10 lags) 
vario.plot(Vario_E1, npairdw=1, npairpt=1, varline=FALSE, title="Isotropic Model (nugget effect + spherical)", reset=TRUE, col="blue", lty=2, xlab="Distance(m)",ylab="Variogram")

# display the model on top of the variogram plot
model.plot(Vario_M1 ,vario=Vario_E1,add=T,col="red")

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

# add legend to the plot
legend("bottomright",legend=c("Empirical Variogram","Fitted Variogram (Model)"),lty=c(2,1),cex=1, col=c("blue","red"))

Note:

• The nugget is the variance between very close observations.
• The range is 873 meters in the above variogram. It indicates the distance after which observations are independent . At this distance, there is maximal variance which corresponds to the sill.
• Here, a combined model (total sill=0.116) with a nugget effect model (sill=0.011) and an spherical model (sill=0.105) is fitted to the experimental variogram. The model can be ragarded as a theoretical variogram which provides a smooth variation of spatial variance as a function of distance between spatial locations.

# Fit an Exponential isotropic model on Vario_E1 and display model characteristics 
Vario_M2 = model.auto(Vario_E1, struct = melem.name(c(1,2)), draw=FALSE)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

print(Vario_M2)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

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

Covariance Part
---------------
Exponential
- Sill         =      0.126
- Range        =   1258.019
- Theo. Range  =    419.937
Total Sill     =      0.126

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

# Fit Gaussian isotropic model on Vario_E1 and display model characteristics
Vario_M3 = model.auto(Vario_E1,struct = melem.name(c(1,4)),draw=FALSE)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

print(Vario_M3)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

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.028
Gaussian
- Sill         =      0.088
- Range        =    753.473
- Theo. Range  =    435.328
Total Sill     =      0.116

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

Note: The argument struct = melem.name() controls which model to fit.

• struct = melem.name(c(1,2) fits combination of nugget effect and exponential models.

• struct = melem.name(c(1,3) fits combination of nugget effect and spherical models.
  

• struct = melem.name(c(1,4) fits combination of nugget effect and gaussian models.

# plot experimental variogram with 10 lags (Vario_E1)
vario.plot(Vario_E1, npairdw=1, npairpt=0, varline=FALSE, title="Isotropic Models", reset=TRUE, col="blue", lty=2, xlab="Distance(m)",ylab="Variogram")

# display all the models on top of the variogram plot
model.plot(Vario_M1 ,vario=Vario_E1,add=T,col="red")

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

model.plot(Vario_M2 ,vario=Vario_E1,add=T,col="green")

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

model.plot(Vario_M3 ,vario=Vario_E1,add=T,col="black")

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

# add legend to the plot
legend("bottomright",legend=c("Experimental variogram", "Spherical model","Exponential model","Gaussian model"), lty=c(2,1,1,1),cex=0.8, col=c("blue", "red","green","black"))

Variogram Models:

• Spherical Models: The most commonly used model, with a somewhat linear behavior at small separation distances near the origin, but flattening out at larger distances and reaching a sill limit.

• Exponential Models: Reach the sill asymptotically. Like the spherical model, the exponential model is linear at small distances near the origin, yet rises more steeply and flattens out more gradually. Erratic data sets can sometimes be fit better with exponential models.

• Gaussian Models: Characterized by parabolic behavior at the origin, then rising to reach its sill asymptotically. This model is used to model extremely continuous phenomena.

• Nugget Effect: In models with the nugget effect, the variance does not go through the origin of the plot, indicating that even at very close distances (indeed, even at a distance of zero) the data points show some degree of variability. This can happen when change occurs over the surface at distances less than the sampling interval.

• Task 4: Directional Variography Analysis

• Compute the experimental variogram in four directions:
- north-south Direction (90°)
- east-west Direction (0°)
- north-west South-East Direction (45°)
- north-east South-West Direction (135°)

# number of directions 
ndir=4 

# first direction
ang0=0 

# directions vector
dir_vect=((seq(1,ndir)-1) * 180 / ndir+ ang0) # or simply dir_vec = c(0,45,90,135)
print(dir_vect)

[1]   0  45  90 135

# number of lags
nlag=10

# compute directional variogram
Vario_dir=vario.calc(db.data, dirvect=dir_vect, lag=D/(2*nlag), nlag=nlag)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.
Warning: Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.
Warning: Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.
Warning: Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

# display  variogram characteristics
print(Vario_dir)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

Variogram characteristics
=========================
Number of variable(s)       = 1
Number of direction(s)      = 4
Space dimension             = 2
Variance-Covariance Matrix     0.098

Direction #1
------------
Number of lags              = 10
Direction coefficients      =      1.000     0.000
Direction angles (degrees)  =      0.000     0.000
Tolerance on direction      =     22.500 (degrees)
Calculation lag             =    239.493
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0    22.000    87.912     0.020
         1   200.000   248.620     0.059
         2   231.000   477.365     0.109
         3   249.000   713.427     0.135
         4   180.000   948.658     0.188
         5   131.000  1188.232     0.189
         6    62.000  1417.947     0.153
         7    27.000  1648.575     0.098
         8     6.000  1899.022     0.036

Direction #2
------------
Number of lags              = 10
Direction coefficients      =      0.707     0.707
Direction angles (degrees)  =     45.000     0.000
Tolerance on direction      =     22.500 (degrees)
Calculation lag             =    239.493
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0    14.000    89.483     0.022
         1   254.000   246.610     0.036
         2   363.000   482.863     0.053
         3   486.000   722.335     0.076
         4   602.000   958.703     0.086
         5   622.000  1198.084     0.087
         6   660.000  1438.580     0.091
         7   676.000  1674.412     0.090
         8   638.000  1907.329     0.082
         9   561.000  2153.920     0.079

Direction #3
------------
Number of lags              = 10
Direction coefficients      =      0.000     1.000
Direction angles (degrees)  =     90.000     0.000
Tolerance on direction      =     22.500 (degrees)
Calculation lag             =    239.493
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0    18.000    94.231     0.051
         1   229.000   254.770     0.047
         2   341.000   481.181     0.087
         3   354.000   720.277     0.117
         4   348.000   956.839     0.142
         5   288.000  1188.486     0.157
         6   249.000  1426.902     0.166
         7   210.000  1672.172     0.156
         8   170.000  1912.897     0.153
         9   147.000  2155.917     0.172

Direction #4
------------
Number of lags              = 10
Direction coefficients      =     -0.707     0.707
Direction angles (degrees)  =    135.000     0.000
Tolerance on direction      =     22.500 (degrees)
Calculation lag             =    239.493
Tolerance on distance       =     50.000 (Percent of the lag value)

For variable 1
      Rank    Npairs  Distance     Value
         0    25.000    85.607     0.038
         1   179.000   247.882     0.073
         2   221.000   476.352     0.138
         3   204.000   711.754     0.166
         4   125.000   946.607     0.183
         5    42.000  1192.718     0.203
         6    22.000  1411.434     0.090
         7    17.000  1699.572     0.041
         8     8.000  1901.688     0.051
         9     3.000  2126.160     0.093

# plot the directional variogram
vario.plot(Vario_dir, npairdw = 1, npairpt=0, varline=FALSE, title="Directional Experimental Variograms", lty=2, xlab="Distance(m)", ylab="Variogram")

# add a legend
legend("bottomright",legend=c("Direction 0°","Direction 45°","Direction 90°","Direction 135°"),lty=c(2,2,2,2),col=c(1,2,3,4), cex=0.8)

• Fit automatically anisotropic variogram model on directional experimental variograms

# fit automatically  anisotropic variogram model 
Vario_M4=model.auto(Vario_dir, struct = melem.name(c(1,2)),auth.aniso=TRUE, auth.rotation=TRUE, auth.locksame=TRUE,draw=FALSE)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

print(Vario_M4)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'

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.011
Exponential
- Sill         =      0.202
- Ranges       =   1833.689 10357.138
- Theo. Ranges =    612.100  3457.298
- Angles       =    -35.689     0.000
- Rotation Matrix
               [,  0]    [,  1]
     [  0,]     0.812     0.583
     [  1,]    -0.583     0.812
Total Sill     =      0.213

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

# plot the directional variogram
vario.plot(Vario_E4,npairdw=1,npairpt=0,varline =FALSE,title="Anisotropic Model",reset=TRUE,lty=2, xlab="Distance(m)",ylab="Variogram")

# plot the anisotropic variogram model
model.plot(Vario_M4 ,vario=Vario_E4,add=T)

Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: is.na() applied to non-(list or vector) of type 'S4'Warning: Recycling array of length 1 in vector-array arithmetic is deprecated.
  Use c() or as.vector() instead.

# add a legend
legend("bottomright", legend=c("Direction 0°","Direction 45°","Direction 90°","Direction 135°"), lty=c(1,1,1,1), col=c(1,2,3,4), cex=0.8)

• Questions:
- What do you observe? - What is the direction of maximum spatial continuity? - can you detect Anisotropy ?

Isotropy and Anisotropy: general speaking, the procedures of distinguishing between isotropy and anisotropy are the same. We can fit a variogram model along all the directions and then obtain different variograms along with different directions. Once we obtain the values of Still, Range, Nugget Effect from all variograms, we can determine if the spatial data is isotropic or anisotropic. If the values of Still, Range, Nugget Effect from the variogram along with all directions are all the same, then it is isotropic; otherwise, it is the anisotropic.

Anisotropy, according to its different structures, could be divided into:

• Geometric Anisotropy: the sill value of the geometric anisotropic variogram only vary along with distances, but not along with directions. In other words, when the distances between samples in various spatial angles are the same, their Sill is the same, but their ranges may vary along with different angles. The direction having biggest range is the direction of variogram. This range is called “Major Range”, and this variogram is called Directional Variogram. The smallest range is called “Minor Range”.

Fig 1: Geometric anisotropy

• Zonal Anisotropy: the Sill value of the zonal anisotropic variogram does not only vary along with all distances, but also varies along with all directions. In other words, when the distances between samples in any spatial angles are the same, their Sill values are not the same, and their ranges also vary along with different angles. Zonal Anisotropic variogram model consists of 2 or more anisotropic variograms. Natural phenomena usually do not have this kind of anisotropy.

Fig 2: Zonal anisotropy

• Mixed Anisotropy: is a mix of the Geometric Anisotropy and Zonal Anisotropy . Its Sill values and ranges vary along with all directions. Natural phenomena often have this kind of spatial variability. For example, geological changes usually have a bigger effective range in the horizon direction than in the vertical direction, and changes in the vertical driection have a bigger range than in the horizon direction.

Fig 3: Mixed anisotropy