Geostatistics

Hello,

This product is a homework assignment for Western Washington Univeristy’s ESCI597 Time Series and Spatial Analysis course, lead by Dr. Andy Bunn. Most of this work prior to the “My Work” section is attributed to him. Please see his RPubs Page Here.

IDW

library(sp)
library(gstat)

d <- 1:100
w <- d^-1
plot(d,w,type="l",xlab="Distance (d)", ylab="Weight (w)")

Next, we select a value for p

p <- 0
w <- d^-p
plot(d,w,type="l",xlab="Distance (d)", ylab="Weight (w)", ylim = c(0,1),
     col="black")
p <- 0.25
w <- d^-p
lines(d,w,col="red")
p <- 0.5
w <- d^-p
lines(d,w,col="blue")
p <- 1
w <- d^-p
lines(d,w,col="purple")
p <- 2
w <- d^-p
lines(d,w,col="green")
legend("topright",legend=c("p=0","p=0.25","p=0.5","p=1","p=2"),
       col=c("black","red","blue","purple","green"),lwd=2)

Meuse River

Let’s perform an IDW on the Meuse river data

data(meuse)
coordinates(meuse) <- ~x+y
proj4string(meuse) <- CRS("+init=epsg:28992")
spplot(meuse,"cadmium",colorkey=TRUE)

data(meuse.grid)
coordinates(meuse.grid) <- ~x+y
proj4string(meuse.grid) <- CRS("+init=epsg:28992")
meuse.grid <- as(meuse.grid,"SpatialPixelsDataFrame")
plot(meuse.grid)
points(meuse,pch=20)

# Here is the idw with a power of 2.5
cdIDW <- idw(formula=cadmium~1,locations=meuse,newdata=meuse.grid,idp=2.5)

## [inverse distance weighted interpolation]

Let’s look at this IDW

head(cdIDW@data)

##   var1.pred var1.var
## 1  6.860300       NA
## 2  8.118605       NA
## 3  7.188141       NA
## 4  6.390504       NA
## 5 10.163382       NA
## 6  8.736704       NA

What we really want here is the variance column.

spplot(cdIDW,"var1.pred")

Plot still looks pretty. So let’s try Krieging, which calculates and fits a statistical model to the variances.

Kriging

  cadVar <- variogram(cadmium~1, meuse)
  plot(cadVar,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Cadmium concentrations (ppm)")

  sph.model <- vgm(psill=12, model="Sph", range=1000, nugget=5)
  sph.fit <- fit.variogram(object = cadVar, model = sph.model)

Above we need to give this sph.model some initial conditions to fit this linear model. We can add these in by eye and hope that it fits welle enough. Gives a starting point for the solution to converge upon.

  plot(cadVar,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Cadmium concentrations (ppm)",model=sph.fit)

  cadSph <- krige(cadmium~1, meuse, meuse.grid, model = sph.fit)

## [using ordinary kriging]

  spplot(cadSph,"var1.pred")

My Work

Varience

Let’s look at the valuable varience plot that is provided by Kriging:

spplot(cadSph,"var1.var")

print(bubble(meuse, xcol="x",ycol="y",zcol="cadmium", maxsize = 2.5, main = "Cadmium concentrations (ppm)"))

No surprise, this map is the warmest in areas with few data points. Ideally you would take additional samples in these areas. Still, due to the high cost of metals testing, if these areas are not of particular interest you could forego this additional testing with the understanding that your results may be unreliable in these areas.

Fitting the variogram model - Spherical

Let’s play with the variogram parameters, try a few, and see the results. First a slight modification

  sph.model <- vgm(psill=10, model="Sph", range=1000, nugget=4)
  sph.fit <- fit.variogram(object = cadVar, model = sph.model)
  plot(cadVar,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Cadmium concentrations (ppm)",model=sph.fit)

  cadSph <- krige(cadmium~1, meuse, meuse.grid, model = sph.fit)

## [using ordinary kriging]

  spplot(cadSph,"var1.pred")

Now a big modification:

  sph.model <- vgm(psill=2, model="Sph", range=1500, nugget=6)
  sph.fit <- fit.variogram(object = cadVar, model = sph.model)
  plot(cadVar,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Cadmium concentrations (ppm)",model=sph.fit)

  cadSph <- krige(cadmium~1, meuse, meuse.grid, model = sph.fit)

## [using ordinary kriging]

  spplot(cadSph,"var1.pred")

Discussion

It appears that this method is fairly insensitive to a poor approximation of the fit parameters, as all inputs resolved to the same linear model, and therefore the same Kriging output.

Different Linear Models