Introduction

Variogram analysis is a technique used in geostatistics to analyze the spatial dependence or correlation of a variable across a geographic area.
In this presentation, we will walk through an R script that performs variogram analysis.

Semivariogram

Semivariogram is a function describing the degree of spatial correlation of a spatial random variable.

In spatial modeling, semivariogram begins with a graph of the empirical semivariogram, which is the half of average squared difference between points separated by a distance.

The semivariogram is calculated as:

>- Where, Nh is the number of pairs separated by vector h, >- vector h is lag distance, >- xi is the starting location and xi+h is the ending location. >- If y is only dependent on the length of lag distance but not its direction.

Data Import and Experimental Variogram

## The legacy packages maptools, rgdal, and rgeos, underpinning this package
## will retire shortly. Please refer to R-spatial evolution reports on
## https://r-spatial.org/r/2023/05/15/evolution4.html for details.
## This package is now running under evolution status 0

## Warning: package 'openxlsx' was built under R version 4.3.1

## Warning: package 'writexl' was built under R version 4.3.1

## Warning: package 'broom' was built under R version 4.3.1

## Loading required package: mvtnorm

## Loading required package: survival

## Loading required package: TH.data

## Loading required package: MASS

## 
## Attaching package: 'TH.data'

## The following object is masked from 'package:MASS':
## 
##     geyser

## Warning: package 'gridExtra' was built under R version 4.3.1

VARIOGRAM ANALYSIS

Data needs to be comprised of 3 columns Column 1 = latitude (x coordinates), Column 2 = longitude (y coordinates), Column 3 = z (desired measurement per point in the space).

#DEVELOPING EXPERIMENTAL VARIOGRAM

Load the required packages

##FILL THESE input_path <- “C:\Users\barbi\Desktop\variogram-example.xlsx” output_path <- “C:\Users\barbi\Desktop\variogram_output.xlsx”

##IMPORTING DATA

data <- read_excel(input_path)

data(meuse)

Convert data to a spatial points data frame (SPDF)

coordinates(meuse) <- ~x + y

##GENERATING EXPERIMENTAL VARIOGRAM

Create experimental variogram

exp_variogram <- variogram(zinc ~ 1, meuse)

Results of experimental variogram

head(exp_variogram)

##    np      dist     gamma dir.hor dir.ver   id
## 1  57  79.29244  37362.96       0       0 var1
## 2 299 163.97367  72718.34       0       0 var1
## 3 419 267.36483  82655.53       0       0 var1
## 4 457 372.73542 111575.91       0       0 var1
## 5 547 478.47670 123080.69       0       0 var1
## 6 533 585.34058 147414.28       0       0 var1

plot(exp_variogram, main = "Variogram - default", xlab = "Separation distance (m)")

Some variogram functions in gstat

show.vgms()

Extract distance and gamma values

exp_variogram_data <- data.frame(distance = exp_variogram$dist, gamma = exp_variogram$gamma)

Plot the experimental variogram using ggplot 2

p=ggplot(exp_variogram_data, aes(x = distance, y = gamma)) +
  geom_point() +
  geom_line() +
  labs(x = "Lag (h)", y = "Semivariance") +
  theme_minimal()

Experimental Variogram

Extraction of parameters

sill <- max(exp_variogram_data$gamma)
nugget <- exp_variogram_data$gamma[1]
range <- exp_variogram_data$distance[which.max(exp_variogram_data$gamma)]

cat("Sill: ", sill, "\n")

## Sill:  186528.2

cat("Nugget: ", nugget, "\n")

## Nugget:  37362.96

cat("Range: ", range, "\n")

## Range:  1117.862

##GENERATING VARIOGRAM MODELS

Generate variogram models

spherical_model <- vgm(psill = sill, model = "Sph", range = range, nugget = nugget)

exponential_model <- vgm(psill = sill, model = "Exp", range = range, nugget = nugget)

linear_model <- vgm(psill = sill, model = "Lin", range = range, nugget = nugget)

gaussian_model <- vgm(psill = sill, model = "Gau", range = range, nugget = nugget)

Compute variogram values for the models (spherical_variogram)

spherical_variogram <- variogramLine(spherical_model, maxdist = max(exp_variogram_data$distance))

Plot variogram values for the models (spherical_variogram)

plot(spherical_variogram)

Plot the variograms (Spherical)

initial_fit <- ggplot(exp_variogram_data, aes(x = distance, y = gamma)) +
  geom_point(color = "blue") +
  geom_line(data = spherical_variogram, aes(x = dist, y = gamma), color = "red") +
  labs(x = "Lag (h)", y = "Semivariance", title = " Spherical Variogram Model")+
  theme_minimal()

Experimental vs Theoretical (Spherical) Varigram model

plot(initial_fit)

Compute variogram values for the models (Exponential_variogram)

exponential_variogram <- variogramLine(exponential_model, maxdist = max(exp_variogram_data$distance))

