#COLLATION
#READER http://www.antarctica.ac.uk/met/READER/
#http://www.antarctica.ac.uk/met/READER/surface/stationpt.html
#10 Australian http://www.antarctica.ac.uk/met/READER/temperature.html
#http://www.antarctica.ac.uk/met/READER/aws/awspt.html
####STATION INFO

  #1. Read Bristish Antarctic Survey READER metadata
	#7 UK: http://www.antarctica.ac.uk/met/READER/metadata/metadata.html
	#10 Russian: http://www.aari.aq/default_en.html
	#4 Australian: http://www.antarctica.ac.uk/met/READER/metadata/Australia.html

   source("http://www.climateaudit.org/data/antarctic_mann/collation.functions.txt")
	Info=list()
	Info$surface=reader.info("http://www.antarctica.ac.uk/met/READER/surface/stationpt.html")
	Info$surface$name[16]=substr(Info$surface$name[16],1,7) #faraday name
	Info$aws=reader.info("http://www.antarctica.ac.uk/met/READER/aws/awspt.html") #65
	Info$si_table2=si_table2f()
	Info$si_table1=si_table1f()
	Info$ocean=reader.info2(url="http://www.antarctica.ac.uk/met/READER/temperature.html")
	temp= (Info$aws$long== -180);	Info$aws$long[temp]=180;	


    #STEIG  #this is 63 series 
     #http://faculty.washington.edu/steig/nature09data/ant_recon_aws.txt
	download.file("http://faculty.washington.edu/steig/nature09data/data/ant_recon_aws.txt","temp.dat",mode="wb")
	grid=scan("temp.dat",n= -1) #  37800
	  #600*63 #[1] 3305400
	recon_aws=ts(t(array(grid,dim=c(63,600))),start=c(1957,1),freq=12)	
	dimnames(recon_aws)[[2]]=Info$aws_grid$id
	lat=scan("http://faculty.washington.edu/steig/nature09data/data/aws_lats.txt")
	long=scan("http://faculty.washington.edu/steig/nature09data/data/aws_lons.txt")

	Info$aws_grid=data.frame(id=1:63,lat,long)
	temp=(Info$aws_grid$long>180);
	Info$aws_grid$long[temp]= -( 360- Info$aws_grid$long[temp])

	Info$aws_grid$id=NA
	for(i in 1:nrow(Info$aws_grid) ) {
	 lat0=Info$aws_grid$lat[i];long0=Info$aws_grid$long[i]
	temp1=abs(Info$aws$lat-lat0)<.01;sum(temp1)
	temp2=abs(Info$aws$long-long0)<.01;sum(temp2)
	if (sum(temp1&temp2)==1) Info$aws_grid$id[i]=Info$aws$id[temp1&temp2]
	}
	Info$aws_grid$id[8]=89810  #match manually (two versions)
	test=match(Info$aws_grid$id,Info$aws$id);temp=!is.na(test);sum(temp)
	Info$aws_grid$name[temp]=Info$aws$name[test[temp]]
   
