[#R] Short tutorial variogram using geoR package

A tutorial how to make variogram: an exercise

Author: Dasapta Erwin Irawan

Just starting to work on the variogram analysis of the Cikapundung dataset using geoR package. This post is connected to these previous posts:

The revised code can be downloaded here: http://goo.gl/7I4cuV

1. #Pairs function short tutorial
2. [#R] How to convert lat-long coordinates to UTM (easting-northing)
3. Data preparation and plotting on #QGIS
4. #R: Timeseries analysis try out
5. #R: Hydromad on Cikapundung-update on monthly analysis

The following is the outline of this code is:

1. loading geoR package and data
2. assigning data frame as geo data and checking for duplicate coordinate
3. making variogram with variog function
4. fitting variogram with variofit function

I'll post the explanation later on.

# ANALYSIS FOR 1997 DATA

# LOADING LIBRARY
require("geoR")

## Loading required package: geoR
## Loading required package: sp
## Loading required package: MASS
## --------------------------------------------------------------
##  Analysis of geostatistical data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-4 (built on 2012-06-29) is now loaded
## --------------------------------------------------------------

# require('gstat') require('lattice') require('knitcitations') --------

# LOADING DATA
Data97 <- read.csv("data97utm.csv")

# ANALYSIS FOR 1ST GROUP PARAMETER: ELEVATION, EC, pH, HARD, TDS, TEMP Pairs
# analysis Insert Pairs Code

## Variogram analysis

#### ??? What is variogram? It's basically similar to autocorrelation but in
#### spatial arrangement.

## Binding necessary columns
Data97.1 <- Data97[, c("x", "y", "elv", "ec", "ph", "hard", "tds", "temp")]

### Assigning EC data as geodata (column 8)
EC97 <- as.geodata(Data97.1, coords.col = 1:2, data.col = 4)  #variogram for EC column 

### Checking any duplicate coordinates
dup.coords(EC97)

## NULL

### Building and plotting variogram
Var1.EC97 <- variog(EC97, trend = "1st")

## variog: computing omnidirectional variogram

Var2.EC97 <- variog(EC97, trend = "2nd")

## variog: computing omnidirectional variogram

#### Plot side by side
plot.new()
par(mfrow = c(2, 1))
plot(Var1.EC97, pch = 19, col = "blue", main = "1st order variogram")
plot(Var2.EC97, pch = 19, col = "red", main = "2nd order variogram")

#### Plot overlaying
plot.new()
par(mfrow = c(1, 1))

plot(Var1.EC97, pch = 19, col = "blue", main = "Variogram EC Data 1997")
par(new = TRUE)
plot(Var2.EC97, pch = 19, col = "red", xaxt = "n", yaxt = "n")

##### ??? Question: How to choose the 1st or 2nd order?

### Fitting 1st EC variogram
ini.vals <- expand.grid(seq(15000, 35000, by = 100), seq(1000, 2000, by = 100))
ols <- variofit(Var1.EC97, ini = ini.vals, fix.nug = TRUE, wei = "equal")

## variofit: covariance model used is matern 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "27500" "1400" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 470520606.080271

summary(ols)

## $pmethod
## [1] "OLS (ordinary least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   27509    1392 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 4170
## 
## $sum.of.squares
##     value 
## 470511572 
## 
## $estimated.pars
## sigmasq     phi 
##   27509    1392 
## 
## $weights
## [1] "equal"
## 
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = ini.vals, fix.nugget = TRUE, 
##     weights = "equal")
## 
## attr(,"class")
## [1] "summary.variomodel"

wls <- variofit(Var1.EC97, ini = ini.vals, fix.nug = TRUE, wei = "npairs")

## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "24800" "1100" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 40933024994.4065

summary(wls)

## $pmethod
## [1] "WLS (weighted least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   24777    1114 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 3338
## 
## $sum.of.squares
##     value 
## 4.093e+10 
## 
## $estimated.pars
## sigmasq     phi 
##   24777    1114 
## 
## $weights
## [1] "npairs"
## 
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = ini.vals, fix.nugget = TRUE, 
##     weights = "npairs")
## 
## attr(,"class")
## [1] "summary.variomodel"

lines(wls)
lines(ols, lty = 2, col = "blue")

init.values <- expand.grid(seq(15000, 35000, by = 100), seq(1000, 2000, by = 100))
olsFit11.EC97 <- variofit(Var1.EC97, ini = init.values, fix.nug = T, wei = "equal")

## variofit: covariance model used is matern 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "27500" "1400" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 470520606.080271

olsFit12.EC97 <- variofit(Var1.EC97, ini = init.values, fix.nug = T, wei = "equal")

## variofit: covariance model used is matern 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "27500" "1400" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 470520606.080271

#### ??? Warning on both Fits: unreasonable initial value for sigmasq + nugget
#### (too low), why?
summary(olsFit11.EC97)

## $pmethod
## [1] "OLS (ordinary least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   27509    1392 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 4170
## 
## $sum.of.squares
##     value 
## 470511572 
## 
## $estimated.pars
## sigmasq     phi 
##   27509    1392 
## 
## $weights
## [1] "equal"
## 
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = init.values, fix.nugget = T, 
##     weights = "equal")
## 
## attr(,"class")
## [1] "summary.variomodel"

summary(olsFit12.EC97)

## $pmethod
## [1] "OLS (ordinary least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   27509    1392 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 4170
## 
## $sum.of.squares
##     value 
## 470511572 
## 
## $estimated.pars
## sigmasq     phi 
##   27509    1392 
## 
## $weights
## [1] "equal"
## 
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = init.values, fix.nugget = T, 
##     weights = "equal")
## 
## attr(,"class")
## [1] "summary.variomodel"

wlsFit11.EC97 <- variofit(Var1.EC97, ini = init.values, fix.nug = T, wei = "npairs")

## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "24800" "1100" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 40933024994.4065

wlsFit12.EC97 <- variofit(Var1.EC97, ini = init.values, fix.nug = T, wei = "npairs")

## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "24800" "1100" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 40933024994.4065

#### ??? Warning on both Fits: unreasonable initial value for sigmasq + nugget
#### (too low), why?  ??? What's the different between equal and npairs
summary(wlsFit11.EC97)

## $pmethod
## [1] "WLS (weighted least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   24777    1114 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 3338
## 
## $sum.of.squares
##     value 
## 4.093e+10 
## 
## $estimated.pars
## sigmasq     phi 
##   24777    1114 
## 
## $weights
## [1] "npairs"
## 
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = init.values, fix.nugget = T, 
##     weights = "npairs")
## 
## attr(,"class")
## [1] "summary.variomodel"

summary(wlsFit12.EC97)

## $pmethod
## [1] "WLS (weighted least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   24777    1114 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 3338
## 
## $sum.of.squares
##     value 
## 4.093e+10 
## 
## $estimated.pars
## sigmasq     phi 
##   24777    1114 
## 
## $weights
## [1] "npairs"
## 
## $call
## variofit(vario = Var1.EC97, ini.cov.pars = init.values, fix.nugget = T, 
##     weights = "npairs")
## 
## attr(,"class")
## [1] "summary.variomodel"

### Fitting 2nd EC variogram
init.values <- expand.grid(seq(15000, 35000, by = 100), seq(1000, 2000, by = 100))
olsFit21.EC97 <- variofit(Var2.EC97, ini = init.values, fix.nug = T, wei = "equal")

## variofit: covariance model used is matern 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "28100" "2000" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 515433457.115543

olsFit22.EC97 <- variofit(Var2.EC97, ini = init.values, fix.nug = T, wei = "equal")

## variofit: covariance model used is matern 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "28100" "2000" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 515433457.115543

#### ??? Warning on both Fits: unreasonable initial value for sigmasq + nugget
#### (too low), why?
summary(olsFit21.EC97)

## $pmethod
## [1] "OLS (ordinary least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   28197    2117 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 6342
## 
## $sum.of.squares
##     value 
## 515162980 
## 
## $estimated.pars
## sigmasq     phi 
##   28197    2117 
## 
## $weights
## [1] "equal"
## 
## $call
## variofit(vario = Var2.EC97, ini.cov.pars = init.values, fix.nugget = T, 
##     weights = "equal")
## 
## attr(,"class")
## [1] "summary.variomodel"

summary(olsFit22.EC97)

## $pmethod
## [1] "OLS (ordinary least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   28197    2117 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 6342
## 
## $sum.of.squares
##     value 
## 515162980 
## 
## $estimated.pars
## sigmasq     phi 
##   28197    2117 
## 
## $weights
## [1] "equal"
## 
## $call
## variofit(vario = Var2.EC97, ini.cov.pars = init.values, fix.nugget = T, 
##     weights = "equal")
## 
## attr(,"class")
## [1] "summary.variomodel"

wlsFit21.EC97 <- variofit(Var2.EC97, ini = init.values, fix.nug = T, wei = "npairs")

## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "23000" "1100" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 34351583362.0137

wlsFit22.EC97 <- variofit(Var2.EC97, ini = init.values, fix.nug = T, wei = "npairs")

## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "23000" "1100" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 34351583362.0137

#### ??? Warning on both Fits: unreasonable initial value for sigmasq + nugget
#### (too low), why?  ??? What's the different between equal and npairs
summary(wlsFit21.EC97)

## $pmethod
## [1] "WLS (weighted least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   22975    1127 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 3375
## 
## $sum.of.squares
##     value 
## 3.434e+10 
## 
## $estimated.pars
## sigmasq     phi 
##   22975    1127 
## 
## $weights
## [1] "npairs"
## 
## $call
## variofit(vario = Var2.EC97, ini.cov.pars = init.values, fix.nugget = T, 
##     weights = "npairs")
## 
## attr(,"class")
## [1] "summary.variomodel"

summary(wlsFit22.EC97)

## $pmethod
## [1] "WLS (weighted least squares)"
## 
## $cov.model
## [1] "matern"
## 
## $spatial.component
## sigmasq     phi 
##   22975    1127 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
## tausq 
##     0 
## 
## $fix.nugget
## [1] TRUE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 3375
## 
## $sum.of.squares
##     value 
## 3.434e+10 
## 
## $estimated.pars
## sigmasq     phi 
##   22975    1127 
## 
## $weights
## [1] "npairs"
## 
## $call
## variofit(vario = Var2.EC97, ini.cov.pars = init.values, fix.nugget = T, 
##     weights = "npairs")
## 
## attr(,"class")
## [1] "summary.variomodel"

plot.new()
par(mfrow = c(1, 2))
plot(Var1.EC97, main = "1st order EC Semivar", pch = 19, col = "blue")
lines(wlsFit11.EC97, lty = 2, col = "green")
lines(olsFit11.EC97, lty = 3, col = "blue")
lines(wlsFit12.EC97, lty = 2, col = "grey")
lines(olsFit12.EC97, lty = 3, col = "red")

plot(Var2.EC97, main = "2nd order EC Semivar", pch = 19, col = "red")
lines(wlsFit21.EC97, lty = 2, col = "green")
lines(olsFit21.EC97, lty = 3, col = "blue")
lines(wlsFit22.EC97, lty = 2, col = "grey")
lines(olsFit22.EC97, lty = 3, col = "red")

# legend('topleft',c('ordinary least squares','weighted least
# squares'),lty=c(2,1), lwd=c(1,1), col=c('blue','black'))

####### REFERENCES

# R Core Team (2014). R: A language and environment for statistical
# computing. R Foundation for Statistical Computing, Vienna, Austria. URL
# http://www.R-project.org/.

# RIBEIRO JR., P.J. and DIGGLE, P.J. (2001) geoR: A package for
# geostatistical analysis. R-NEWS Vol 1, No 2. ISSN 1609-3631.