##CRN GRIDDING

	source("http://www.climateaudit.org/scripts/station/collation.functions.txt")
	library(akima);library(fields)

##MASK FUNCTIONS on 40 x 40 
   gulfmask=function(Z){
		Z[21:31,1:2]=NA;Z[21:30,3]=NA;Z[22:30,4]=NA;
		Z[23:30,5]=NA;Z[c(23:24,27:30),6]=NA; Z[27:28,7]=NA; Z}	

   atlmask=function(Z){
		Z[40,c(1:26,29:31)]=NA;  Z[38:39,1:26]=NA; Z[37,1:23]=NA;
		Z[36,1:16]=NA;Z[35,1:15]=NA; Z[34,1:13]=NA;Z[33,1:11]=NA;
		 Z}	

   canadamask=function(Z){
		Z[40,37:40]=NA;  Z[37:39,35:40]=NA; Z[36,34:40]=NA;
		Z[34:35,32:40]=NA;Z[32:33,29:40]=NA; Z[31,c(28:29,33:40)]=NA;
		Z[30,36:40]=NA;Z}

   ushcn_wt=function(index_good,ushcn=ushcn_tobs) {
  	good=ushcn[,index_good] #2544 53
	M=ncol(good)  #anomalies
	good.anom=good
	for(i in 1:M) good.anom[,i]=anom(good[,i]) #calculate anomaly
	good.ann=ts( array(NA,dim=c(nrow(good)/12,M)),start=tsp(good)[1]) #make array
 	for(i in 1:M){ good.ann[,i]=ts.annavg(good.anom[,i])}  #calculate annual averages
	temp=(c(time(good.ann))>=1890); good.ann=ts(good.ann[temp,],start=1890) #restrict to 1890 and later
	time0=c(time(good.ann)); N=length(time0) #1890:2006
	for(i in 1:M) good.ann[N,i]=mean(good.anom[(nrow(good.anom)-11):nrow(good.anom),i],na.rm=T)  #extrapolate for final year from available
     
	annual=rep(NA,N)
	for (j in 1:N){
	   X=data.frame(details[index_good,c("long","lat")],good.ann[j,]) #locations of "good" stations
	   names(X)=c("long","lat","anom")
	   temp=!is.na(X$anom)
	   gridded <- interp(X$long[temp], X$lat[temp], X$anom[temp])
 	   gridded$z=gulfmask(gridded$z)
	   gridded$z=atlmask(gridded$z)
	   gridded$z=canadamask(gridded$z)
	   annual[j]=weighted.mean(gridded$z,rep(cos(gridded$y*pi/180),40),na.rm=T)
	 }
	ushcn_wt=ts(annual,start=1890)
	ushcn_wt
	}


##WATTS RATINGS
  	url="http://www.climateaudit.org/data/ushcn/details.dat"  #USHCN station details
	details=read.table(url,sep="\t",header=TRUE,fill=TRUE)	

	loc="http://www.climateaudit.org/data/station/watts/stationlist.txt"#
		## "d:/climate/data/station/watts"
		#downloaded from http://www.surfacestations.org/downloads/USHCN_stationlist.xls 
		#some columns deleted
		## replace ' with _	
	watts=read.table(loc,sep="\t",fill=TRUE,nrow=1221,header=TRUE)
		#Surveyed. Active. lat long CRN.Rating. USHCN.ID Station.name
		#1        NA      NA       NA          NA    11084 BREWTON 3SSE
	details$watts=NA
	test=match(details$id,watts[,"USHCN.ID"]);temp=!is.na(test); sum(temp) #1221
	details$watts[temp]=watts[  test[temp],5]
		#  0   1   2   3   4   5 
		#  9  17  36  69 215  58 
	details$watts[(details$watts==0)&!is.na(details$watts)]=NA
		#save(details,file="d:/climate/data/station/ushcn/details.tab")
		#  USHCN_ID Station_name  GHCN_ID Surveyed CRN.Rating
		#1    11084 BREWTON 3SSE 74777006       NA         NA

##LOAD USHCN TOBS
	#load("d:/climate/data/station/ushcn/ushcn.collation.tab")
	#name0=dimnames(ushcn.collation[[1]])[[2]]
	#ushcn.tobs=ushcn.collation[["tobs"]]
	#save(ushcn.tobs,file="d:/climate/data/station/ushcn/ushcn.tobs.tab")

	url="http://www.climateaudit.org/data/ushcn/ushcn_tobs.tab" 
  	download.file("http://www.climateaudit.org/data/ushcn/ushcn_tobs.gz","ushcn_tobs.gz",mode="wb"); 
	handle = gzfile("ushcn_tobs.gz");
	ushcn_tobs=read.table(handle,sep="\t",header=TRUE)
     	close(handle);
	dim(ushcn_tobs)  #2544 1222
	ushcn_tobs=ts(ushcn_tobs[,2:ncol(ushcn_tobs)],start=c(ushcn_tobs[1,1],1),freq=12) 
	name0=dimnames(ushcn_tobs)[[2]];dimnames(ushcn_tobs)[[2]]=substr(name0,2,nchar(name0))


##IDENTIFY CRN1-2
	temp=(details$watts<=2)&!is.na(details$watts);sum(temp)
	details[temp,1:10]
		#temp=(details$watts<=2)&!is.na(details$watts)&((details$lights=="dark")|(details$bright==0));sum(temp)
	index_good=(1:1221)[temp]
    		#good=ushcn_tobs[,index_good] #2544 53
		#M=ncol(good);M #53

##RECONSTRUCTION
	recon.crn12=ushcn_wt(index_good)
	ts.plot(recon.crn12)
	X=cbind(time(recon.crn12),recon.crn12);dimnames(X)[[2]]=c("time","recon.crn12")
	write.table(X,file="d:/climate/data/station/watts/recon.crn12.dat",sep="\t",row.names=FALSE)