date()
## [1] "Wed Apr 23 12:26:40 2014"
The objective of this problem set is to have you demonstrate skill at using functions in the geoR (or gstat) package for creating an interpolated surface using the methods of kriging. #“{r data} require(geoR) library(geoR) library(rgdal) library(maptools) library(sp) library(gstat) data(parana) summary(parana) par(mar = c(3,3,1,1)) #v = variogram(data ~ 1, data = parana) ## cannot work in my computer plot(variog(parana, option='cloud'),pch='.',cex=0.5) a<-variog(parana, option='smooth', band=100) lines(a$u,a$v)
summary(a) plot(parana)
par(las = 1) plot(a$u, sqrt(a$v), xlim = c(0, 500), ylim = c(0, 200), xlab = "Lagged distance (h) [m]”, ylab = expression(paste(“Semivariance (”, v, “) [”, cm2, “]”)), las = 1, pch = 19) grid() points(a$u, a$v, pch=10) text(a$u/10, a$v, pos=1, labels=as.character(a$n.data), cex=.5) vmi = vgm(model = “Gau”, psill = 80, range = 500, nugget = 25)
#v.fit = fit.variogram(a, vmi) #v.fit parana.bin<-variog(parana,option=“bin”,estimator.type=“classical”,trend='cte') plot(parana.bin$u,parana.bin$v,pch=20,ylab=“variogramme”,xlab='h',ylim=c(0,max(parana.bin$v))) ini<-c(6000,600) cov.model<-“exp” cov.model2<-“gaussian” wls <- variofit(parana.bin, ini = ini, cov.model=cov.model,weights=“cressie”) wls <- variofit(parana.bin, ini = ini, cov.model=cov.model2) #ols <- variofit(parana.bin, ini = ini, cov.model=cov.model,weights=“equal”) #ml <- likfit(parana, method=“ML”,trend=“1st”, ini = ini, cov.model=cov.model) #lines(ols) lines(wls,lty=2) #lines(ml,lty=3) #legend(400, 400, legend=c(“MCO”,“MCP”,“MV”), lty=c(1,2,3)) ###choose gaussian model, cressie weight fix.kappa<-TRUE kappa<-0.5 ml2 <- likfit(parana,cov.model=cov.model2,nugget=25,ini=c(6000,600),method=“ML”,trend=“2nd”,fix.kappa=fix.kappa) ##kriging library(fields) ngridx<-200 ngridy<-200 xgrid<-seq(min(parana$borders[,1]),max(parana$borders[,1]),l=ngridx) ygrid<-seq(min(parana$borders[,2]),max(parana$borders[,2]),l=ngridy) cols<-gray(seq(1,0.1,l=30)) #Universal Kriging pred.grid <- expand.grid(x=xgrid,y=ygrid) # defines a grid, with (min,max,length of vector) of x coordinate, and (min, max, length of vector) of y coordinate.
krige.ml <- krige.conv(parana,loc=pred.grid,krige=krige.control(obj.m=ml2,type.krige='ok',trend.d=“2nd”,trend.l=“2nd”),borders=NULL) image.plot(xgrid,ygrid,matrix(krige.ml$predict,ngridx,ngridy),xlab='x',ylab='y',col =cols) points(parana, pch=20,cex.max=1,pt.divide='equal',add.to.plot=TRUE, lty=2) contour(krige.ml,add=T,nlevels=10) #“ ## cannot plot the kriging in my computer
library(geoR) #loads geoR into your R session
## Warning: package 'geoR' was built under R version 3.0.3
## Loading required package: sp
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 3.0.3
## --------------------------------------------------------------
## 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
## --------------------------------------------------------------
data(parana) #loads the data parana into your R session.
par(mar = c(3, 3, 1, 1))
plot(variog(parana, option = "cloud"), pch = ".", cex = 0.5)
## variog: computing omnidirectional variogram
a <- variog(parana, option = "smooth", band = 100)
## variog: computing omnidirectional variogram
lines(a$u, a$v)
par(mar = c(3, 3, 1, 1))
a1 <- variog(parana, option = "smooth", direction = 0, band = 100)
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
a2 <- variog(parana, option = "smooth", direction = pi/4, band = 100)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
a3 <- variog(parana, option = "smooth", direction = pi/2, band = 100)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
a4 <- variog(parana, option = "smooth", direction = 3 * pi/4, band = 100)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
n.o <- c("0", "45", "90", "135")
ymax <- max(c(a1$v, a2$v, a3$v, a4$v), na.rm = TRUE)
plot(a1$u, a1$v, type = "l", xlab = "distance", ylab = "semivariance", ylim = c(0,
ymax))
lines(a2$u, a2$v, lty = 2)
lines(a3$u, a3$v, lty = 3, lwd = 2)
lines(a4$u, a4$v, lty = 4)
legend(0, ymax, legend = c(substitute(a * degree, list(a = n.o[1])), substitute(a *
degree, list(a = n.o[2])), substitute(a * degree, list(a = n.o[3])), substitute(a *
degree, list(a = n.o[4])), expression()), lty = c(1:4), lwd = c(1, 1, 2,
1))
par(mar = c(3, 3, 1, 1), pch = 20)
a1 <- variog(parana, option = "smooth", direction = 0, band = 100, trend = "1st")
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
a2 <- variog(parana, option = "smooth", direction = pi/4, band = 100, trend = "1st")
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
a3 <- variog(parana, option = "smooth", direction = pi/2, band = 100, trend = "1st")
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
a4 <- variog(parana, option = "smooth", direction = 3 * pi/4, band = 100, trend = "1st")
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
n.o <- c("0", "45", "90", "135")
ymax <- max(c(a1$v, a2$v, a3$v, a4$v), na.rm = TRUE)
plot(a1, type = "l", xlab = "distance", ylab = "semivariance", ylim = c(0, ymax))
lines(a2$u, a2$v, lty = 2)
lines(a3$u, a3$v, lty = 3, lwd = 2)
lines(a4$u, a4$v, lty = 4)
legend(0, ymax, legend = c(substitute(a * degree, list(a = n.o[1])), substitute(a *
degree, list(a = n.o[2])), substitute(a * degree, list(a = n.o[3])), substitute(a *
degree, list(a = n.o[4])), expression()), lty = c(1:4), lwd = c(1, 1, 2,
1))
par(mar = c(3, 3, 1, 1), pch = 20)
a1 <- variog(parana, option = "smooth", direction = 0, band = 100, trend = "2nd")
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
a2 <- variog(parana, option = "smooth", direction = pi/4, band = 100, trend = "2nd")
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
a3 <- variog(parana, option = "smooth", direction = pi/2, band = 100, trend = "2nd")
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
a4 <- variog(parana, option = "smooth", direction = 3 * pi/4, band = 100, trend = "2nd")
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
n.o <- c("0", "45", "90", "135")
ymax <- max(c(a1$v, a2$v, a3$v, a4$v), na.rm = TRUE)
plot(a1, type = "l", xlab = "distance", ylab = "semivariance", ylim = c(0, ymax))
lines(a2$u, a2$v, lty = 2)
lines(a3$u, a3$v, lty = 3, lwd = 2)
lines(a4$u, a4$v, lty = 4)
legend(0, ymax, legend = c(substitute(a * degree, list(a = n.o[1])), substitute(a *
degree, list(a = n.o[2])), substitute(a * degree, list(a = n.o[3])), substitute(a *
degree, list(a = n.o[4])), expression()), lty = c(1:4), lwd = c(1, 1, 2,
1))
par(mfrow = c(1, 1), mar = c(3, 3, 1, 1))
parana.cloud <- variog(parana, option = "cloud")
## variog: computing omnidirectional variogram
parana.bin <- variog(parana, option = "bin", estimator.type = "classical", trend = "cte")
## variog: computing omnidirectional variogram
plot(parana.cloud, ylab = "semi-variogramme", xlab = "h", pch = ".", cex = 0.3)
points(parana.bin$u, parana.bin$v, pch = 20, col = 1, cex = 1.4)
abline(v = parana.bin$bins.lim)
text(parana.bin$uvec, parana.bin$v, label = parana.bin$n, pos = 3, offset = 1,
cex = 0.5)
par(mfrow = c(1, 1), mar = c(3, 3, 1, 1))
parana.bin <- variog(parana, option = "bin", estimator.type = "classical", trend = "cte")
## variog: computing omnidirectional variogram
parana.rob <- variog(parana, option = "bin", estimator.type = "modulus", trend = "cte")
## variog: computing omnidirectional variogram
plot(parana.bin$u, parana.bin$v, pch = 20, col = 1, ylab = "semi-variogramme",
xlab = "h", ylim = c(0, max(c(parana.bin$v, parana.rob$v))))
points(parana.rob$u, parana.rob$v, pch = 21, col = 1)
legend(70, 6000, legend = c("classique", "robuste"), col = c(1, 1), pch = c(20,
21))
par(mfrow = c(1, 1), mar = c(3, 3, 1, 1))
parana.bin <- variog(parana, option = "bin", estimator.type = "classical", trend = "1st")
## variog: computing omnidirectional variogram
plot(parana.bin$u, parana.bin$v, pch = 20, ylab = "semi-variogramme", xlab = "h",
ylim = c(0, max(parana.bin$v)))
ini <- c(1000, 70)
cov.model <- "exp"
wls <- variofit(parana.bin, ini = ini, cov.model = cov.model, weights = "cressie")
## variofit: covariance model used is exponential
## variofit: weights used: cressie
## variofit: minimisation function used: optim
ols <- variofit(parana.bin, ini = ini, cov.model = cov.model, weights = "equal")
## variofit: covariance model used is exponential
## variofit: weights used: equal
## variofit: minimisation function used: optim
ml <- likfit(parana, method = "ML", trend = "1st", ini = ini, cov.model = cov.model)
## kappa not used for the exponential correlation function
## Warning: argument "method" has changed and is now used as an argument to
## be passed to optim(). Use "lik.method" to define the likelihood method
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
lines(ols)
lines(wls, lty = 2)
lines(ml, lty = 3)
legend(400, 400, legend = c("MCO", "MCP", "MV"), lty = c(1, 2, 3))
cov.model <- "exp"
fix.kappa <- TRUE
kappa <- 0.5
ml1 <- likfit(parana, cov.model = cov.model, nugget = 25, ini = c(1000, 70),
method = "ML", trend = "1st", fix.kappa = fix.kappa, kappa = kappa)
## kappa not used for the exponential correlation function
## Warning: argument "method" has changed and is now used as an argument to
## be passed to optim(). Use "lik.method" to define the likelihood method
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
ml1.aniso <- likfit(parana, cov.model = cov.model, nugget = 25, ini = c(1000,
70), method = "ML", trend = "1st", fix.psiA = FALSE, fix.psiR = FALSE, fix.kappa = fix.kappa,
kappa = kappa)
## kappa not used for the exponential correlation function
## Warning: argument "method" has changed and is now used as an argument to
## be passed to optim(). Use "lik.method" to define the likelihood method
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
ml2 <- likfit(parana, cov.model = cov.model, nugget = 25, ini = c(100, 20),
method = "ML", trend = "2nd", fix.kappa = fix.kappa, kappa = kappa)
## kappa not used for the exponential correlation function
## Warning: argument "method" has changed and is now used as an argument to
## be passed to optim(). Use "lik.method" to define the likelihood method
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
xv.ml1 <- xvalid(parana, model = ml1)
## xvalid: number of data locations = 143
## xvalid: number of validation locations = 143
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
## xvalid: end of cross-validation
xv.ml1.aniso <- xvalid(parana, model = ml1.aniso)
## xvalid: number of data locations = 143
## xvalid: number of validation locations = 143
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
## xvalid: end of cross-validation
xv.ml2 <- xvalid(parana, model = ml2)
## xvalid: number of data locations = 143
## xvalid: number of validation locations = 143
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
## xvalid: end of cross-validation
eqnm.ml1 <- mean(xv.ml1$std.error^2)
eqnm.ml1.aniso <- mean(xv.ml1.aniso$std.error^2)
eqnm.ml2 <- mean(xv.ml2$std.error^2)
par(mar = c(3, 3, 0, 1))
par(mgp = c(2, 1, 0))
plot(xv.ml2, ask = FALSE, error = FALSE, std.error = TRUE, data.predicted = FALSE,
pp = FALSE, map = TRUE, histogram = FALSE, error.predicted = FALSE, error.data = FALSE,
pch = 20)
par(mar = c(3, 3, 0, 1))
par(mgp = c(2, 1, 0))
plot(xv.ml2, ask = FALSE, error = FALSE, std.error = TRUE, data.predicted = FALSE,
pp = FALSE, map = FALSE, histogram = FALSE, error.predicted = TRUE, error.data = FALSE,
pch = 20)
par(mar = c(3, 3, 1, 1), pch = 20)
nsim <- 40
parana.vario <- variog(parana, trend = "2nd")
## variog: computing omnidirectional variogram
parana.env <- variog.model.env(parana, obj.v = parana.vario, model.pars = ml2,
nsim = nsim)
## Warning: the condition has length > 1 and only the first element will be
## used
## variog.env: generating 40 simulations (with 143 points each) using the function grf
## Warning: .Random.seed not initialised. Creating it with by calling
## runif(1)
## variog.env: adding the mean or trend
## variog.env: computing the empirical variogram for the 40 simulations
## variog.env: computing the envelops
parana.vario <- variog(parana)
## variog: computing omnidirectional variogram
plot(parana.vario, env = parana.env)
par(mar = c(3, 3, 1, 1), pch = 20)
ngridx <- 100
ngridy <- 100
xgrid <- seq(min(parana$borders[, 1]), max(parana$borders[, 1]), l = ngridx)
ygrid <- seq(min(parana$borders[, 2]), max(parana$borders[, 2]), l = ngridy)
cols <- gray(seq(1, 0.1, l = 30))
# Universal Kriging
pred.grid <- expand.grid(x = xgrid, y = ygrid) # defines a grid, with (min,max,length of vector) of x coordinate, and (min, max, length of vector) of y coordinate.
krige.ml <- krige.conv(parana, loc = pred.grid, krige = krige.control(obj.m = ml2,
type.krige = "ok", trend.d = "2nd", trend.l = "2nd"), borders = NULL)
## krige.conv: model with mean given by a 2nd order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood
image(krige.ml, loc = pred.grid, col = gray(seq(1, 0.1, l = 30)), xlab = "O-E",
ylab = "S-N", bor = parana$borders)
contour(krige.ml, loc = pred.grid, values = krige.ml$pred, xlab = "O-E", ylab = "S-N",
bor = parana$borders, add = TRUE)
points(parana$coords, pch = 20, cex = 0.7)
par(mar = c(3, 3, 1, 1), pch = 20)
image(krige.ml, loc = pred.grid, col = gray(seq(1, 0.1, l = 30)), xlab = "O-E",
ylab = "S-N", bor = parana$borders)
contour(krige.ml, loc = pred.grid, values = sqrt(krige.ml$krige.var), xlab = "O-E",
ylab = "S-N", bor = parana$borders, add = TRUE, nlevels = 15)
points(parana$coords, pch = 20, cex = 0.7)
par(mar = c(3, 3, 1, 1), pch = 20)
library(fields)
## Warning: package 'fields' was built under R version 3.0.3
## Loading required package: spam
## Warning: package 'spam' was built under R version 3.0.3
## Loading required package: grid
## Spam version 0.41-0 (2014-02-26) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
##
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
##
## Loading required package: maps
ngridx <- 200
ngridy <- 200
xgrid <- seq(min(parana$borders[, 1]), max(parana$borders[, 1]), l = ngridx)
ygrid <- seq(min(parana$borders[, 2]), max(parana$borders[, 2]), l = ngridy)
cols <- gray(seq(1, 0.1, l = 30))
# Universal Kriging
pred.grid <- expand.grid(x = xgrid, y = ygrid) # defines a grid, with (min,max,length of vector) of x coordinate, and (min, max, length of vector) of y coordinate.
krige.ml <- krige.conv(parana, loc = pred.grid, krige = krige.control(obj.m = ml2,
type.krige = "ok", trend.d = "2nd", trend.l = "2nd"), borders = NULL)
## krige.conv: model with mean given by a 2nd order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood
image.plot(xgrid, ygrid, matrix(krige.ml$predict, ngridx, ngridy), xlab = "x",
ylab = "y", col = cols)
points(parana, pch = 20, cex.max = 1, pt.divide = "equal", add.to.plot = TRUE,
lty = 2)
contour(krige.ml, add = T, nlevels = 10)
par(mar = c(3, 3, 1, 1), pch = 20)
image.plot(xgrid, ygrid, matrix(krige.ml$predict, ngridx, ngridy), xlab = "x",
ylab = "y", col = cols)
points(parana, pch = 20, cex.max = 1, pt.divide = "equal", add.to.plot = TRUE,
lty = 2)
contour(krige.ml, values = sqrt(krige.ml$krige.var), add = T, nlevels = 10)
