# Computing Issues with Kriging

I am having some issues with Kriging in R, and I was looking for some idea where I am going wrong.

From what I can tell, I done a decent job removing the trend, and I believe my transformed data is more or less isotropic. So I cannot figure out why my prediction SE is so high (see plot below). There are few places where I think I may have tripped up...

1. I didn't understand the R-code, and I just put some argument in the wrong place.
2. I am miss informed, and I am not Kriging inappropriately.
3. I took the log transformation of the response and that is somehow messed up my mapping of the SE.

Thanks for the help!

library(geoR)
library(gstat)
library(scatterplot3d)
library(rgl)
library(sp)
library(lattice)
library(akima)
library(grid)

############ Begin Functions ###############
# This function is written by Dr. Jon Graham
#(http://www.math.umt.edu/graham/stat544/)

rotateaxis <- function(xy, theta)
{
pimult <- (theta * 2 * pi)/360
newx <- c(cos(pimult), sin(pimult))
newy <- c( - sin(pimult), cos(pimult))
XY <- as.matrix(xy) %*% cbind(newx, newy)
as.data.frame(XY)
}

# function written by Dr. Jon Graham (http://www.math.umt.edu/graham/stat544/)

rose <- function(dat.var,critgam){
numcases <- length(dat.var$dist) numdirec <- length(unique(dat.var$dir.hor))
rose.len <- rep(0,numdirec)
rose.azi <- rep(0,numdirec)
j <- 1
for (i in 1:(numcases-1)){
prod <- (dat.var$gamma[i] - critgam)*(dat.var$gamma[i+1] - critgam)
if (dat.var$dir.hor[i] != dat.var$dir.hor[i+1]) prod <- 0
if (j > 1){
if (rose.azi[j-1] == dat.var$dir.hor[i]) prod <- 0 } if (prod < 0){ rose.len[j] <- dat.var$dist[i] + (critgam - dat.var$gamma[i])* (dat.var$dist[i+1] - dat.var$dist[i])/ (dat.var$gamma[i+1] - dat.var$gamma[i]) rose.azi[j] <- as.numeric(unique(dat.var$dir.hor)[j])
j <- j + 1
}
}
rose.azi <- rose.azi*pi/180
dscale <- max(rose.len)
rose.len <- rose.len/dscale
plot(0,0,type="n",axes=F,xlim=c(-1.1,1.1),ylim=c(-1.1,1.1),
xlab="E-W Direction",ylab="",cex.lab=1.6)
for (j in 1:(length(rose.len))){
segments(0,0,(rose.len[j]*sin(rose.azi[j])),(rose.len[j]*cos(rose.azi[j])))
text((rose.len[j]*sin(rose.azi[j]) + (.1*sin(rose.azi[j]))),
(rose.len[j]*cos(rose.azi[j]) + (.1*cos(rose.azi[j]))),
(format(round(rose.len[j]*dscale,2))),cex=1.5)
segments(0,0,-(rose.len[j]*sin(rose.azi[j])),
-(rose.len[j]*cos(rose.azi[j])))
}
maxx <- max(rose.len*sin(rose.azi))
text(-(maxx+.25),0,"N-S Direction",srt=90,crt=90,cex=1.6)
title("Rose Diagram",srt=0,crt=0,cex.main=2.0)
print("Rose Lengths")
rose.len*dscale
}

## written by Dr. Jon Graham, website: http://www.math.umt.edu/graham/stat544/panel.gamma0.r

