library( fields)



# Back to the colorado data.
  data(COmonthlyMet)
  good<- !is.na(CO.tmin.MAM.climate)
  xHW1<- as.matrix(CO.loc[good,])
  yHW1<- CO.tmin.MAM.climate[good]
  elevHW1<- CO.elev[good]
## estimated value of theta from previous fit
# see lecture 7 
 theta2<- 0.5343372
 sm2 <-  1.25
 
  fit2<- Krig( xHW1, yHW1, Z= elevHW1,Covariance="Matern",
              smoothness=sm2, theta=theta2, method="REML")
  sur2<- predict.surface( fit2, drop.Z=TRUE, nx=50, ny=50)
  surcheck<- predict.surface( fit2, drop.Z=TRUE, nx=150, ny=150)

  xg<- make.surface.grid( list( x= sur2$x, y= sur2$y))
  ghat<- predict( fit2, drop.Z=TRUE, x=xg)
  
# this takes about 1.5 minutes
  system.time(bigV<- predict.se( fit2, xg, drop.Z=TRUE, cov=TRUE))
  V.chol<- chol( bigV)
  N<- nrow( xg)
  set.seed( 233)
  gstar<- matrix( ghat, ncol=25, nrow=N) +
              t( V.chol)%*%matrix( rnorm(N*25), ncol=25, nrow=N)
  gstar.se<- sqrt(diag( bigV))

#  image.plot( as.surface( xg, ghat), axes=FALSE, xlab="", ylab="")
  set.panel(3,3)
  fields.style()
  zr<- range( gstar)
  ctab<- designer.colors(256,col=c("brown", "wheat2","seagreen", "blue4"))
  par( oma=c(0,0,0,4))
  par( mar=c(1,1,0,0))
  image( as.surface( xg, ghat), axes=FALSE, xlab="", ylab="", zlim=zr,
             col=ctab)
  contour( as.surface( xg, ghat), levels=c(9.5), col="magenta", lwd=2,
               drawlabels=FALSE, add=TRUE)
  box( lwd=6, col="grey30")
  US( add=TRUE)    
  for ( k in 1:8){
      image( as.surface( xg, gstar[,k]), axes=FALSE, xlab="", ylab="", zlim=zr,
           col= ctab)
      contour( as.surface( xg, gstar[,k]), levels=c(9.5), col="magenta", lwd=2,
               drawlabels=FALSE, add=TRUE)
     US( add=TRUE, lty=3)
    }
    par( oma=c( 0,0,0,2))  
    image.plot( legend.only=TRUE, zlim=zr, col=ctab)
##  dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture8/pix/figCOtemp.pdf")

## CV residuals

    SE<- predict.se( fit2, Z= fit2$Z)
    Aweights<- diag( Krig.Amatrix(fit2, Z= fit2$Z))
    CVresiduals<- fit2$residuals/ ( 1- Aweights)
    CV.std<- CVresiduals/ sqrt( fit2$shat.MLE^2 + SE^2)

    fields.style()
    qqnorm(CV.std, cex=1.5)
    abline( 0,1, col="green")
    
##  dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture8/pix/figCOqqplot.pdf")
    
    


#   points( x,y)

### bivariate normal

Sigma<- rbind(  c(1,.9),
                c(.9,1))
Sigma.chol<- chol( Sigma)
set.seed( 233)
U<- matrix( rnorm( 50000*2), ncol=50000, nrow=2)
Z<- t(Sigma.chol)%*%U
Z<- t( Z)
#############
#  elliptical contour
R<- 2.32
theta<- seq( 0, 2*pi,,400)
u<- R*cbind(sin( theta), cos( theta))
x<- t(Sigma.chol)%*%t(u) 
xcontour<- t( x)
print(100*sum( in.poly(Z,xcontour))/50000)
# rectangular box
xrect<- rbind(   c(-2,-2),
                  c(-2,2),
                  c( 2,2),
                  c( 2,-2),
                  c(-2,-2))
print(100*sum( in.poly(Z,xrect))/50000)
######
fields.style()
plot( Z, xlab= "z1", ylab="z2", pch=15, cex=.15, col="green4",
     xlim=c(-4,4), ylim=c(-4,4))
  xline( c( -2,2), col=1, lwd=.5)
  yline( c( -2,2), col=1, lwd=.5)
  lines( xrect, lwd=2)
  lines( xcontour,lwd=2)
## dev.copy2pdf( file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture8/pix/fig1a.pdf")