library( fields) load( "ozmax8.rda") # sanity check hold2<- stats( ozmax8$PPB)[2,] quilt.plot( ozmax8$lon.lat, hold2) US( add=TRUE, col="grey50", lwd=2) # years variable includes the fraction of year fracyear<- ozmax8$years - trunc(ozmax8$years) X<- cbind( sin( 2*pi*fracyear), cos( 2*pi*fracyear), sin( 3*pi*fracyear), cos( 3*pi*fracyear), sin( 4*pi*fracyear), cos( 4*pi*fracyear)) X<- scale( X, center=TRUE, scale=FALSE) coef.hold<- matrix( NA, ncol=7, nrow=513) for( k in 1:513){ cat(k, " ") station.data<- ozmax8$PPB[,k] ind<- !is.na( station.data) coef.hold[k,]<-lsfit( X[ind,], station.data[ind])$coef } coef.svd<- svd( coef.hold) # take a look at the leading mean coefficient # largely the spatial mean: quilt.plot( ozmax8$lon.lat,coef.hold[,1] ) US(add=TRUE, col="grey50", lwd=2) # retain 3 PCs coef.smooth<- coef.svd$u%*% diag( c( coef.svd$d[1:3], rep(0,4)))%*%t(coef.svd$v) # subtracting seasonality big.season<- cbind( rep( 1,920), X) %*% t( coef.smooth) # mean of seasonal efects being subtracted hold3<- colMeans(big.season) quilt.plot( ozmax8$lon.lat,hold3 ) US(add=TRUE, col="grey50", lwd=2) ozmax8.residual<- ozmax8$PPB - big.season # simple check that seasonality has been removed set.panel( 2,1) bplot.xy( rep(fracyear, 513), c( ozmax8$PPB)) bplot.xy( rep(fracyear, 513), c( ozmax8.residual))