Plot variogram values for the models (Exponential_variogram)

plot(exponential_variogram)

Plot the variograms (Eponential)

initial_fit <- ggplot(exp_variogram_data, aes(x = distance, y = gamma)) +
  geom_point(color = "blue") +
  geom_line(data = exponential_variogram, aes(x = dist, y = gamma), color = "green") +
  labs(x = "Lag (h)", y = "Semivariance", title = " Eponential Variogram Model") +
  theme_minimal()

Experimental vs Theoretical (Exponential) Varigram model

plot(initial_fit)

Compute variogram values for the models (linear_variogram)

linear_variogram <- variogramLine(linear_model, maxdist = max(exp_variogram_data$distance))

Compute variogram values for the models (Linear_variogram)

plot(linear_variogram)

Plot the variograms (Linear Varioram)

Experimental vs Theoretical (Linear) Varigram model

plot(initial_fit)

Compute variogram values for the models (Gaussian_variogram)

gaussian_variogram <- variogramLine(gaussian_model, maxdist = max(exp_variogram_data$distance))

Compute variogram values for the models (Gaussian_variogram)

plot(gaussian_variogram)

Plot the variograms (Gaussian)

Experimental vs Theoretical (Gaussian) Varigram model

plot(initial_fit)

Calculate RMSE for each model

spherical_rmse <- sqrt(mean((spherical_variogram$gamma - exp_variogram_data$gamma)^2))

## Warning in spherical_variogram$gamma - exp_variogram_data$gamma: longer object
## length is not a multiple of shorter object length

exponential_rmse <- sqrt(mean((exponential_variogram$gamma - exp_variogram_data$gamma)^2))

## Warning in exponential_variogram$gamma - exp_variogram_data$gamma: longer
## object length is not a multiple of shorter object length

linear_rmse <- sqrt(mean((linear_variogram$gamma - exp_variogram_data$gamma)^2))

## Warning in linear_variogram$gamma - exp_variogram_data$gamma: longer object
## length is not a multiple of shorter object length

gaussian_rmse <- sqrt(mean((gaussian_variogram$gamma - exp_variogram_data$gamma)^2))

## Warning in gaussian_variogram$gamma - exp_variogram_data$gamma: longer object
## length is not a multiple of shorter object length

RMSE for each model

rmse=cbind(spherical_rmse,exponential_rmse,linear_rmse,gaussian_rmse)
rmse

##      spherical_rmse exponential_rmse linear_rmse gaussian_rmse
## [1,]       81064.01          57278.9    77362.37      70296.15

Create a data frame for the legend labels

legend_df <- data.frame(model = c("Spherical", "Exponential", "Linear", "Gaussian"),
                        color = c("red", "green", "purple", "orange"))

Plot all the four the variograms

initial_fit <- ggplot(exp_variogram_data, aes(x = distance, y = gamma)) +
  geom_point(color = "blue") +
  geom_line(data = spherical_variogram, aes(x = dist, y = gamma), color = "red") +
  geom_line(data = exponential_variogram, aes(x = dist, y = gamma), color = "green") +
  geom_line(data = linear_variogram, aes(x = dist, y = gamma), color = "purple") +
  geom_line(data = gaussian_variogram, aes(x = dist, y = gamma), color = "orange") +
  labs(x = "Lag (h)", y = "Semivariance", title = "Variogram Models") +
  scale_color_manual(values = legend_df$color, labels = legend_df$model)+
  theme_minimal()

Varioram Models

initial_fit

Variogram Modeling

Zulfiqar Ali, Ph.D (Statistics, PAK/NL)

November 5, 2023

Introduction

Semivariogram

Data Import and Experimental Variogram

VARIOGRAM ANALYSIS

Load the required packages

Convert data to a spatial points data frame (SPDF)

Create experimental variogram

Results of experimental variogram

Some variogram functions in gstat

Extract distance and gamma values

Plot the experimental variogram using ggplot 2

Experimental Variogram

Extraction of parameters

Generate variogram models

Compute variogram values for the models (spherical_variogram)

Plot variogram values for the models (spherical_variogram)

Plot the variograms (Spherical)

Experimental vs Theoretical (Spherical) Varigram model

Compute variogram values for the models (Exponential_variogram)

Plot variogram values for the models (Exponential_variogram)

Plot the variograms (Eponential)

Experimental vs Theoretical (Exponential) Varigram model

Compute variogram values for the models (linear_variogram)

Compute variogram values for the models (Linear_variogram)

Plot the variograms (Linear Varioram)

Experimental vs Theoretical (Linear) Varigram model

Compute variogram values for the models (Gaussian_variogram)

Compute variogram values for the models (Gaussian_variogram)

Plot the variograms (Gaussian)

Experimental vs Theoretical (Gaussian) Varigram model

Calculate RMSE for each model

RMSE for each model

Create a data frame for the legend labels

Plot all the four the variograms

Varioram Models