# 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