Outline

Semivariance
Variogram
Ordinary Kriging
Cokriging

Semivariance

Regionalized variable theory uses a related property called the semivariance to express the degree of relationship between points on a surface.
The semivariance is simply half the variance of the differences between all possible points spaced a constant distance apart.

Semivariance is a measure of the degree of spatial dependence between samples (elevation)

##Meuse River Dataset

library(sp)
data(meuse)
class(meuse)

## [1] "data.frame"

coordinates(meuse) <- c("x", "y")
class(meuse)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

str(meuse)

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  155 obs. of  12 variables:
##   .. ..$ cadmium: num [1:155] 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
##   .. ..$ copper : num [1:155] 85 81 68 81 48 61 31 29 37 24 ...
##   .. ..$ lead   : num [1:155] 299 277 199 116 117 137 132 150 133 80 ...
##   .. ..$ zinc   : num [1:155] 1022 1141 640 257 269 ...
##   .. ..$ elev   : num [1:155] 7.91 6.98 7.8 7.66 7.48 ...
##   .. ..$ dist   : num [1:155] 0.00136 0.01222 0.10303 0.19009 0.27709 ...
##   .. ..$ om     : num [1:155] 13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
##   .. ..$ ffreq  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ soil   : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
##   .. ..$ lime   : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
##   .. ..$ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
##   .. ..$ dist.m : num [1:155] 50 30 150 270 380 470 240 120 240 420 ...
##   ..@ coords.nrs : int [1:2] 1 2
##   ..@ coords     : num [1:155, 1:2] 181072 181025 181165 181298 181307 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:155] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] 178605 329714 181390 333611
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

plot(meuse, asp = 1, pch = 1)
data(meuse.riv)
lines(meuse.riv)

Display a postplot of the untransformed Zn values, that is, plot the sample locations (as above) and represent the data value by the size of the symbol.

plot(meuse, asp = 1,
cex = 4*meuse$zinc/max(meuse$zinc),
pch = 1)
lines(meuse.riv)

Compute the number of point-pairs

meuse$logZn<-log10(meuse$zinc)
n <- length(meuse$logZn)
n * (n - 1)/2

## [1] 11935

Compute the distance and semivariance between the first two points in the data set

dim(coordinates(meuse))

## [1] 155   2

coordinates(meuse)[1, ]

##      x      y 
## 181072 333611

coordinates(meuse)[2, ]

##      x      y 
## 181025 333558

(sep <- dist(coordinates(meuse)[1:2, ]))

##          1
## 2 70.83784

(gamma <- 0.5 * (meuse$logZn[1] - meuse$logZn[2])^2)

## [1] 0.001144082

Plot the experimental variogram of the log-Zn concentrations

library(gstat)
v <- variogram(logZn ~ 1, meuse, cutoff = 1300, width = 90)

Plot the experimental variogram of the log-Zn concentrations

print(plot(v, plot.numbers = T))

The trend in decreasing semivariance with decreasing separation seems to intersect the yaxis (i.e., at 0 separation) at about 0.01 log(mg kg-1)2; this is the nugget.

Variogram

The semivariance at a distance d = 0 should be zero, because there are no differences between points that are compared to themselves.
However, as points are compared to increasingly distant points, the semivariance increases.

Display the variogram model forms which can be used in gstat

print(show.vgms())

Fit a spherical variogram model

vm <- vgm(psill = 0.12,
          model = "Sph",
          range = 850,
          nugget = 0.01)
print(plot(v, pl = T,
      model = vm))

Adjust the model with gstat automatic fit

vm

##   model psill range
## 1   Nug  0.01     0
## 2   Sph  0.12   850

vmf = fit.variogram(v,vm)
print(plot(v, pl = T, model = vmf))

Ordinary Kriging

The theory of regionalised variables leads to an ”optimal“ prediction method, in the sense that the kriging variance is minimized.

What is so special about kriging?

Predicts at any point as the weighted average of the values at sampled points
Weights given to each sample point are optimal, given the spatial covariance structure as revealed by the variogram model (in this sense it is “best”)
The kriging variance at each point is automatically generated as part of the process of computing the weights

Ordinary kriging on a regular grid

Load the 40 m x 40 m interpolation grid covering the sample area and convert it to a spatial object

data(meuse.grid)
coordinates(meuse.grid) <- c("x", "y")
gridded(meuse.grid) <- T

