library( chron)

frac.of.year <-
function( obj){
mon<- cbind( cumsum(c(0, 31,28,31,30,31,30,31,31,30,31,30)),
 cumsum(c( 0,31,29,31,30,31,30,31,31,30,31,30)))
leap<- ifelse( obj[,1]%%4 ==0 , 2,  1)
doy<- mon[cbind(obj[,2],leap)] + obj[,3]
foy<- ifelse( leap==1, doy/365, doy/366)
foy
}



 source("ozmax8.info.q")
 hold<-scan("ozmax8.dat")
 N<- length(ozmax8.info$stat.no)
 M<- length(ozmax8.info$dates)
 ozmax8<- list()
 ozmax8$days<- ozmax8.info$dates
 ozmax8$lon.lat<- ozmax8.info$loc
 dd<- chron(ozmax8.info$dates, origin=c( 1,1,1960))
   ytemp<- years(dd)
   ytemp<- as.numeric( levels( ytemp))[ ytemp]

 ozmax8$years<- ytemp + frac.of.year( cbind(ytemp, months(dd), days(dd)) )
 ozmax8$dates <- dd

 PPB<- matrix(hold, ncol=N, nrow=M)
 dimnames( PPB)<- list( as.character(dd), ozmax8.info$stat.no)
 ozmax8$PPB<- PPB

 save("ozmax8", file="ozmax8.rda")

#################
################# testing

library( fields)
load( "ozmax8.rda")

# !is.na(ozmax8$PPB) -> temp
# t(temp)%*% rep( 1, nrow( temp))-> look
 
ind<- c(125,  149)
 ind2<- years(ozmax8$dates) ==1997 

 stats( ozmax8$PPB[ind2,])-> look

png( "ozmax8.png", width=960, height=480)
set.panel( 1,2)

 quilt.plot( ozmax8$lon.lat, look[2,],axes=FALSE, xlab="", ylab="",
 legend.shrink=.7)
 points( ozmax8$lon.lat[ind,], pch="o", cex=2, col="magenta")
points(  ozmax8$lon.lat, pch=".", col="black", cex=2)

title( "Seasonal mean 1997")
US( add=TRUE, col="grey")

 ind2<- years(ozmax8$dates) ==1997 


 matplot( ozmax8$years[ind2], ozmax8$PPB[ind2,ind], 
   type="l", col=c(1:3), lty=1,
   ylab="PPB max 8 hour average", xlab="days")
title("Two individual stations")

dev.off()