#Read files - see www.rankexploits.com/musings for the ant_recon_names.txt file
awsnames=read.table("ant_recon_names.txt", stringsAsFactors=FALSE)
full_recon=read.table("ant_recon.txt")
full_recon=ts(data=full_recon, start=1957, end=(2006+11/12), deltat=1/12, ts.eps=.03)
aws_recon=read.table("ant_recon_aws.txt")
aws_recon=ts(data=aws_recon, start=1957, end=(2006+11/12), deltat=1/12, ts.eps=.03, names=awsnames[[1]])
temp=read.table("awsreconcorrected.txt")
temp=ts(data=temp, start=1957, end=(2006+11/12), deltat=1/12, ts.eps=.03, names=awsnames[[1]])
aws_recon=cbind(aws_recon, temp)	#Store Steig's corrected recon as columns 64-126

i=1;while(i<64){aws_recon=cbind(aws_recon, full_recon[,awsnames[i,2]]);i=i+1}

#Entering arguments for the linear fit functions.
#*num.trends*:  8 total trends.
#*trend.years*:  row 1 are start years, row 2 are end years (not inclusive)
num.trends=matrix(1:8, nrow=1);trend.years=matrix(c(1957, 1980, 1957, 2007, 1969, 2000, 1979, 2004, 1982, 1994.5, 1994.5, 2007, 1982, 1999, 1982, 2007), nrow=2)

#Define function to create logical array for the trend times - trend times in columns
map=function(x) {(time(aws_recon)>=I(x[1]) & time(aws_recon)<I(x[2]))}
temp=apply(trend.years, 2, map)	#Run it.

#Define generic function to generate linear fits.
#1:length(x)/12 used to return slope in degC/year,
#since x is in monthly buckets.
trendf=function(x) {fm=lm(x~I((1:length(x))/12));trendf=unlist(fm$coef[2]);trendf}

#Define function to repeat the linear fit for each of the 8
#trend periods.  Calls the linear fit function for each column of the
#logical array *temp*
repeatf=function(y) {cbind(apply(aws_recon[temp[,y],],2,trendf))}
aws_trends=10*(apply(num.trends, 2, repeatf)) #factor of 10 to convert to C/decade

#Make names that make sense
atnames=as.vector(paste(trend.years[1,], "-",(trend.years[2,]), sep=""));rownames(aws_trends)=c(paste("ORIG_", awsnames[[1]], sep=""),paste("CORR_", awsnames[[1]], sep=""),paste("TIR_", awsnames[[1]], sep=""))
colnames(aws_trends)=atnames
atrnames=c("Original Reconstruction", "Corrected Reconstruction", "AVHRR Reconstruction")

#Define function to get raw READER files, convert to anomalies, convert to time series
num.files=matrix(1:63, ncol=1);mode.vector=matrix()	#total number of files to load; matrix to store row modes for anomaly calcs
getREADER=function(x) {

	#Use the one below if you saved the READER files to your HD (remove the # in front of it)
	#filename=paste(awsnames[x,1], ".All.temperature.txt", sep="")	#Make filenames

	#Use the one below to get the READER files on line (remove the # in front of it)
	filename=paste("http://www.antarctica.ac.uk/met/READER/aws/", awsnames[x,1], ".All.temperature.txt", sep="")	#Make filenames

	temp=read.table(filename, sep="", skip=1, na.strings="-", stringsAsFactors=FALSE)	#Read 'em!
	datemap=(temp[,1]>1956 & temp[,1]<2007); temp=temp[datemap,]	#Make a logical array to throw away pre-1957 and post-2006
	starty=temp[1,1];	temp=temp[,-1]	#Save the starting year and get rid of the year column to prepare for matrixing
	temp=t(temp); temp=matrix(temp, ncol=1)	#Transpose the matrix (so the rows become columns for stacking) and set columns=1

	#Make anomalies.  Because Steig did not use all available months for the anomalies sometimes,
	#making anomalies using monthly means of the actual data results in offsets for each month.  To
	#prevent, anomalies will be formed by finding the mode of the RECON - ACTUAL
	#for each month.  This ensures that ACTUAL - RECON for all ACTUAL data used equals
	#zero.  Any points different than zero indicate ACTUAL data NOT USED for the
	#reconstruction.
	temp2=aws_recon[,x]-ts(temp, start=starty, end=2006+11/12, deltat=1/12, ts.eps=0.03);temp2=matrix(temp2, nrow=12)
	i=1;while(i<13){
	#sort each row for the mode and store in mode.vector
	mode.vector[i]=ifelse(length(temp2[i,][!is.na(temp2[i,])])>0, names(sort(-table(temp2[i,1:(length(temp2)/12)])))[1], 0)
	i=i+1}
	temp=matrix(temp, nrow=12);temp=temp+as.numeric(mode.vector);temp=matrix(temp, ncol=1)	#Make anomalies

	#Convert to time series and return.
	temp=ts(temp, start=starty, end=2006+11/12, deltat=1/12, ts.eps=0.03)
	temp} #Return our time series of anomalies
raw_data=apply(num.files, 1, getREADER)

i=1;while(i<64){aws_recon=cbind(aws_recon, raw_data[[i]]);i=i+1}
colnames(aws_recon)=c(paste("ORIG_", awsnames[[1]], sep=""), paste("CORR_", awsnames[[1]], sep=""), paste("TIR_", awsnames[[1]], sep=""), paste("ACT_", awsnames[[1]], sep=""))	#Get rid of cbind's autonames

#Cleanup
rm(full_recon, getREADER, i, map, num.files, num.trends, repeatf, trend.years, trendf)

#x defines the station 1-63 in alphabetical order, y gives the trend period to look at
x=1;y=1
barplot(aws_trends[(y-1)*63+(1:63), x], ylab="deg C/decade", ylim=c(-.4, 1.0), xlab="AWS Stations in Alphabetical Order (Right to Left)", main=atnames[x], sub=atrnames[y], axisnames=FALSE, col=rainbow(63))	

#Same as above, but only using the stations from S1
stns=c(2,5,6,7,9,10,17,19,21,22,25,30,31,33,35,39,40,41,44,45,47,51,55,59,61,63) #Table S1 stations
barplot(aws_trends[(y-1)*63+stn, x], ylab="deg C/decade", ylim=c(-.4, 1.0), xlab="AWS Stations (Table S1 Only) in Alphabetical Order (Right to Left)", main=atnames[x], sub=atrnames[y], axisnames=FALSE, col=rainbow(26))	

#Scatterplot
plot(aws_recon[1:600, 1:63], aws_recon[1:600,127:189], type="p", xlab="AVHRR", ylab="AWS", xlim=c(-10,10), ylim=c(-10,10))

#xlimit gives the period to graph, l defines whether you use the original, corrected, or AVHRR recons, z gives the stations 1-63 in alphabetical order
xlimit=c(1982, 2007);l=1;z=1
temp=(aws_recon[,((l-1)*63)+z]-aws_recon[,189+z]);lmf=lm(temp~I(time(temp)));plot.ts(temp, xlim=xlimit, type="p", xlab=paste(atrnames[l], " - READER (Actual) Data", sep=""), main=awsnames[z,1], sub=paste("Slope = ", round(as.numeric(lmf$coefficients[2])*10, 4), " deg C/Decade", sep=""));abline(lmf)