panelgamma0 <- function(x,y,gamma0=gamma0,span=2/3, ...){
lofit <- loess.smooth(x,y,span=span)
panel.xyplot(x,y,...)
panel.xyplot(lofit$x,lofit$y,type="l")
nump <- round(length(lofit$x)/2) dist0 <- approx(lofit$y[1:nump],lofit$x[1:nump],xout=gamma0)$y
panel.segments(0,gamma0,dist0,gamma0)
panel.segments(dist0,0,dist0,gamma0)
parusr <- par()$usr panel.text(max(x),min(y),paste("d0=",format(round(dist0,4)),sep=""),adj=1) } ############### End Functions ########## ### Read data ### air <-structure(list(easting = c(256L, 256L, 256L, 275L, 275L, 275L, 275L, 275L, 275L, 275L, 275L, 295L, 295L, 295L, 295L, 295L, 295L, 295L, 295L, 295L, 295L, 295L, 295L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 314L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 333L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 353L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 372L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 392L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 411L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 430L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 450L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 469L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 489L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 508L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 527L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 547L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 566L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 586L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 605L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 624L, 644L, 644L, 644L, 644L, 644L, 644L, 644L, 644L, 644L, 644L, 663L, 663L, 663L, 663L, 663L, 663L, 663L, 683L, 683L, 683L, 683L, 683L, 683L, 683L, 683L, 683L, 683L, 683L, 702L, 702L, 702L, 702L, 702L, 702L, 702L, 702L, 721L, 721L, 721L, 721L, 721L, 721L, 721L, 741L, 741L, 741L, 741L, 741L, 741L, 760L, 760L, 760L, 760L, 760L), northing = c(258L, 276L, 293L, 258L, 276L, 293L, 311L, 329L, 346L, 539L, 557L, 276L, 293L, 311L, 329L, 346L, 434L, 452L, 469L, 487L, 504L, 539L, 575L, 135L, 241L, 258L, 276L, 293L, 311L, 329L, 399L, 416L, 434L, 452L, 469L, 504L, 522L, 539L, 575L, 592L, 100L, 118L, 135L, 153L, 223L, 241L, 258L, 276L, 293L, 311L, 381L, 399L, 416L, 434L, 469L, 487L, 504L, 522L, 539L, 557L, 575L, 592L, 65L, 82L, 100L, 118L, 135L, 153L, 170L, 223L, 241L, 258L, 276L, 293L, 311L, 381L, 399L, 416L, 434L, 452L, 469L, 487L, 504L, 522L, 539L, 557L, 575L, 65L, 82L, 100L, 135L, 153L, 170L, 223L, 241L, 258L, 276L, 293L, 311L, 329L, 364L, 381L, 399L, 416L, 434L, 452L, 469L, 487L, 504L, 522L, 539L, 698L, 65L, 82L, 118L, 135L, 153L, 170L, 188L, 205L, 223L, 241L, 258L, 276L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 434L, 452L, 487L, 504L, 522L, 539L, 575L, 592L, 610L, 645L, 65L, 82L, 100L, 135L, 153L, 170L, 188L, 205L, 223L, 241L, 258L, 276L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 434L, 487L, 504L, 522L, 539L, 557L, 575L, 627L, 663L, 680L, 65L, 82L, 118L, 135L, 153L, 170L, 188L, 205L, 223L, 241L, 258L, 276L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 452L, 469L, 487L, 504L, 522L, 539L, 557L, 680L, 65L, 100L, 118L, 135L, 153L, 170L, 188L, 205L, 223L, 241L, 258L, 276L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 434L, 452L, 469L, 487L, 504L, 522L, 539L, 610L, 645L, 663L, 82L, 153L, 223L, 241L, 258L, 276L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 452L, 469L, 487L, 504L, 522L, 539L, 557L, 575L, 592L, 680L, 698L, 82L, 100L, 258L, 276L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 434L, 452L, 469L, 487L, 504L, 522L, 557L, 575L, 645L, 663L, 680L, 698L, 135L, 258L, 276L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 452L, 469L, 487L, 504L, 522L, 539L, 557L, 575L, 592L, 610L, 627L, 645L, 663L, 680L, 698L, 135L, 205L, 223L, 258L, 276L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 434L, 452L, 469L, 487L, 522L, 539L, 557L, 575L, 592L, 610L, 627L, 645L, 663L, 680L, 698L, 188L, 205L, 258L, 276L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 434L, 469L, 487L, 522L, 539L, 557L, 575L, 610L, 627L, 645L, 663L, 680L, 698L, 258L, 276L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 434L, 452L, 469L, 539L, 557L, 627L, 645L, 663L, 680L, 698L, 170L, 205L, 223L, 258L, 276L, 311L, 346L, 364L, 381L, 399L, 416L, 434L, 452L, 469L, 487L, 504L, 522L, 539L, 592L, 610L, 627L, 645L, 698L, 223L, 241L, 258L, 276L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 434L, 452L, 469L, 487L, 504L, 522L, 539L, 557L, 575L, 592L, 645L, 680L, 698L, 153L, 223L, 241L, 258L, 293L, 311L, 346L, 364L, 381L, 399L, 416L, 469L, 487L, 504L, 522L, 539L, 557L, 575L, 698L, 153L, 311L, 329L, 346L, 364L, 381L, 399L, 469L, 522L, 539L, 153L, 311L, 364L, 381L, 504L, 610L, 645L, 293L, 311L, 329L, 346L, 364L, 381L, 399L, 416L, 434L, 504L, 539L, 346L, 364L, 381L, 399L, 452L, 522L, 557L, 627L, 416L, 434L, 452L, 469L, 522L, 539L, 592L, 188L, 241L, 364L, 452L, 504L, 557L, 487L, 539L, 575L, 592L, 610L), ammonia = c(103L, 175L, 95L, 297L, 370L, 369L, 121L, 101L, 137L, 5L, 1L, 67L, 13L, 299L, 163L, 37L, 13L, 13L, 9L, 10L, 6L, 2L, 2L, 12L, 28L, 230L, 112L, 160L, 346L, 6L, 31L, 19L, 22L, 28L, 15L, 1L, 1L, 2L, 6L, 6L, 4L, 60L, 143L, 112L, 248L, 119L, 102L, 394L, 331L, 108L, 25L, 23L, 17L, 32L, 5L, 6L, 10L, 4L, 4L, 3L, 10L, 1L, 9L, 82L, 92L, 78L, 84L, 119L, 7L, 220L, 376L, 190L, 169L, 158L, 42L, 73L, 24L, 48L, 16L, 8L, 8L, 2L, 24L, 17L, 23L, 19L, 20L, 29L, 66L, 78L, 96L, 91L, 78L, 373L, 486L, 342L, 271L, 404L, 196L, 291L, 193L, 500L, 95L, 73L, 39L, 42L, 4L, 28L, 58L, 33L, 20L, 1L, 27L, 62L, 80L, 97L, 63L, 58L, 42L, 99L, 285L, 174L, 124L, 182L, 371L, 751L, 745L, 458L, 380L, 185L, 96L, 48L, 34L, 10L, 60L, 48L, 29L, 24L, 70L, 5L, 5L, 5L, 40L, 63L, 47L, 45L, 61L, 77L, 195L, 178L, 279L, 275L, 208L, 181L, 264L, 177L, 279L, 228L, 326L, 175L, 49L, 15L, 18L, 60L, 17L, 88L, 58L, 70L, 25L, 10L, 5L, 5L, 3L, 20L, 43L, 41L, 49L, 82L, 114L, 192L, 240L, 259L, 244L, 215L, 171L, 151L, 269L, 383L, 317L, 239L, 105L, 117L, 50L, 136L, 24L, 50L, 75L, 75L, 5L, 13L, 18L, 40L, 12L, 5L, 32L, 85L, 19L, 79L, 132L, 138L, 312L, 322L, 348L, 251L, 380L, 269L, 242L, 250L, 200L, 46L, 224L, 164L, 137L, 150L, 110L, 75L, 10L, 10L, 10L, 20L, 33L, 27L, 116L, 297L, 290L, 156L, 334L, 280L, 293L, 283L, 324L, 319L, 276L, 223L, 211L, 182L, 150L, 100L, 100L, 75L, 50L, 50L, 10L, 10L, 20L, 30L, 178L, 446L, 247L, 113L, 288L, 261L, 301L, 260L, 274L, 280L, 248L, 199L, 199L, 199L, 166L, 150L, 100L, 100L, 25L, 25L, 50L, 25L, 10L, 108L, 303L, 145L, 151L, 188L, 252L, 244L, 237L, 274L, 272L, 203L, 199L, 199L, 183L, 150L, 150L, 150L, 100L, 150L, 75L, 75L, 75L, 75L, 75L, 50L, 30L, 125L, 36L, 94L, 111L, 122L, 185L, 261L, 232L, 205L, 236L, 240L, 247L, 235L, 209L, 180L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 100L, 100L, 75L, 50L, 17L, 17L, 150L, 73L, 30L, 126L, 132L, 208L, 288L, 261L, 241L, 245L, 245L, 245L, 245L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 100L, 50L, 113L, 81L, 120L, 142L, 182L, 214L, 224L, 202L, 250L, 242L, 245L, 245L, 150L, 150L, 150L, 150L, 150L, 75L, 50L, 20L, 14L, 14L, 70L, 101L, 70L, 181L, 131L, 194L, 185L, 247L, 248L, 245L, 245L, 245L, 245L, 216L, 150L, 150L, 150L, 150L, 150L, 100L, 83L, 24L, 75L, 43L, 15L, 100L, 183L, 177L, 252L, 281L, 252L, 257L, 219L, 196L, 245L, 245L, 245L, 245L, 233L, 150L, 150L, 100L, 50L, 50L, 100L, 5L, 33L, 16L, 9L, 38L, 118L, 190L, 202L, 202L, 313L, 85L, 215L, 184L, 184L, 202L, 216L, 150L, 50L, 50L, 5L, 117L, 96L, 116L, 172L, 184L, 39L, 175L, 200L, 150L, 10L, 83L, 87L, 86L, 150L, 5L, 5L, 66L, 66L, 19L, 37L, 133L, 96L, 124L, 86L, 91L, 30L, 100L, 56L, 124L, 129L, 129L, 219L, 91L, 70L, 5L, 118L, 146L, 183L, 211L, 13L, 116L, 40L, 5L, 10L, 70L, 138L, 250L, 133L, 157L, 212L, 128L, 128L, 28L), nitroxides = c(3L, 3L, 1L, 3L, 57L, 20L, 54L, 35L, 4L, 6L, 2L, 2L, 3L, 4L, 3L, 5L, 11L, 6L, 3L, 1L, 3L, 2L, 1L, 36L, 1L, 22L, 108L, 395L, 376L, 1L, 37L, 3L, 1L, 9L, 3L, 1L, 2L, 2L, 25L, 10L, 1L, 30L, 6L, 201L, 17L, 8L, 4L, 96L, 275L, 80L, 5L, 4L, 38L, 11L, 8L, 3L, 5L, 5L, 5L, 3L, 38L, 15L, 2L, 147L, 9L, 1L, 9L, 50L, 4L, 18L, 36L, 15L, 96L, 145L, 7L, 4L, 5L, 41L, 6L, 7L, 8L, 4L, 9L, 13L, 8L, 7L, 2L, 1L, 1L, 5L, 7L, 1L, 11L, 57L, 26L, 125L, 18L, 169L, 255L, 36L, 21L, 49L, 40L, 30L, 24L, 32L, 1L, 20L, 13L, 11L, 7L, 26L, 15L, 15L, 15L, 74L, 3L, 49L, 10L, 49L, 34L, 32L, 207L, 86L, 242L, 557L, 358L, 136L, 159L, 84L, 24L, 13L, 15L, 50L, 25L, 27L, 11L, 5L, 5L, 2L, 9L, 24L, 20L, 5L, 30L, 3L, 10L, 2L, 34L, 21L, 21L, 26L, 28L, 24L, 69L, 328L, 325L, 163L, 52L, 10L, 9L, 3L, 1L, 45L, 30L, 28L, 2L, 2L, 2L, 24L, 30L, 2L, 2L, 14L, 5L, 1L, 151L, 4L, 24L, 43L, 14L, 78L, 43L, 63L, 163L, 378L, 64L, 284L, 184L, 108L, 2L, 2L, 11L, 11L, 15L, 18L, 2L, 5L, 2L, 4L, 5L, 43L, 85L, 30L, 98L, 7L, 12L, 53L, 49L, 74L, 106L, 109L, 124L, 139L, 503L, 283L, 2L, 25L, 25L, 3L, 5L, 20L, 1L, 5L, 40L, 1L, 2L, 2L, 2L, 1L, 10L, 54L, 58L, 94L, 121L, 39L, 139L, 24L, 180L, 307L, 34L, 161L, 36L, 12L, 3L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 15L, 14L, 2L, 10L, 117L, 175L, 88L, 38L, 58L, 17L, 225L, 541L, 88L, 108L, 8L, 5L, 15L, 15L, 2L, 1L, 12L, 12L, 15L, 2L, 24L, 14L, 1L, 68L, 114L, 46L, 40L, 45L, 71L, 60L, 111L, 59L, 15L, 2L, 100L, 15L, 14L, 2L, 12L, 149L, 24L, 7L, 2L, 2L, 2L, 5L, 19L, 63L, 2L, 27L, 14L, 85L, 48L, 26L, 55L, 71L, 106L, 112L, 4L, 9L, 1L, 2L, 5L, 2L, 1L, 36L, 36L, 24L, 43L, 9L, 12L, 2L, 24L, 2L, 63L, 4L, 3L, 97L, 83L, 3L, 8L, 11L, 17L, 35L, 24L, 5L, 9L, 1L, 7L, 1L, 2L, 15L, 38L, 5L, 2L, 24L, 5L, 24L, 2L, 58L, 58L, 43L, 7L, 7L, 32L, 33L, 6L, 43L, 5L, 5L, 3L, 76L, 33L, 7L, 24L, 24L, 5L, 12L, 12L, 2L, 9L, 8L, 35L, 30L, 7L, 20L, 6L, 1L, 6L, 4L, 1L, 1L, 1L, 1L, 2L, 22L, 1L, 24L, 2L, 1L, 2L, 2L, 58L, 16L, 32L, 29L, 5L, 4L, 3L, 17L, 11L, 95L, 93L, 7L, 38L, 7L, 76L, 76L, 5L, 65L, 5L, 2L, 5L, 5L, 1L, 7L, 12L, 2L, 23L, 6L, 4L, 2L, 3L, 64L, 16L, 5L, 7L, 2L, 5L, 2L, 24L, 413L, 38L, 2L, 38L, 2L, 4L, 47L, 7L, 3L, 41L, 4L, 5L, 1L, 5L, 2L, 7L, 2L, 3L, 5L, 15L, 5L, 1L, 27L, 30L, 2L, 8L, 3L, 4L, 37L, 11L, 1L, 5L, 2L, 35L, 1L, 3L, 3L, 1L, 2L, 2L, 1L, 42L, 1L, 3L, 3L, 3L, 2L, 1L, 4L, 2L, 14L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 12L)), .Names = c("easting", "northing", "ammonia", "nitroxides"), class = "data.frame", row.names = c(NA, -496L)) air.coord <- air[c("easting", "northing")] coordinates(air.coord) <- ~ easting + northing coordinates(air) <- ~ easting + northing vario.emperical <- variogram(nitroxides ~ easting + northing, data = air) plot(vario.emperical) scatterplot3d(air$easting, air$northing, air$nitroxides,
main="3D Scatterplot",
zlab="Nitroxides",
xlab="Easting",
ylab="Northing",
type="h")

