[R-sig-Geo] Kriging - model selection

Tue May 19 21:32:09 CEST 2015

Dear all,

I would like to perform interpolation of my data, using kriging. I am
pretty new to this field, I would appreciate any advice.

After reading various sources about how to do this kind of analysis in R, i
came out with the code below. What I am struggling with, is the selection
of the "correct" variogram model and parameters which would best fit my
data.

So far, I am most convinced with the linear model using "gstat" package (as
there is a linear correlation between value and longitude), although I am
not sure with the "range" parameter.

I tried two packages, "kriging" and "gstat", here is the code I used.

I welcome any recommendation about which model and parameters to use.
Thank you,

Zuzana

#####################
###### KRIGING ######
#####################

> dput(ddd)
 structure(list(value = c(-15.63, -15.85, -15.66, -15.59, -15.75,
                         -15.67, -15.92, -16.1, -15.94, -15.86, -16.34,
-15.79, -15.74,
                         -16.35, -16.32, -16.34, -16.93, -16.33, -16.09,
-16.15, -16.07,
                         -16.32, -16.14, -16.11, -16.42, -16.06, -16.09,
-15.67), lat = c(7.1840870481035,

7.99259507438396, 11.9797581011895, 7.49274296840517, 9.67854016040156,

12.2741332643094, 8.21606349541263, 7.79878658231986, 7.2680142040734,

7.79600014665774, 6.88856759081784, 2.71471364119849, 3.98276218734202,

9.14251390940468, 7.40818160773457, 6.74255907653643, 17.2606233919273,

13.5775247347178, 15.1123203521451, 10.5186702866744, 7.49287138877741,

11.6415159877501, 7.37141449045903, 10.4795385344797, 11.7311433605838,

7.69341033361702, 5.83759526736964, 6.99439199170752), lon =
c(-37.9137449718093,

-35.3295205743433, -35.7977158093952, -41.4601218055409, -38.6729750935733,

-39.6048692841313, -34.5682210607579, -36.3718007962673, -33.1684288611563,

-40.9943844182636, -30.76067153812, -34.1185676916864, -37.1053996472096,

-32.2766293664428, -33.7449877248086, -34.5910032927336, -23.8938102400732,

-37.0736715552604, -33.9073649895551, -34.0622219110843, -34.656850099908,

-32.4980341080834, -40.5035367177989, -34.3580298956927, -33.1636537431894,

-39.8857782758382, -36.3770793693717, -45.0208499705794)), .Names =
c("value",

"lat", "lon"), row.names = c(21L, 53L, 98L, 131L, 171L, 232L,

263L, 312L, 344L, 383L, 413L, 445L, 511L, 538L, 568L, 629L, 669L,

708L, 729L, 800L, 837L, 853L, 882L, 943L, 988L, 1037L, 1057L,

1095L), class = "data.frame")

##############################
library(kriging)
library(fields)
##############################

kriged.ddd <- kriging(ddd$lon, ddd$lat, ddd$value,  pixels=300)

plot(kriged.ddd$semivariogram)

r <- range(kriged.ddd$map$pred)

## plot
image(kriged.ddd, xlim = extendrange(ddd$lon), ylim = extendrange(ddd$lat),
      zlim = r, col=topo.colors(100), add=F)
points(ddd$lon, ddd$lat, pch=21, col="white", bg="grey20" )
image.plot(kriged.ddd, col=topo.colors(100), legend.only=TRUE, zlim=r,
horizontal=T, legend.shrink=0.3)

##############################
library(gstat)
##############################

coordinates(ddd) = ~lon+lat
plot(ddd)

x.range <- as.numeric(range(ddd at coords[,1])) # lon
y.range <- as.numeric(range(ddd at coords[,2])) # lat

## expand grid
grd <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=0.5),
                   y=seq(from=y.range[1], to=y.range[2], by=0.5))

# grd2 <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=1),
#                    y=seq(from=y.range[1], to=y.range[2], by=1))

## convert grid to SpatialPixel class
coordinates(grd) <- ~ x+y
gridded(grd) <- TRUE

# coordinates(grd2) <- ~ x+y
# gridded(grd2) <- TRUE

## test the grid and points
plot(grd, cex=1.5)
points(ddd, pch=1, col="red", cex=1)

## construction of a semivariogram model
variogcloud<-variogram(value~1, locations=ddd, data=ddd, cloud=TRUE)
plot(variogcloud)

semivariog <- variogram(value~1, locations=ddd, data=ddd)
plot(semivariog)

