library(fields)
# set up the data set to be easier to type
  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]
# finds some variograms
  vg0<-vgram( xHW1, yHW1, N=15)
  res1<-  lm( yHW1~xHW1)$residuals
  res2<-  lm( yHW1~xHW1 +elevHW1)$residuals
  vg1<- vgram( xHW1, res1, N=15)
  vg2<- vgram( xHW1, res2, N=15)
  fit1<- MLE.Matern( xHW1,yHW1, smoothness=1.0, m=1,Dmax=6 )
  fit2<- MLE.Matern( xHW1,yHW1, smoothness=1.0, m=2, Dmax=6)
  fit3<- MLE.Matern( xHW1,yHW1, Z=elevHW1, smoothness=1.0, m=2, Dmax=6)
#plot'em
  fields.style()
  plot( vg0$centers, vg0$stats[2,], ylim=c(0,20), ylab="Variogram", xlab="degrees", col=1)
  points(vg1$centers, vg1$stats[2,], col=2 )
  points(vg2$centers, vg2$stats[2,], col=3 )
  lines( fit1$vgram, lwd=3, col=1)
  lines( fit2$vgram, lwd=3, col=2)
  lines( fit3$vgram, lwd=3, col=3)
#dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture4/pix/vario1.pdf")

# simulate a random surface and add a spatial trend.
   xg<- make.surface.grid( list(x= seq( 0,1,,50), y = seq( 0,1,,50)) )
   K<- exp( -rdist(xg,xg))
   set.seed(123)
   g<- t(chol(K)) %*%rnorm(nrow(K))
   spatial.drift<- 10 + 2*xg[,1] - xg[,2]
# plot these fields
   fields.style()
   set.panel(1,3)
   par( mar=c(1,1,2,6))
   image.plot( as.surface( xg, spatial.drift), axes=FALSE, xlab="", ylab="")
   title("spatial trend")
   image.plot( as.surface( xg, g) , axes=FALSE, xlab="", ylab=""            )
   title("random process")
   image.plot( as.surface( xg, spatial.drift+ g), axes=FALSE, xlab="", ylab="")
   title( "spatial trend + random process")
#dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture4/pix/drift1.pdf")

# plot variogram of trend and random components
  dev.off()
  set.panel()
  fields.style()
  vg1<- vgram( xg, spatial.drift, N=15)
  vg2<- vgram( xg,  g, N=15)
  vg3<- vgram( xg,  spatial.drift + g, N=15)
  matplot( vg1$centers, cbind( vg1$stats[2,], vg2$stats[2,], vg3$stats[2,]), col=1:3, type="p",
  pch=16, xlab="distance", ylab="Variogram") 
#dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture4/pix/drift2.pdf")

# not talked about in lecture 
# a variogram of a fixed but complex function 
 g<-  sin(xg[,1]*10)*sin(xg[,2]*10)
 g <- g* exp( -( (xg[,1]-.5)^2 + (xg[,1]-.5)^2)/.3)
 look<- vgram( xg, g, N=30)
 plot( look$centers, look$stats[2,], xlab="Distance", ylab="variogram")
#dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture4/pix/counter1.pdf")
 fields.style()
 image.plot( as.surface(xg, g))
#dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture4/pix/counter2.pdf")

# different spatial predictions for the climate data
# m=1 , m=2 and  m=2 with elevation
obj0<- mKrig( xHW1, yHW1, Covariance="Matern", theta=1.1,
                        smoothness=1.0,  lambda= (1.26^2)/10.1, m=1 )
obj1<- mKrig( xHW1, yHW1, Covariance="Matern", theta=1.13,
                        smoothness=1.0,  lambda= (1.26^2)/10.5, m=2 )
obj2<- mKrig( xHW1, yHW1, Z= elevHW1, Covariance="Matern", theta=1.13,
                        smoothness=1.0,  lambda= (1.10^2)/1.09, m=2 )
# evaluate on a grid (default grid id 80X80)
out0<- predict.surface( obj0)
out1<- predict.surface( obj1)
out2<- predict.surface( obj2, drop.Z=TRUE)
# plot m=1 and m=2 cases
set.panel( 1,2)
fields.style()
par( mar=c(1,1,2,2))
image.plot( out0, axes=FALSE, xlab="", ylab="")
US( add=TRUE)
title( "constant (m=1) ")
image.plot( out1, axes=FALSE, xlab="", ylab="")
US( add=TRUE)
title("linear trend (m=1)")
#dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture4/pix/mfits1.pdf")

# smooth surface for model with elevation.
out2<- predict.surface( obj2, drop.Z=TRUE, nx=200, ny=200)
grid.list<- list( x= out2$x, y=out2$y)
# approximate elevations on this grid (this is cheating a bit!)
data(RMelevation)
# take the Rocky Mountain data set (at 4km resolutoin) and evaluate on
# the grid assumed by predict.surface using bilinear interpolation.
elev.grid<- interp.surface.grid( RMelevation, grid.list)
# highlight the d[4] estimated coefficient
elev.coef<- obj2$d[4] # where the coefficients are kept!

# add two surfaces together
# note we are adding the matrices of values but _not_ adding the lists out2 and elev.grid
#
full.surface<-  out2$z + elev.grid$z*elev.coef
set.panel( 1,3)
fields.style()
par( mar=c(1,1,1,3))
image.plot( grid.list$x, grid.list$y,out2$z, axes=FALSE, xlab="", ylab="")
US( add=TRUE)
title("linear (m=2) + g(x)")
image.plot( grid.list$x, grid.list$y, elev.grid$z*elev.coef, axes=FALSE, xlab="", ylab="")
US( add=TRUE)
title("elevation only")
image.plot( grid.list$x, grid.list$y, full.surface, axes=FALSE, xlab="", ylab="")
US( add=TRUE)
title(" linear + elevation + g(x)")
#dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture4/pix/mfits2.pdf")