# some air quality data, daily surface ozone measurements for 
# the Midwest:
data(ozone2)
s<-ozone2$lon.lat
y<- ozone2$y[16,] # June 18, 1987 
# quick plot of spatial data with map
bubblePlot( s,y)
US( add=TRUE) # add US map
# fitting a thin plate spline surface (always a good place to 
# start). Here the  default smoothing (aka lambda) found by cross-validation
  fit0<- Tps(s,y)
# fits a GCV thin plate smoothing spline surface to ozone measurements.
# Hey, it does not get any easier than this!
  summary(fit0) #diagnostic summary of the fit 
  set.panel(2,2)
  plot(fit0) # four diagnostic plots of fit and residuals.
# quick plot of predicted surface
  set.panel()
  surface(fit0) # contour/image plot of the fitted surface
# see also predictSurface for more control over the evaluation grid
#
  US( add=TRUE, col="magenta", lwd=2) # US map overlaid
  title("Daily max 8 hour ozone in PPB,  June 18th, 1987")
####
  fit2<- spatialProcess( s,y)
# a "Kriging" model. The covariance defaults to a Matern 
# with smoothness 1.0.
# the nugget, sill and range parameters are found by maximum likelihood
# summary, plot, and surface also work for \code{fit2} !
  surface(fit2) # contour/image plot of the fitted surface
  US( add=TRUE, col="magenta", lwd=2) # US map overlaid
  title("Daily max 8 hour ozone in PPB,  June 18th, 1987")
if (FALSE) {
# And 20 approximate conditional draws of the spatial field on a grid
# with uncertainty in the 120PPB contour 
   look<- simLocal.spatialProcess(fit2, M=20)
for( k in 1:20){
contour( look$x, look$y, look$z[,,k], add=TRUE, level=c(120),
  col="white",  drawlabels=FALSE)
}
}  
Run the code above in your browser using DataLab