##CALCULATION OF URBAN-RURAL DIFFERENCE IN PETERSON 2003 NETWORK

    	peter=read.table("http://www.climateaudit.org/data/station/peterson/station.dat",sep="\t",header=TRUE,fill=TRUE)
	source("http://www.climateaudit.org/scripts/station/read.ghcnd.txt"
	source("http://www.climateaudit.org/scripts/gridcell/collation.functions.txt"

##STEPS
   #collate GHCND data for PEterson network (282 of 289 series)

   ##COLLATE PETERSON NETWORK
	id=peter$ghcnd
	index=(1:length(id))[!is.na(id)]
	peterson=rep(list(NA),289)
	names(peterson)=id

	for (i in index)){
	peterson[[i]]=read.ghcnd(usid=id[i])
	}
	save(peterson,file="d:/climate/data/station/peterson/peterson.collation.tab")
		#this yielded 282 of 289 series

	test=sapply(peterson,dim);test=sapply(test,length); temp=(test==0)  #7 NA
	index=(1:289)[!temp]

    ##ORGANIZE IN TIME SERIES
	chron=NULL
	for( i in index){
	  test=ts(peterson[[i]]$tmean,start=c( floor(peterson[[i]]$time[1]/100),1),freq=12)	
          chron=ts.union(chron,test)
	  }
	dimnames(chron)[[2]]=id[index]

    ##MAKE 1961-1990 ANOMALIES
	anomaly=chron
	for(i in 1:ncol(chron)){anomaly[,i]=anom(chron[,i])}

    ##MAKE AVERAGES FOR "URBAN" and "RURAL" NETWORKS
	f=function(x){f=filter.combine.pad(x,truncated.gauss.weights(36))[,2];f}
 
     #all series 
   	average0=ts(apply(anomaly,1,mean,na.rm=TRUE),start=c(tsp(chron)[1],1),freq=12)
	ts.plot(average0,col="grey80")
	lines(c(time(average0)),f(average0),lwd=2)

     #urban and rural     
	temp=(peter$urban[index]=="u")
	average1=ts(apply(anomaly[,temp],1,mean,na.rm=TRUE),start=c(tsp(chron)[1],1),freq=12)
	temp=(peter$urban[index]=="r")
	average2=ts(apply(anomaly[,temp],1,mean,na.rm=TRUE),start=c(tsp(chron)[1],1),freq=12)

	
	nf=layout(array(1:3,dim=c(3,1)),heights=c(1.1,1,1.3));ylim0=c(-1.6,1.6)
	par(mar=c(0,4,2,1))
	plot(c(time(average1)),average1,col="grey80",type="l",ylim=ylim0,xlab="",ylab="dge C",axes=FALSE)
	axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2)
	lines(c(time(average1)),f(average1),lwd=2)
	text(1880,1.6,"'Urban'",pos=4,font=2)

	par(mar=c(0,4,0,1))
	plot(c(time(average2)),average2,col="grey80",type="l",ylim=ylim0,xlab="",ylab="dge C",axes=FALSE)
	axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2)
	lines(c(time(average2)),f(average2),lwd=2)
	text(1880,1.6,"'Rural'",pos=4,font=2)

	par(mar=c(3,4,0,1))
	plot(c(time(average2)),average1-average2,col="grey80",type="l",ylim=ylim0,xlab="",ylab="dge C",axes=FALSE)
	axis(side=1,labels=FALSE);axis(side=2,las=1);box();abline(h=0,lty=2)
	lines(c(time(average2)),f(average1-average2),lwd=2)
	text(1880,1.6,"Difference",pos=4,font=2)
	
	fm=lm(I(average1-average2)~c(time(average1)) )
	coef(fm)
#      (Intercept) c(time(average1)) 
#    -11.797121267       0.005986827 


	temp=!is.na(average1-average2)
	lines(c(time(average1))[temp],fm$fitted.values)