symbols(air$easting, air$northing, circles=air$nitroxides, bg = adjustcolor("blue", alpha=0.2)) trend<-lm(log(nitroxides) ~ easting + northing, data=air) summary(trend) plot(trend$residuals,  ylim = c(-4,4))

air3 <-as.data.frame(cbind(easting = air$easting, northing = air$northing,
ammonia = air$ammonia, resid = trend$residuals))
coordinates(air3) <- ~ easting + northing

plot(air3$northing, air3$resid, ylim = c(-4,4))
plot(air3$easting, air3$resid, ylim = c(-4,4))
hist(air3$resid) scatterplot3d(air3$easting, air3$northing, air3$resid,
main="3D Scatterplot Air3",
zlab="Response",
xlab="Easting",
ylab="Northing",
type="h")

### Semi-variograms ###
vario <- variogram(resid ~ easting + northing, data = air3)
vario
plot(vario, cex=1.3, pch=16, col=1, xlab="Distance", ylab="Semivariogram")

az <- seq(0, 165, 15)
vario.direction <- variogram(resid ~ easting + northing,
data = air3, alpha = az,
cutoff = 250, width = 12)
vario.direction
plot(vario.direction)

## Initial Rose Diagram ##
rose(vario.direction, 2.5)

### Rotate Axis ###
coord <- coordinates(air3)
air3.rot <- cbind(rotateaxis(coord, -90), resid = air3$resid) coordinates(air3.rot) <- ~ newx + newy az <- seq(0, 165, 15) vario.direction2 <- variogram(resid ~ newx + newy, data = air3.rot, alpha = az, cutoff = 225, width = 12) rose(vario.direction2, 2.5) xyplot(gamma ~ dist|as.factor(dir.hor), data = vario.direction2, layout=c(4,3), panel = panelgamma0, gamma0=2.1, pch=16, main="Interpolated Directional Ranges",col=1) ## Stretching Rose Diagram ## ratio <- (130/97) coords.new <- coords.aniso(coordinates(air3.rot), aniso.pars=c(0, ratio)) head(coords.new) air4 <- data.frame(anix=coords.new[,1],aniy=coords.new[,2], resid = air3.rot$resid )
coordinates(air4) <- ~anix + aniy

