########################## ## An example of kriging## ########################## par(mfrow=c(1,1)) vcloud(y=paws\$Soil8, locs=paws[,c("X", "Y")]) par(mfrow=c(1,2)) plot(paws[, c("X", "Y")]) hist(paws\$Soil8, main="", xlab="Soil temperature") par(mfrow=c(1,2)) y2=dist(paws\$Soil8) d=dist(paws[,c("X", "Y")]) plot(d, (y2^2)/2, xlab="distance", ylab="Semivariance", main="(a)") plot(d[d<100], (y2[d<100]^2)/2, xlab="distance", ylab="Semivariance", main="(b)") plot(variog(coords=paws[, c("X", "Y")], data=paws[, "Soil8"], max.dist=100)) lines(soil.likfit) # REML estimates soil.likfit=likfit(coords = paws[, c("X", "Y")], data = paws\$Soil8, ini.cov.pars=c(13,4.4), nugget=0.3, method="REML") #> soil.likfit\$nugget #[1] 1.279161 #> soil.likfit\$cov.par #[1] 10.593613 6.238187 par(mfrow=c(1,2)) newx=seq(min(paws\$X), max(paws\$X),length=50) newy=seq(min(paws\$Y), max(paws\$Y), length=50) predloc=cbind(rep(newx, length=length(newy)), rep(newy, each=length(newx))) plot(paws[,c("X", "Y")], main="(a)") points(predloc[,1], predloc[,2], pch="+", cex=0.5) ## Ordinary Kriging soil.k=krige.conv(coords=paws[,c("X", "Y")], data=paws[, "Soil8"], locations=predloc, krige=krige.control(cov.model="exponential", cov.pars=c(sigmasq=10.59, phi=6.24), nugget=1.27 ))\$predict soil.mat=matrix(soil.k, ncol=length(newy)) persp(x=newx,y=newy,z=soil.mat, xlab="X", ylab="Y", zlab="Kriging surface", theta=45, phi=35, r=5, expand=0.5, axes=T, ticktype="detailed") soil.k=krige.conv(coords=paws[,c("X", "Y")], data=paws[, "Soil8"], locations=predloc, krige=krige.control(cov.model="exponential", cov.pars=c(sigmasq=10.59, phi=6.24), nugget=1.27 ))\$krige.var soil.mat=matrix(soil.k, ncol=length(newy)) persp(x=newx,y=newy,z=soil.mat, xlab="X", ylab="Y", zlab="Kriging surface", theta=45, phi=35, r=5, expand=0.5, axes=T, ticktype="detailed") filled.contour(x=newx,y=newy,z=soil.mat, xlab="X", ylab="Y", color.palette=gray.colors) filled.contour(x=seq(min(newx), max(newx), length=48),y=seq(min(newy), max(newy), length=48),z=soil.mat, xlab="X", ylab="Y", color.palette=terrain.colors) vertices=paws[chull(paws[,c("X", "Y")]), c("X", "Y")] plot(paws[,c('X', 'Y')]) lines(vertices) lines(paws[c(2,34), c("X", "Y")]) # You need to load geoR or sp to call the function point.in.polygon pip=point.in.polygon(predloc[,1], predloc[,2], vertices[,1], vertices[,2]) newloc=predloc[pip>0,] plot(paws[,c("X", "Y")], main="(b)") points(newloc[,1], newloc[,2], pch="+", cex=0.5) soil.krige=rep(NA, nrow(predloc)) soil.krige[pip>0]=krige.conv(coords=paws[,c("X", "Y")], data=paws[, "Soil8"], locations=newloc, krige=krige.control(cov.model="exponential", cov.pars=c(sigmasq=10.59, phi=6.24), nugget=1.27 ))\$predict soil.mat=matrix(soil.krige, ncol=length(newy)) persp(x=newx,y=newy,z=soil.mat, xlab="X", ylab="Y", zlab="Kriging surface", theta=45, phi=35, r=5, expand=0.5, axes=T, ticktype="detailed") filled.contour(x=newx,y=newy,z=soil.mat, xlab="X", ylab="Y", color.palette=gray.colors) plot(paws[,c("X", "Y")]) vertices=locator(type="l") vertices=cbind(vertices\$x, vertices\$y) > vert2 #\$x # [1] 275.9409 265.8782 267.3157 275.3659 292.3288 310.1543 327.9797 332.0048 # [9] 331.7173 317.0545 301.8166 276.2284 # #\$y # [1] 5133.105 5120.694 5098.186 5080.095 5095.872 5091.034 5091.875 5096.713 # [9] 5107.442 5119.222 5131.422 5133.315 plot(paws[,c("X", "Y")], main="(a)") lines(vert2\$x, vert2\$y) newloc2=predloc[point.in.polygon(predloc[,1], predloc[,2], vert2\$x, vert2\$y)>0,] plot(paws[,c("X", "Y")], main="(b)") points(newloc2[,1], newloc2[,2], pch="+", cex=0.5)