################################################### ### chunk number 1: Rsettings ################################################### rm(list=ls(all=T)) require(geoR) par.ori <- par(no.readonly=T) ################################################### ### chunk number 2: ################################################### data(elevation) ################################################### ### chunk number 3: elevation-hist ################################################### with(elevation, hist(data, main="", xlab="elevation")) ################################################### ### chunk number 5: elevation-coords ################################################### par(mar=c(3.5,3,0.5,0.5), mgp=c(2,0.8,0), mfrow=c(1,2)) with(elevation, plot(coords[,1], data, xlab="W-E", ylab="elevation data", pch=20, cex=0.7)) lines(lowess(elevation$data~elevation$coords[,1])) with(elevation,plot(coords[,2], data, xlab="S-N", ylab="elevation data", pch=20, cex=0.7)) lines(with(elevation, lowess(data~coords[,2]))) par(par.ori) ################################################### ### chunk number 7: elevation-points ################################################### par(mar=c(2,2,0.2,0.2), mgp=c(8,0.6,0), mfrow=c(1,3)) points(elevation, cex.max=2.5) points(elevation, trend="1st", pt.div=2, abs=T, cex.max=2.5) points(elevation, trend="2nd", pt.div=2, abs=T, cex.max=2.5) par(par.ori) ################################################### ### chunk number 10: elevation-variog ################################################### par(mar=c(3.5,3,0.4,0.4), mgp=c(1.8,0.7,0), mfrow=c(1,2)) plot(variog(elevation, uvec=seq(0,5,by=0.5)), xlab="u", ylab=expression(bar(v)), type="b") res1.v <- variog(elevation, trend="1st", uvec=seq(0,5,by=0.5)) plot(res1.v,xlab="u", ylab=expression(bar(v)), type="b") res2.v <- variog(elevation, trend="2nd", uvec=seq(0,5,by=0.5)) lines(res2.v, type="b", lty=2) par(par.ori) ################################################### ### chunk number 13: mc ################################################### set.seed(231) mc1 <- variog.mc.env(elevation, obj=res1.v) par(mar=c(3.5,3,0.4,0.4), mgp=c(1.8,0.7,0), mfrow=c(1,2)) plot(res1.v, env=mc1, xlab="u", ylab=expression(bar(v))) mc2 <- variog.mc.env(elevation, obj=res2.v) plot(res2.v, env=mc2,xlab="u", ylab=expression(bar(v))) par(par.ori) ################################################### ### chunk number 14: ml ################################################### ml0 <- likfit(elevation, ini=c(3000, 2), cov.model="matern", kappa=1.5) ml0 ml1 <- likfit(elevation, trend="1st", ini=c(1300, 2), cov.model="matern", kappa=1.5) ml1 ################################################### ### chunk number 15: ml-others ################################################### #ml1 <- likfit(elevation, trend=~coords[,2], ini=c(1300, 2), cov.model="matern", kappa=1.5) #ml1 <- likfit(elevation, trend=~coords[,2]+coords[,1]+I(coords[,1]^2), ini=c(1300, 2), cov.model="matern", kappa=1.5) ml2 <- likfit(elevation, trend="2nd", ini=c(1000, 1.5), cov.model="matern", kappa=1.5) ################################################### ### chunk number 16: elevation-kriging-cte ################################################### locs <- pred_grid(c(0,6.3), c(0,6.3), by=0.1) #KC <- krige.control(type="sk", beta=827, cov.model="mat", kappa=1.5, cov.pars=c(5000, 1.5), nug=0) KC <- krige.control(type="sk", obj.mod=ml0) sk <- krige.conv(elevation, krige=KC, loc=locs) ################################################### ### chunk number 17: elevation-kriging-trend ################################################### #KC <- krige.control(type="sk", beta=c(908, -25), cov.model="mat", kappa=1.5, cov.pars=c(1300, 1), nug=0, trend.d=~elevation$coords[,2], trend.l=~locs[,2]) KCt <- krige.control(type="sk", obj.mod=ml1, trend.d="1st", trend.l="1st") skt <- krige.conv(elevation, krige=KCt, loc=locs) ################################################### ### chunk number 18: elevation-kriging-trend2 ################################################### KCt2 <- krige.control(type="sk", obj.model=ml2, trend.d="2nd", trend.l="2nd") skt2 <- krige.conv(elevation, krige=KCt2, loc=locs) ################################################### ### chunk number 20: elevation-kriging-fig ################################################### pred.lim <- range(c(sk$pred,skt$pred)) sd.lim <- range(sqrt(c(sk$kr,skt$kr))) par(mar=c(2,2,0.4,0.4), mgp=c(2.5,0.7,0), mfrow=c(2,2)) image(sk, col=gray(seq(1,0.1,l=51)), zlim=pred.lim) contour(sk, add=T, nlev=6) points(elevation, add=TRUE, cex.max=2) image(skt, col=gray(seq(1,0.1,l=51)), zlim=pred.lim) contour(skt, add=T, nlev=6) points(elevation, add=TRUE, cex.max=2) image(sk, value=sqrt(sk$krige.var), col=gray(seq(1,0.1,l=51)), zlim=sd.lim) contour(sk, value=sqrt(sk$krige.var), levels=seq(10,27,by=2), add=T) points(elevation$coords, pch='+') image(skt, value=sqrt(skt$krige.var), col=gray(seq(1,0.1,l=51)), zlim=sd.lim) contour(skt, value=sqrt(skt$krige.var), levels=seq(10,27,by=2), add=T) points(elevation$coords, pch='+') par(par.ori) ################################################### ### chunk number 29: plot-elevation-res ################################################### plot(elevation, low = TRUE, trend = ~coords[, 2], qt.col = 1) ################################################### ### chunk number 31: plotca20-area ################################################### data(ca20) plot(ca20, trend = ~area, qt.col = 1)