#GISS STATION data download
   ##SAVED AS "http://www.climateaudit.org/scripts/ghcn.txt" 	
   ###GIS station data id's can be obtained from http://data.giss.nasa.gov/gistemp/station_data/

   ###the form of GIS station data is as follows for Niwot Ridge CO
	#http://data.giss.nasa.gov/work/gistemp/STATIONS//tmp.425724690080.2.1/station.txt
	#"http://data.giss.nasa.gov/work/gistemp/STATIONS//tmp.425724690080.0.1/station.txt"
	#425724690060.0.1 #raw GHCN + UHCN correction
	#425724690080.1.1 #combining sources
	#425724690080.2.1 after homogeneity adjustments 
	#for Niwot Ridge, all 3 series are identical
	##other types #combined; adjusted

##download station data from GISS as time series
	#freq - "annual" or "monthly"
	#type "UHCN" or "combined" or "adjusted"
  ghcn<-function(id,type="UHCN",freq="annual") {
	if (type=="UHCN") K<-0; if(type=="combined") K<-1; if (type=="adjusted") K<-2
	url<- paste( "http://data.giss.nasa.gov/work/gistemp/STATIONS//tmp",id,K,"1/station.txt",sep=".")
	test4<-scan(url,skip=1)
	N<-length(test4) ##N should be divisible by 18  
	test4<-array(test4,dim=c(18,N/18))
	test4<-t(test4)
	start0<-test4[1]
	temp<- (test4 == 999.9)
	test4[temp]<-NA
	if (freq=="monthly") { ghcn<- c(t(test4[,2:13])); ghcn<-ts(ghcn,start=c(start0,1),freq=12)}
	if (freq=="annual") {ghcn<-  ts( test4[,14],start=start0)}
	ghcn
	}

#convert monthly time series to an anomaly series centered on average of series
   #index - default is centered on period of series; otherwise input start and end as c(1960,1990) 
   anomaly<-function(x,index=0) {
	N<-length(x)  #presume div by 12
	y<-array(x,dim=c(12,N/12))
	y<-t(y)
	if (length(index)== 1) y<-scale(y,scale=FALSE) else {
		temp<- ((index[1]-tsp(x)[1]+1):(index[2]-tsp(x)[2]+1)); 
		m0<-apply(y[temp,],2,mean,na.rm=TRUE);
		y<-scale(y,scale=FALSE,center=m0)}
	anomaly<-c(t(y))
	anomaly<-ts(anomaly,start=c(tsp(x)[1],1),freq=12)
	anomaly
	}