library("zoo")
library("chron")

source("gsod.common.R")

# mostly just so I can test on small data sets
fdir <- "../../gsod/data/ftproot"

# get a list of stations
ishwids=c(6,1,5,1,28,1,8,1,4,3,6,1,7,1,6,4,8,1,8)
ishcols=c("USAF","a","WBAN","b","STATION NAME","c","CCST","f","CALL","g","LAT","h","LON","i","ELEV","j","BEG","k","END")
ish <- read.fwf("../../gsod/data/ftproot/ish-history.txt-rb",widths=ishwids,col.names=ishcols)

# get just northern hemisphere (30N-70N)
ish <- ish[as.numeric(ish[,11])/1000>=lat0 & as.numeric(ish[,11]/1000)<=lat1,]
ish <- ish[as.numeric(ish[,13])/1000>=lon1 & as.numeric(ish[,13]/1000)<=lon0,]
ish <- ish[!is.na(ish[,11]),]
ish <- ish[(ish[,11]!=-99999),]
ish <- ish[!is.na(ish[,13]),]
ish <- ish[(ish[,13]!=-999999),]
ish <- ish[!is.na(ish[,15]),]
ish <- ish[(ish[,15]!=-99999),]
ish <- ish[!is.na(ish[,1]),]
ish <- ish[!is.na(ish[,3]),]

stncnt <- nrow(ish)
stncnt

# column names for data file 
cn <- c("USAF","WBAN","YEARMODA","TEMP","A",
        "DEWP","B","SLP","C","STP","D",
        "VISIB","E","WDSP","F","MXSPD",
        "GUST","MAX","MIN","PRCP","SNDP","FRSHTT")

# clear table
tabcol <- cbind("stnid","yearmoda","temp")
tabcol2 <- cbind("stnid","lat","lon","elev","yearmoda","temp","snow")
#write.table(tabcol2,stn_da2.txt,sep=" ",row.names=F,col.names=F,append=F,quote=F)

# used to pad the year data 
for (c in 1:stncnt) { 
	elev <- as.numeric(ish[c,15])
	lat <- as.numeric(ish[c,11])
	lon <- as.numeric(ish[c,13])
	stntab <- NULL
	stntab2 <- NULL
	name <- paste(sprintf("%06d", ish[c,1]),"-",sprintf("%05d", ish[1,3]),sep="")
	newfile <- paste("stndat/",name,".daily.txt",sep="")
	ycnt <- 0
	for (year in y0:y1) {
		stnfile <- paste(name,"-",year,".op",sep="")
    		stnfile <- paste(fdir,"/",year,"/",stnfile,sep="")
		if (file.exists(stnfile)) {
print(stnfile)
        		gsod <- read.table(stnfile,header=F,skip=1,col.names=cn)
			m1 <- nrow(gsod[gsod[,3]<(year*10000)+0132,]) # jan
			m2 <- nrow(gsod[gsod[,3]>(year*10000)+0132 & gsod[,3]<(year*10000)+0232,]) # feb
			m3 <- nrow(gsod[gsod[,3]>(year*10000)+0232 & gsod[,3]<(year*10000)+0332,]) # mar
			m4 <- nrow(gsod[gsod[,3]>(year*10000)+1032 & gsod[,3]<(year*10000)+1132,]) # nov
			m5 <- nrow(gsod[gsod[,3]>(year*10000)+1132,]) # dec
			# if enough data in each of			4 quarters
			if (m1>tau & m2>tau & m3>tau & m4>tau & m5>tau) { 
				# convert to C and adjust for height
				t <- as.numeric(gsod[,4])
				s <- as.numeric(gsod[,21])
				t <- as.integer((5/9) * (t - 32.0) * 10.0)
				stntab2 <- rbind(stntab2,cbind(name,lat,lon,elev,gsod[,3],t,s))
				ycnt <- ycnt + 1
			}
		}
	}

	if (!is.null(stntab2)){
	  # write individual dat file
	  colnames(stntab2) <- tabcol2
	  write.table(stntab2,newfile,sep="\t",row.names=F,col.names=T,append=F,quote=F)

	  }
}