# gridded elevations on CO grid.
 data( RMelevation)
 COelevation<- crop.image( RMelevation, x)

# attach Colorado climate data
 attach( "COmonthlyMet.rda")

 CO.years[96] #should be 1990

 yr<-96
 y<- (CO.tmax.MAM[yr,]+ CO.tmin.MAM[yr,])/2 
 x<- CO.loc

 Tps( x,y)-> look

# diagnositc plots
 set.panel( 2,2)
 plot( look)

# summary
 summary( look)

# reset panel 
 set.panel()
 surface( look)
 US( add=TRUE)

# superimpose on elevations
 image.plot( COelevation, col=topo.colors(256))
 contour( predict.surface( look), add=TRUE)


# fit additive model with elevation
#strip out NAs

 ind<- !is.na( y)
 y<- y[ind]
 x<- x[ind,]
 elev<-CO.elev[ind]

hold <- pred.lonlat<- rep( 0 ,length(y))
pred.elev<- rep( 0 ,length(y))

NBACK<- 25

for (  k in 1:NBACK){
########### first back fit
# subtract fitted Tps and estimate elevation by OLS

      ytemp<- y- pred.lonlat
      pred.elev<- sreg( elev,ytemp)$fitted.values

###########  second back fit
# fit Tps to residuals from elev OLS
      ytemp<- y - pred.elev
      pred.lonlat<- Tps( x,ytemp)$fitted.values 

########## informal check of convergence for the Tps part
########## relative to the MAD
     cat(k,  mean( abs( hold- pred.lonlat))/mean( abs(y-median(y))),
                                 fill=TRUE)
#
# save new Tps predictions
      hold<-pred.lonlat
}
# Done with back fitting compute final estimates

 ytemp<- y - pred.elev
 tps.obj<- Tps( x,ytemp)

 ytemp<- y- pred.lonlat
 sreg.obj<-  sreg(elev,ytemp)

   grid.list<- list(  x=COelevation$x, y=COelevation$y)
   make.surface.grid(grid.list)-> xg
#
# predict on grid
   look1<- predict( sreg.obj, c( COelevation$z) ) 
   look2<- predict( tps.obj, xg)
   look3<- as.surface(xg, look1+look2)


# coerce to image type objects
   look1<- as.surface(xg, look1)
   look2<- as.surface(xg, look2)

   set.panel( 2,2)
   coltab<- tim.colors(256)

   quilt.plot( x,y, axes=FALSE,col=coltab)
   contour( COelevation, levels=c(1500, 3000),
                col="grey",labex=0,add=TRUE)
   US( add=TRUE, col="darkgreen")
 

   image.plot( look3, zlim=range(y), col=coltab, axes=FALSE)
   image( look1, zlim=range(y),col=coltab, axes=FALSE)
   image( look2, zlim=range(y),col=coltab, axes=FALSE)