## estimation of semivariogram parameters and fit the model

# linear
model.variog.lin <- vgm(psill=1, model="Lin", nugget=0.02) ## range???
model.variog.lin
plot(semivariog, model.variog.lin, plot.numbers=T)
fit.variog.lin <- fit.variogram(semivariog, model.variog.lin)
fit.variog.lin
plot(semivariog, fit.variog.lin, plot.numbers=T)

# spherical
model.variog.sph <- vgm(psill=0.05, model="Sph", nugget=0.01, range=4)
model.variog.sph
plot(semivariog, model.variog.sph, plot.numbers=T)
fit.variog.sph <- fit.variogram(semivariog, model.variog.sph)
fit.variog.sph
plot(semivariog, fit.variog.sph, plot.numbers=T)

# exponencial
model.variog.exp <- vgm(psill=0.05,model="Exp", nugget=0.01, range=3)
model.variog.exp
plot(semivariog, model.variog.exp, plot.numbers=T)
fit.variog.exp <- fit.variogram(semivariog, model.variog.exp)
fit.variog.exp
plot(semivariog, fit.variog.exp, plot.numbers=T)

####### plots

## linear semivariogram model
# krig.lin <- krige(formula=value ~ 1, locations=ddd, newdata=grd,
model=model.variog.lin)
krig.lin <- krige(formula=value ~ 1, locations=ddd, newdata=grd,
model=fit.variog.lin)
krig.output.lin=as.data.frame(krig.lin)
names(krig.output.lin)[1:3] <- c("long","lat","var1.pred.lin")
plot.lin <- ggplot(data=krig.output.lin,aes(x=long,y=lat))
layer.lin <- c(geom_tile(data=krig.output.lin, aes(fill=var1.pred.lin)))
plot.lin+
  layer.lin+
    scale_fill_gradientn(colours=topo.colors(100),
limits=range(range(krig.output.lin$var1.pred.lin), range(ddd$value)))+
  coord_equal()+
  theme_bw()+
  geom_point(data=data.frame(ddd), aes(x=lon, y=lat),  fill="black",
colour="black", shape=21, size=4)+
  geom_point(data=data.frame(ddd), aes(x=lon, y=lat,  fill=value),
colour="black", shape=21, size=3.5)

## spherical semivariogram model
krig.sph <- krige(formula=value ~ 1, locations=ddd, newdata=grd,
model=model.variog.sph)
# krig.sph <- krige(formula=value ~ 1, locations=ddd, newdata=grd,
model=fit.variog.sph)
krig.output.sph=as.data.frame(krig.sph)
names(krig.output.sph)[1:3] <- c("long","lat","var1.pred.sph")
plot.sph <- ggplot(data=krig.output.sph,aes(x=long, y=lat))
layer.sph <- c(geom_tile(data=krig.output.sph, aes(fill=var1.pred.sph)))
plot.sph+
  layer.sph+
  scale_fill_gradientn(colours=topo.colors(100),
limits=range(range(krig.output.sph$var1.pred.sph), range(ddd$value)))+
  coord_equal()+
  theme_bw()+
  geom_point(data=data.frame(ddd), aes(x=lon, y=lat),  fill="black",
colour="black", shape=21, size=4)+
  geom_point(data=data.frame(ddd), aes(x=lon, y=lat,  fill=value),
colour="black", shape=21, size=3.5)

## exponential semivariogram model
krig.exp <- krige(formula=value ~ 1, locations=ddd, newdata=grd,
model=model.variog.exp)
# krig.exp <- krige(formula=value ~ 1, locations=ddd, newdata=grd,
model=fit.variog.exp)
krig.output.exp=as.data.frame(krig.exp)
names(krig.output.exp)[1:3] <- c("long","lat","var1.pred.exp")
plot.exp <- ggplot(data=krig.output.exp,aes(x=long,y=lat))
layer.exp <- c(geom_tile(data=krig.output.exp,aes(fill=var1.pred.exp)))
plot.exp+
  layer.exp+
 scale_fill_gradientn(colours=topo.colors(100),
limits=range(range(krig.output.exp$var1.pred.exp), range(ddd$value)))+
  coord_equal()+
  theme_bw()+
  geom_point(data=data.frame(ddd), aes(x=lon, y=lat),  fill="black",
colour="black", shape=21, size=4)+
  geom_point(data=data.frame(ddd), aes(x=lon, y=lat,  fill=value),
colour="black", shape=21, size=3.5)

<zuzanazajkova at ub.edu>

	[[alternative HTML version deleted]]