##GISS

	load("d:/climate/data/station/giss/giss.info.tab")
		#from 7234 to 7350

 #antarc1 are 44 READER surface 
	url="d:/climate/data/station/giss/code/200802/STEP0/input_files"
	target=read.table(file.path(url,"antarc1.list")) #44 rows
	names(target)=c("id","name","lat","long","alt","code")
	test=match(Info$surface$name,target$name);temp=!is.na(test);sum(temp)	#46
	Info$surface$gissid=NA;
	Info$surface$gissid[temp]=target$id[test[temp]]
	Info$surface$gissid[c(16,28  )]= target$id[c(23,40)] #manual match of Faraday; Zucchelli
	Info$surface[is.na(Info$surface$gissid),]
	#      id               name   lat  long gissid
	#18 68906              Gough -40.4  -9.9     NA
	#27 68994             Marion -46.8  37.8     NA
	test=match(target$id,Info$surface$gissid);temp=!is.na(test);sum(temp)	#44 #all matched
	Info$surface$gissid[18]=giss.info$id[grep("GOUGH",giss.info$site)]
	  #giss.info[grep("MARION",giss.info$site),1:6]
	Info$surface$gissid[27]= giss.info$id[521]
	sum(is.na(Info$surface$gissid)) #0

 #antarc2 southern Ocean
	url="d:/climate/data/station/giss/code/200802/STEP0/input_files"
	target=read.fwf(file.path(url,"antarc2.list"),widths=c(11,1,30,7,8,5)) #44 rows
	target=target[,c(1,3:6)]
	names(target)=c("id","name","lat","long","alt");dim(target)
	target$name=gsub(" +$","",target$name)
	target$name=gsub("ISLAND","IS.",target$name)
	target$name=gsub("IS)","IS.)",target$name)

	Info$ocean$gissid=NA;
	test=match(Info$ocean$name,target$name);temp=!is.na(test);sum(temp)	#46
	Info$ocean$gissid[temp]=target$id[test[temp]]
	Info$ocean[!temp,]
	#   reader_id                      name    id gissid
	#2       brow                     BROWN 88971     NA
	#9       sana                S.A.N.A.E. 89001     NA
	#21      pasc ISLAS PASCUA (EASTER IS.) 85469     NA
	#23      kerg             KERGUELEN IS. 61988     NA
	#35      tris          TRISTAN DA CUNHA 68902     NA
	Info$ocean$gissid[c(21,23,35,9,2)]=target$id[c(8,6,2,30,29)] #manual 
	test=match(target$id,Info$ocean$gissid);temp=!is.na(test);sum(!temp) #0 #all 35 matched

  #antarc3 
	url="d:/climate/data/station/giss/code/200802/STEP0/input_files"
	target=read.table(file.path(url,"antarc3.list")) #66 rows
	names(target)=c("id","name","lat","long","alt","code")
	test=match(Info$aws$name,target$name);temp=!is.na(test);sum(temp)	#65
	Info$aws$gissid=NA;
	Info$aws$gissid[temp]=target$id[test[temp]]  #65 assigned
	test=match(target$id,Info$aws$gissid);temp=!is.na(test);sum(!temp) #1 #all 35 matched
	target[!temp,]
	#         id name   lat  long  alt    code
            id           name   lat  long  alt    code
	#23 70089269000 Bonapart_Point -64.8 -64.1    8 antarc3
	#54 70089833900           D_17 -66.7 139.7 -999 antarc3
	grep("17",giss.info$site) #7336 match
 	Info$aws[is.na(Info$aws$gissid),]
	#     id            name   lat  long gissid
	#1 89269 Bonaparte_Point -64.8 -64.1     NA
	grep("Bonap",giss.info$site) #7290 match
 	Info$aws$gissid[1]=giss.info$id[7290] #manual 


		#this is in http://www.antarctica.ac.uk/met/READER/aws/D_17.All.temperature.html
		#but only a few values in 1980

	#save(Info,file="d:/climate/data/steig/Info.tab")

 	sapply(Info,nrow)
 	#surface       aws si_table2 si_table1  aws_grid     ocean 
        #    46        65        46        26        63        35 













	save(Info,file="d:/climate/data/antarctic_mann/Info.tab")

	Data=list()
	Data$aws=chronf("aws")
	  dim(Data$aws) #360 65
	  tsp(Data$aws) #  1980.000 2009.917   12.000

	Data$surface=chronf("surface")
	  dim(Data$surface) #1284   43
	  tsp(Data$surface) #  1903.000 2009.917   12.000
	#save(Data,file="d:/climate/data/mann_antarc/Data.tab")
	#Data.old=Data

	apply(	

	save(Data,file="d:/climate/data/steig/20090210.Data.tab")

	sapply(Info, function(A) A[A$wa,] )


##UPDATES
	id="aws"
	i=28
        loc=paste("http://www.antarctica.ac.uk/met/READER/",id,"/",Info$aws$name[i],".All.temperature.txt",sep="")
         test=try(read.table(loc,skip=1,na.strings="-"))
         if(! (class(test)=="try-error") )test=ts( c(t(test[,2:13])),start=test[1,1],freq=12) else test= ts(rep(NA,12),start=c(2000,1),freq=12)
       	test;
	length(test)
	temp=(time(Data$aws)>=tsp(test)[1]) & (time(Data$aws)<=tsp(test)[2]); sum(temp)
	Data$aws[,28]=NA; Data$aws[temp,28]=test

	ts.union(recon_aws[,paste(id0)],anom(Data$aws[,28],index=1980:2000))

       
 dimnames(chron)[[2]]=info$id
        #count=ts( apply(!is.na(chron),1,sum),start=c(tsp(chron)[1],1),freq=12)
        count=apply(!is.na(chron),2,sum) # 3 are zero
        chronf=chron[,count>0]
 



##RECOLLATE FEB 10
 Info.old=Info
    chronf=function(id) {	
    id="aws"
	info=Info[[paste(id)]]
	chron=NULL
	for (i in 1:nrow(info)) {
	 loc=paste("http://www.antarctica.ac.uk/met/READER/",id,"/",info$name[i],".All.temperature.txt",sep="")
	 test=try(read.table(loc,skip=1,na.strings="-"))
	 if(! (class(test)=="try-error") )test=ts( c(t(test[,2:13])),start=test[1,1],freq=12) else test= ts(rep(NA,12),start=c(2000,1),freq=12)
	 chron=ts.union(chron,test)
	}
	dimnames(chron)[[2]]=info$id
	#count=ts( apply(!is.na(chron),1,sum),start=c(tsp(chron)[1],1),freq=12)
	count=apply(!is.na(chron),2,sum) # 3 are zero
	chronf=chron[,count>0]
	chronf
	}

	info=reader.info("http://www.antarctica.ac.uk/met/READER/aws/awspt.html") #65
	test=ts.union(anom.aws[,1],recon_aws[,1])
	ts.plot(test[,2:1],col=1:2)
	fred=anom(Data$aws[,1],index=1980:2000)	


	test=ts.union(Data.old$aws[,1],chron[,1],recon_aws[,1])
	   #vhanged
	for(i in 1:2) test[,i]=anom(test[,i],index=1980:2000)

   GDD(file="d:/climate/images/2009/steig/bonapart.gif",type="gif",w=480,h=360)
	ts.plot(test[,3],ylab="deg C anomaly")
	points(c(time(test)),test[,2],pch=19,cex=.5,col=4)
	points(c(time(test)),test[,1],pch=19,cex=.5,col=2)
	lines(c(time(test)),test[,1],col=2)
	title("Bonapart/e Point AWS")
	legend("topleft",fill=c(1,2,4),legend=c("Recon Old","First Scrape","Second Scrape"),cex=.9)
    dev.off()
		

  GDD(file="d:/climate/images/2009/steig/bonapart.gif",type="gif",w=480,h=360)
	ts.plot(recon_corr[,1],ylab="deg C anomaly")
	points(c(time(test)),test[,2],pch=19,cex=.5,col=4)
	points(c(time(test)),test[,1],pch=19,cex=.5,col=2)
	lines(c(time(test)),test[,1],col=2)
	title("Bonapart/e Point AWS")
	legend("topleft",fill=c(1,2,4),legend=c("Recon Old","First Scrape","Second Scrape"),cex=.9)
    dev.off()


	test1=apply(!is.na(chron),2,sum)
	test2=apply(!is.na(Data.old$aws),2,sum)

	x=  !(chron==Data.old$aws) &!is.na(chron)&!is.na(Data.old$aws) 
	test3= apply(x  ,2,sum)



	Data=list()
	Data$aws=chronf("aws")
	info=Info[[paste("aws")]]
	  dim(Data$aws) #360 65
	  tsp(Data$aws) #  1980.000 2009.917   12.000

	Data$surface=chronf("surface")
	  dim(Data$surface) #1284   43
	  tsp(Data$surface) #  1903.000 2009.917   12.000
	#save(Data,file="d:/climate/data/mann_antarc/Data.tab")
	#Data.old=Data

	apply(!is.na(Data


	WA=data.frame( lat=-c( 60, 85, 85,   85, 85, 85,  60,  60,  60,  60,  60,  70),
	       long=c(178,178,-150,-120,-90,-56, -56, -56, -90,-120,-150, 178) )
	for (i in 1:4){
	temp1=(Info[[i]]$lat> -85)
	temp2=(Info[[i]]$long < -56) | (Info[[i]]$long>=178)
	Info[[i]]$wa=temp1&temp2
	}
	for(i in 1:4) {
	 temp1=(Info[[i]]$lat> -85)&(Info[[i]]$lat<= -79)
	 temp2=(Info[[i]]$long> 160)
	 Info[[i]]$wa[temp1&temp2]=TRUE
	}

#concordance Wash, READER


    #READER
	download.file("http://www.climateaudit.org/data/antarctic_mann/Info.tab","temp.dat",mode="wb");load("temp.dat")
	download.file("http://www.climateaudit.org/data/antarctic_mann/Data.tab","temp.dat",mode="wb");load("temp.dat")