#HANSEN STEP 2 
#An adjusted urban record is defined only if there are at least three rural
#neighbors for at least two thirds of the period being adjusted. All rural stations within 1000 km are
#used to calculate the adjustment, with a weight that decreases linearly to zero at distance 1000 km.
#The function of the urban adjustment is to allow the local urban measurements to define short-term
#variations of the adjusted temperature while rural neighbors define the long-term change. The break
#in the adjustment line at 1950 allows some time dependence in the rate of growth of the urban
#influence


#This classification was used to identify which stations would be corrected for possible urban warming (adjustments were made only in “urban” areas,
# i.e., those with a population over 50,000) ... the maximum distance used for rural neighbors is 500 km provided that
#sufficient stations are available, and “small-town” (population 10,000 to 50,000) stations are also adjusted.

#The urban adjustment of Hansen et al. [1999] consisted of a two-legged linear adjustment such that 
#the linear trend of temperature before and after 1950 was the same as the mean trend of rural neighboring 
#stations. In the new GISS analysis the hinge year is a variable chosen to be that which allows the adjusted 
#urban record to fit the mean of its neighbors most precisely... The hinge date is now also chosen to minimize the difference between the adjusted urban record and the mean of its
#neighbors.

#in the global analysis we find that the homogeneity adjustment changes the urban record to a
#cooler trend in only 58% of the cases, while it yields a warmer trend in the other 42% of the urban stations.


##FUNCTIONS
source("http://www.climateaudit.org/scripts/station/hansen/step1.txt")
		#loads hansenq
source("http://www.climateaudit.org/scripts/station/hansen/step2/functions.txt")


##LOAD INPUT DATA
 #1. Benchmark Station Data 
	loc.giss="http://www.climateaudit.org/data/giss"
	download.file( file.path(loc.giss,"giss.dset1.tab"),"temp.dat",mode="wb"); load("temp.dat")
  	download.file( file.path(loc.giss,"giss.dset2.tab"),"temp.dat",mode="wb"); load("temp.dat")
 	download.file( file.path(loc.giss,"giss.dset0.tab"),"temp.dat",mode="wb"); load("temp.dat")
 
 #. Station Infor
	url= file.path(loc.giss,"giss.info.dat") 
	stations=read.table(url,sep="\t",header=TRUE,fill=TRUE,quote="")
	names(stations)[3]="name"
	pi180=pi/180;R= 6372.795 #earth's radius in km

 ##WELLINGTON NZ
	id0="50793436001"
	#rural=ruralstations(id0)
	#write.table(rural,file="d:/climate/scripts/station/giss/gistemp.emulation/step2/rural.dat",sep="\t")
	#W=cbind(c(time(Y)),test$chron);dimnames(W)[[2]]=c("year",dimnames(test$chron)[[2]])
	#write.table(W,file="d:/climate/scripts/station/giss/gistemp.emulation/step2/compare.dat",sep="\t")
	#W=ts.union(test$dset1,test$dset2)
	#W=cbind(c(time(W)),W);  dimnames(W)[[2]]=c("year","dset1","dset2")
	#write.table(W,file="d:/climate/scripts/station/giss/gistemp.emulation/step2/target.dat",sep="\t")

	test=emulate_dset2(id0,method="read")
	names(test) 


##PLOT  dset1 versus reference ADJUSTED VALUES

  plotf=function(test) {
	Y=test$emulation
	nf=layout(array(1:3,dim=c(3,1)),heights=c(1.1,1,1.3))
	ylim0=range(Y[,1],na.rm=T)
	par(mar=c(0,3,2,1))
	plot(c(time(Y)),Y[,1],col=1:1,type="l",xlab="",ylab="",axes=FALSE,xlim=c(1880,2007))	
	lines(c(time(Y)),Y[,3],col="green")
	axis(side=1,labels=F);axis(side=2,las=1);box()
	title(main=test$info$name)
	legend(1880,ylim0[2],fill=c(1,3),legend=c("dset1","rural 'reference'"))

	par(mar=c(0,3,0,1))
	ylim1=range(c(Y[,2]-Y[,1],Y[,4]),na.rm=T)
	plot(c(time(Y)),Y[,2]-Y[,1],xlab="",ylab="",axes=FALSE,type="l",ylim=ylim1)
	lines(c(time(Y))[temp],Y[temp,3]-Y[temp,1],lty=2)
	axis(side=1,labels=FALSE);axis(side=2,las=1);box()
 	abline(h=0,lty=2)
	abline(v=range(time(test$dset2)),lty=3,col="red") #1939 1988
	lines(c(time(Y))[temp&!is.na(Y[,1])],fitted(adj0),col=3)
	legend(1880,ylim1[2],fill=c(1,3),legend=c("dset2 adjustment","emulation"))

	ylim0=range(Y[,2],na.rm=T)
	par(mar=c(3,3,0,1))
	plot(c(time(Y)),Y[,2],col=2,type="l",xlab="",ylab="",axes=FALSE,xlim=c(1880,2007))	
	axis(side=1);axis(side=2,las=1);box()
	lines(c(time(Y)),Y[,5],col="green")
	legend(1880,ylim0[2],fill=c(2,3),legend=c("dset2","emulation"))
         }

plotf(test)