##COMPARE NASA AND NOAA SERIES
	anom= function(x) {
                K=floor(tsp(x)[1])
                if( is.null(dim(x))) {test=t(array(x,dim=c(12,length(x)/12)) )
                m0=apply(test[(1961:1990)-K+1,],2,mean,na.rm=TRUE)
                test=scale(test,center=m0,scale=FALSE)
                anom=ts(c(t(test)),start=c(K,1),freq=12)
                anom} else {anom=x;
                for(i in 1:ncol(x)) {test=t(array(x[,i],dim=c(12,nrow(x)/12)) )
                m0=apply(test[(1961:1990)-K+1,],2,mean,na.rm=TRUE)
                test=scale(test,center=m0,scale=FALSE)
                anom[,i]=c(t(test))}}
                anom
                }

  ##COLLATE NOAA VERSION
      #data comes in an anooying format
	url="http://www1.ncdc.noaa.gov/pub/data/cirs/drd964x.tmpst.txt"
	fred=readLines(url)
	state=as.numeric(substr(fred,1,3))
	temp=(state==110);sum(temp) #113
		#[6553] "1100022006  39.50  35.40  43.50  56.10  63.10  71.60  77.20  74.20  63.90  53.10  44.60  36.80   "
		#[6554] "1100022007  31.50  32.90  48.20  51.30  63.20  70.60  75.50  75.40 -99.90 -99.90 -99.90 -99.90   "
	x=fred[temp];
	widths0=cumsum(c(3,4,4,rep(7,12)));M=length(widths0)
	widths0=cbind( c(0,widths0[1:(M-1)]),widths0) 
	X=array(NA,dim=c(sum(temp),M) );
	X=data.frame(X)
	f=function(i){ f=substr(x,widths0[i,1],widths0[i,2]);f}

	for(i in 1:M){X[,i]=f(i)}
	for(i in 1:M) {X[,i]=as.numeric(X[,i])}
	temp=(X== -99.9);X[temp]=NA

	#apply(X[112:113,4:11],1,mean)
		#57.575 56.075 
	X[113,12:15]=X[112,12:15]
	x= ts( c(t(X[,4:15])),start=c(X[1,3],1),freq=12)
	lines(c(time(x)),filter(x,rep(1/24,24)),lwd=2)
	x=round( (x-32)*5/9,2)	
	noaa=anom(x)

	tsp(noaa) # [1] 1895.000 2007.917   12.000
		#ts.plot(noaa,col="grey80")
		#y=ts.annavg(noaa)
		#lines(c(time(y)),y,lwd=2)

##GISS COMPARISON
	giss=read.table("http://www.climateaudit.org/data/giss/giss_us.dat")
	giss=ts(giss[,2:3],start=giss[1,1])
		#giss[1,]
		#Jul07 Aug07 


#COMPARE NASA AND NOAA #OCt review
	Y=ts.union(ts.annavg(noaa),giss[,"Aug07"])
	Y=cbind(Y,Y[,1]-Y[,2])
	temp=(time(Y)>=1890);Y=ts(Y[temp,],start=1890)
	par(mar=c(3,4,2,1))
	plot(c(time(Y)),Y[,3],type="l",ylab="deg C");
	title(main="NOAA minus NASA : Aug 2007 Versions")
	abline(h=0,lty=2)
	index=c(time(Y))
	temp=(time(Y)>=1940)&!is.na(Y[,3])
	fm=lm(Y[temp,3]~index[temp]);coef(fm)
	lines(index[temp],fm$fitted.value,col="red")
	text(locator(1),paste("Trend: ",round(100*coef(fm)[2],2),  " deg/century"),col="red",font=2, pos=4)