# Rose Diagram After Transformation #
az <- seq(0, 165, 15)
air4.vario<- variogram(resid ~ anix + aniy, data = air4,
width = 12, cutoff = 250, alpha = az)
rose(air4.vario, 2.5)

air4.vario2<- variogram(resid ~ anix + aniy, data = air4,
width = 12, cutoff = 250)

model.vario <- fit.variogram(air4.vario2,
vgm(70000,"Sph",40,20000, anis= c(0,1/ratio)),
fit.method=2)
plot(air4.vario2, model.vario,
cex=1.3,main="Final Isotropic Variogram",
pch=16, ylab = "Semivariogram")

plot(vario.emperical, model.vario,
cex = 1.3, main = "Final Isotropic Variogram",
pch = 16, col = 1, ylab = "Semivariogram")

coords <- coordinates(air)
plot(coords, type="n", xlab="Easting", ylab="Northing",
main= "Plot of Sample and Prediction Sites")
poly <- chull(coords)
poly <- c(poly, poly[1])
lines(coords[poly,])

poly.in <- polygrid(seq(256, 760, 12), seq(64, 700, 12), coords[poly,])
summary(poly.in)

coordinates(poly.in) <- ~ x+y
points(poly.in)
points(coords, pch=16)
krige.out <- krige(air$nitroxides ~ 1, air, poly.in, model= model.vario) northpred <- seq(min(air$northing), max(air$northing)) eastpred <- seq(min(air$easting), max(air$easting)) int.obs <- interp(air$easting,
air$northing, air$nitroxides,
xo=eastpred, yo=northpred)

int.pred <- interp(coordinates(poly.in)[,1],
coordinates(poly.in)[,2],
krige.out$var1.pred, xo=eastpred, yo=northpred) int.se <- interp(coordinates(poly.in)[,1], coordinates(poly.in)[,2], sqrt(krige.out$var1.var),
xo=eastpred, yo=northpred)

zmin <- min(int.obs$z[!is.na(int.obs$z)],
int.pred$z[!is.na(int.pred$z)])
zmax <- max(int.obs$z[!is.na(int.obs$z)],
int.pred$z[!is.na(int.pred$z)])

image(int.obs, ylab="Northing", xlab="Easting",
main = "Observed Concentrations",
zlim=c(zmin, zmax),
col = rev(heat.colors(24)))

image(int.pred, ylab="Northing", xlab="Easting",
main="Kriging Predicted Values",
zlim=c(zmin, zmax), col=rev(heat.colors(24)))

image(int.se, ylab="Northing", xlab="Easting",
main= "Kriging Standard Errors",
col=rev(heat.colors(24)))


Also, the call to krige uses the original data air, on which no detrending has occurred, and it assumes a constant mean (~ 1), whereas the variogram model is based on the residuals in air4 in which any linear trend in $x$ and $y$ trend has been removed.