##COLLATION GISS ANOM
   #this script uses giss.dset1 and giss.info to make two new objects:
	#1) a new list giss.anom of annual anomaly histories
	#2) a new info table station_versions of station versions meeting some minimal standards.
	  #more than one version per station is possible  
	#there are 7716 versions in giss.dset1, some of which have no data whatever
	#100 versions have no data whatever
	#1473 have less than 20 years 

  ##INPUT INFORMATION
	download.file( file.path(loc.giss,"giss.dset1.tab"),"temp.dat",mode="wb"); load("temp.dat")
	length(giss.dset1) #7364

	url= file.path(loc.giss,"giss.info.dat") 
	stations=read.table(url,sep="\t",header=TRUE,fill=TRUE,quote="") #watch for quote option
	names(stations)[3]="name"
	dim(stations) #[1] 7364   32
	pi180=pi/180;R= 6372.795 #earth's radius in km

	len1=sapply(giss.dset1,length); M=sum(len1);M #[1] 7716
	tapply(!is.na(len1),len1,sum)
		#   1    2    3 
		#7026  324   14 
		#7026+2*324+3*14 #7716

  ##MAKE A DATA FRAME OF STATIONS VERSIONS TO KEEP TRACK OF NA AND ODDBALL FORMATS
	temp7=array(NA,dim=c(M,5)); temp7=data.frame(temp7)
	names(temp7)=c("station","version","nrow","i","j")
	fred=cumsum(len1)
	for (i in 1:length(giss.dset1)){
	  for (j in 1:len1[i]) { 
	    k=fred[i]-len1[i]+j
	    temp7$station[k]= names(test)[i]
	    temp7$i[k]= i
	    temp7$j[k]= j
	    if (!is.null( try(names(giss.dset1[[i]])[j])) &  !(class(giss.dset1[[i]][[j]])=="logical")   ) {
		   temp7$version[k]= names(giss.dset1[[i]])[j]
	           temp7$nrow[k]= sum(!is.na( giss.dset1[[i]][[j]][,18]))}
		}
		}
	   #this condition deals with some annoying NA components that are in the GISS data set for no purpose	
	sum(temp7$nrow<20,na.rm=TRUE) # 1473
 	sum(is.na(temp7$nrow)) #[1] 100
	sum(temp7$nrow>=20,na.rm=TRUE) #6143
	test=match(temp7$station,stations$id);temp=!is.na(test);sum(temp) #7716
	temp7$name=NA;temp7$lat=NA;temp7$long=NA;temp7$alt=NA;temp7$urban=NA;temp7$lights=NA
	temp7$lat[temp]= stations$lat[test[temp]]
	temp7$long[temp]= stations$long[test[temp]]
	temp7$name[temp]= as.character(stations$name)[test[temp]]
	temp7$urban[temp]= as.character(stations$urban)[test[temp]]
	temp7$lights[temp]= as.character(stations$lights)[test[temp]]
	temp7$alt[temp]= as.character(stations$altitude)[test[temp]]
	station_versions=temp7[	(temp7$nrow>=20)&!is.na(temp7$nrow),]
   	nrow(station_versions) # 6143

   ##MAKE ANOMALY VERSION FOR QUALIFYING SERIES
 	giss.anom=rep(list(NA),nrow(station_versions));names(giss.anom)=station_versions$version
	N=nrow(station_versions)
	for (k in 1:N) {
	  i= station_versions$i[k];j = station_versions$j[k]
	  monthly=ts( c(t(giss.dset1[[i]][[j]][,2:13])) ,
		start=c(giss.dset1[[i]][[j]][1,1],1)  ,freq=12)
	 anom.monthly=anom(monthly,method="hansen")[[1]]
	 anom.qtr=hansenq(anom.monthly,method="ts")
	 giss.anom[[k]]=anom.qtr[,5]
	}
	save(giss.anom,file="d:/climate/data/station/giss/giss.anom.tab")
	save(station_versions,file="d:/climate/data/station/giss/station_versions.tab")
	
	length(giss.anom) #6143
	nrow(station_versions) #6143


	#both files uploaded to http://www.climateaudit.org/data//giss