Predict the attribute value at all grid points using Ordinary Kriging

k40 <- krige(logZn ~ 1, locations = meuse, newdata = meuse.grid, model = vmf)

## [using ordinary kriging]

#[using ordinary kriging]

Display the structure of the kriging prediction object

str(k40)

## Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
##   ..@ data       :'data.frame':  3103 obs. of  2 variables:
##   .. ..$ var1.pred: num [1:3103] 2.83 2.88 2.83 2.78 2.94 ...
##   .. ..$ var1.var : num [1:3103] 0.0596 0.0471 0.0509 0.055 0.0335 ...
##   ..@ coords.nrs : int [1:2] 1 2
##   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
##   .. .. ..@ cellcentre.offset: Named num [1:2] 178460 329620
##   .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
##   .. .. ..@ cellsize         : Named num [1:2] 40 40
##   .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
##   .. .. ..@ cells.dim        : Named int [1:2] 78 104
##   .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
##   ..@ grid.index : int [1:3103] 69 146 147 148 223 224 225 226 300 301 ...
##   ..@ coords     : num [1:3103, 1:2] 181180 181140 181180 181220 181100 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:3103] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] 178440 329600 181560 333760
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

Display the map of predicted values

print(spplot(k40, "var1.pred", 
             asp=1, 
             main="OK prediction, log-ppm Zn"))

Display the map of kriging prediction variances

print(spplot(k40, "var1.var", 
             col.regions=cm.colors(64), 
              asp=1, 
              main="OK prediction variance,
              log-ppm Zn^2"))

Describe the variances map:

Where is the prediction variance lowest?
Does this depend on the data value?

Show the post-plot: value proportional to circle size

pts.s <- list("sp.points", meuse, col="white",
pch=1, cex=4*meuse$zinc/max(meuse$zinc))
print(spplot(k40, "var1.pred", asp=1,
col.regions=bpy.colors(64),
main="OK prediction, log-ppm Zn",
sp.layout = list(pts.s)))

Show the observation locations on the krigging prediction variance

pts.s <- list("sp.points", meuse, col="black",
pch=20)
print(spplot(k40, zcol="var1.var",
col.regions=cm.colors(64),
asp=1,main="OK prediction variance,
log-ppm Zn^2",sp.layout = list(pts.s)))

Evaluating the model

zn.okcv<-krige.cv(log(zinc)~1, meuse, meuse.grid, model = vmf)
RMSE.ok<-sqrt(sum(zn.okcv$residual^2)/length(zn.okcv$residual))
RMSE.ok

## [1] 0.3973315

Co-Kriging

Co-kriging allows samples of an auxiliary variable (also called the co-variable), besides the target value of interest, to be used when predicting the target value at unsampled locations. The co-variable may be measured at the same points as the target (co-located samples), at other points, or both.
The most common application of co-kriging is when the co-variable is cheaper to measure, and so has been more densely sampled, than the target variable.

Eksplorasi Data

data(meuse.riv)# outline of the river
meuse.lst <- list(Polygons(list(Polygon(meuse.riv)), "meuse.riv"))
meuse.sr <- SpatialPolygons(meuse.lst)
image(meuse.grid["dist"])# one of the variables in meuse.grid
plot(meuse.sr, add=TRUE)
title("distance to river")

zn.lm <- lm(log10(zinc) ~ sqrt(dist), meuse)
meuse$fitted.s <- predict(zn.lm,meuse)
-mean(predict(zn.lm,meuse))

## [1] -2.55616

meuse$residuals <- residuals(zn.lm)
spplot(meuse,c("fitted.s","residuals"))

Generate a new kriged surface with the covariate - distance to river

vm.ck <- variogram(log10(zinc)~sqrt(dist), meuse)
plot(vm.ck, plot.numbers = TRUE, pch = "+")

Secara visual, berapakah nilai sill, range, dan nugget berdasarkan variogram di samping?

m.ck <- fit.variogram(vm.ck,vgm(.03,
"Sph", 800,.015))
plot(vm.ck, model=m.ck)

ko.ck <- krige(log(zinc)~ sqrt(dist), meuse,
meuse.grid,model=m.ck)

## [using universal kriging]

pts <- list("sp.points",meuse,pch=3,col="black")
meuse.layout <- list(pts)
spplot(ko.ck["var1.pred"], sp.layout=meuse.layout, main="cokriging predictions-Zn/distance river ")

ko.ck$sek <- sqrt(ko.ck$var1.var)
spplot(ko.ck,zcol='sek',
sp.layout=meuse.layout, main =
"co-kriging se-Zn(covariate)")

summary(ko.ck)

## Object of class SpatialPixelsDataFrame
## Coordinates:
##      min    max
## x 178440 181560
## y 329600 333760
## Is projected: NA 
## proj4string : [NA]
## Number of points: 3103
## Grid attributes:
##   cellcentre.offset cellsize cells.dim
## x            178460       40        78
## y            329620       40       104
## Data attributes:
##    var1.pred        var1.var            sek        
##  Min.   :4.455   Min.   :0.01903   Min.   :0.1379  
##  1st Qu.:5.212   1st Qu.:0.02187   1st Qu.:0.1479  
##  Median :5.590   Median :0.02322   Median :0.1524  
##  Mean   :5.702   Mean   :0.02450   Mean   :0.1561  
##  3rd Qu.:6.134   3rd Qu.:0.02595   3rd Qu.:0.1611  
##  Max.   :7.477   Max.   :0.03983   Max.   :0.1996

Evaluasi

zn.ckcv<-krige.cv(log(zinc)~ sqrt(dist), meuse,meuse.grid,model=m.ck)
RMSE.ck<-sqrt(sum(zn.ckcv$residual^2)/length(zn.ckcv$residual))
RMSE.ck

## [1] 0.3752662

Back Transform

ko.ck$predt <- 10^(ko.ck$var1.pred)
spplot(ko.ck["predt"], sp.layout=meuse.layout,
main = "Co-Kriging predictions-Meuse
zinc log/backtrans(Zn)")

Target Variable

We select lead (abbreviation “Pb”) as the target variable, i.e. the one we want to map. This metal is a serious human health hazard.
It can be inhaled as dust from disturbed soil or taken up by plants and ingested.
The critical value for Pb in agricultural soils, according to the Berlin Digital Environmental Atlas2, is 600 mg kg-1 for agricultural fields: above this level grain crops can not be grown for human consumption.

Selecting the co-variables

Candidates for co-variables must have:

a feature-space correlation with the target variable;
a spatial structure (i.e. be modelled as a regional variable);
a spatial co-variance with the target variable.

Two main ways to select a co-variable:

theoretically, from knowledge of the spatial process that caused the observed spatial (co-)distribution;
empirically, by examining the feature-space correlations (scatterplots) and then the spatial co-variance (cross-correlograms or crossvariograms).

Kandidat Peubah Penjelas

organic matter content (OM)
zinc content (Zn)

Exercise 1

Lakukan prediksi Pb dengan menyertakan co-variables (OM atau Zn), diskusikan hasil yang Anda peroleh.

Exercise 2

Perhatikan data yg tersedia pada alamat berikut: (https://raw.githubusercontent.com/raoy/SpatialStatistics/master/database_Nickel.csv)
Keterangan variabel:
XCOLLAR: Longitude
YCOLLAR: Latitude
ZCOLLAR: Kedalaman
Ni: Kandungan Nikel (target variable)

Coba lakukan interpolasi dengan beberapa metode yang telah Anda pelajari, bandingkan hasilnya, manakah interpolasi yang paling baik menurut Anda?

Kriging dan CoKriging

rahmaanisa@apps.ipb.ac.id

Praktikum 5 | Statistika Spasial

Outline

Semivariance

Compute the number of point-pairs

Compute the distance and semivariance between the first two points in the data set

Plot the experimental variogram of the log-Zn concentrations

Plot the experimental variogram of the log-Zn concentrations

Variogram

Display the variogram model forms which can be used in gstat

Fit a spherical variogram model

Adjust the model with gstat automatic fit

Ordinary Kriging

What is so special about kriging?

Ordinary kriging on a regular grid

Display the structure of the kriging prediction object

Display the map of predicted values

Display the map of kriging prediction variances

Show the post-plot: value proportional to circle size

Show the observation locations on the krigging prediction variance

Evaluating the model

Co-Kriging

Eksplorasi Data

Generate a new kriged surface with the covariate - distance to river

Evaluasi

Back Transform

Target Variable

Selecting the co-variables

Kandidat Peubah Penjelas

Exercise 1

Exercise 2