#COLLATION FUNCTIONS

    reader.info=function(url){	
	fred=readLines(url)
	fred=fred[substr(fred,1,13)=="</tr><tr><th>"]
	n=nchar(fred)
	fred=substr(fred,14,n)
	#id=as.numeric(substr(fred,1,5));
	#n=nchar(fred)
	#fred=substr(fred,15,n)
	n=nchar(fred)
	temp=(substr(fred,n-4,n)=="</th>")
	fred[temp]=substr(fred[temp],1,n[temp]-5)
	fred=strsplit(fred,"</th><th>")
	fred=t(sapply(fred,rbind))
	fred=data.frame(fred)
	names(fred)=c("id","name","lat","long")
	for(i in 1:4) fred[,i]=as.character(fred[,i])
	fred$id=as.numeric(fred$id)
	n=nchar(fred$lat)
	fred$lat= -as.numeric( substr(fred$lat,1,n-1))
	n=nchar(fred$long)
	temp=substr(fred$long,n,n)=="E"
	fred$long[temp]=substr(fred$long[temp],1,n[temp]-1)
	fred$long=as.numeric(fred$long)
	fred$long[!temp]= - fred$long[!temp]
	temp=(fred$long== -180);
	fred$long[temp]=180
	reader.info=fred
	reader.info
	}


  si_table2f=function() {
	info=read.table("http://www.climateaudit.org/data/antarctic_mann/si_table2.dat")
	dim(info) #46 6 
	temp= info== -999; info[temp]=NA
	names(info)=c("name","lat","long","r","RE","CE")
	info$name=as.character(info$name);n=nchar(info$name)
	info$subset=NA
	temp= substr(info$name,n,n)=="*"
	info$subset[temp]=TRUE
	info$name[temp]=substr(info$name[temp],1,n[temp]-1)
	info$lat= -info$lat
	temp= info$long>180
	info$long[temp]= -( 360- info$long[temp])
	
	test=match(info$name,Info$surface$name);temp=!is.na(test);sum(temp) #40
	info$id=NA
	info$id[temp]=Info$surface$id[test[temp]]
	info[is.na(info$id),]
#             name  lat  long   RE   CE   r2 subset id
#15        Faraday 65.4 295.6   NA   NA   NA   TRUE NA  #89063 
#40 Terra_Nova_Bay 74.7 164.1 0.60 0.33 0.33     NA NA   #89662	Mario_Zucchelli
#43   Mt_Siple_AWS 73.2 232.9 0.58 0.31 0.31     NA NA
#44      Siple_AWS 75.9 276.0 0.70 0.44 0.44     NA NA
#45       Byrd_AWS 80.0 240.0 0.64 0.39 0.39     NA NA
#46      Harry_AWS 73.2 232.9 0.54 0.27 0.27     NA NA
	info$id[c(15,40)]=c(89063,89662) #by inspection ALL   are identified
	si_table2f=info
	si_table2f
	}

 si_table1f=function() {
	info=read.table("http://www.climateaudit.org/data/antarctic_mann/si_table1.dat",header=TRUE)
	#info=read.table("d:/climate/data/mann_antarc/si_table1.dat",header=TRUE)
	dim(info) #26 7
	temp= info== -999; info[temp]=NA
	info$name=as.character(info$name);n=nchar(info$name)
		info$lat= -info$lat
	temp= info$long>180
	info$long[temp]= -( 360- info$long[temp])

	test=match(info$name,Info$aws$name);temp=!is.na(test);sum(temp) #40
	info$id=NA
	info$id[temp]=Info$aws$id[test[temp]]
	info[is.na(info$id),]
#             name   lat  long    r   RE   CE trend id
#6             D10 -66.7 139.8 0.58 0.37 0.27  0.06 NA
#9          Ferrel -77.9 170.8 0.88 0.77 0.75  0.25 NA
#12         LGB120 -73.8  55.7 0.80 0.62 0.62  0.06 NA
#18         Marily -80.0 165.1 0.77 0.59 0.57  0.26 NA
#21   Pegaus_South -78.0 166.6 0.84 0.70 0.69  0.26 NA
#25 Tourm._Plateau -74.1 163.4 0.62 0.09 0.15  0.08 NA
	info$id[c(6,9,12,18,21,25)]=c(89832,89872,89757,89869,8937,7356)
		#all match: stupid spelling errors in Mann SI
	si_table1f=info
	si_table1f
	}


    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
	}


  reader.info2=function(url){	
	fred=readLines(url)
	fred=fred[grep("HREF=temp",fred)]
	n=nchar(fred)
	id=substr(fred,60,63)
	info=data.frame(reader_id=id)	
	fred=substr(fred,70,n)	
	test=strsplit(fred,"<BR>")	
	info$name=sapply(test, function(x) x[1])
	test=sapply(test, function(x) x[2])
	test=gsub("</A>","",test)
	info$id=as.numeric(test)
	reader.info2=info
	reader.info2
	}


#NORMAIZE
	anom=function(x,reference=1961:1990) {
		K=floor(tsp(x)[1])
		if( is.null(dim(x))) {
	 	 month= 1+round(12*  c(time(x))%%1)
		 x=c( rep(NA,month[1]-1),x ,rep(NA,12- month[length(month)]))
		 test=t(array(x,dim=c(12,length(x)/12)) )
		 m0=apply(test[reference-K+1,],2,mean,na.rm=TRUE)
		 test=scale(test,center=m0,scale=FALSE)
		 anom=round(ts(c(t(test)),start=c(K,1),freq=12),2)
		 anom} else {anom=x;
		for(i in 1:ncol(x)) {test=t(array(x[,i],dim=c(12,nrow(x)/12)) )
		m0=apply(test[reference-K+1,],2,mean,na.rm=TRUE)
		test=round(scale(test,center=m0,scale=FALSE),2)
		anom[,i]=c(t(test))}}
		anom
		}
			

###
##LOAD DATA 
###
   #this uses BAS scrape for surface and reverse-engineered Steig recon_aws for AWS. Two surface series need to be deducted (Marion, Gough)
   #version here is not quite file compatible with Steig. See Nic L on Dec 2002 problem
   # data from Terra Nova Bay appears to be same as Mario_Zucchelli #

    make.anomalies=function(method="january",loc="http://www.climateaudit.info/data/steig",download_method="online",mario="dup") {
      if(download_method=="offline") {
 	load("e:/OneDrive/climate/data/steig/Info.tab");
	load("e:/OneDrive/climate/data/steig/Data.tab");
	load("e:/OneDrive/climate/data/steig/recon_aws.tab");
	if(method=="original")  { 
	 test=read.table("e:/OneDrive/climate/data/steig/aws_Steigorigdata.dat",sep="\t",skip=1)
	 names(test)=scan("temp.dat",n=ncol(test),what="")
	 Data$aws=ts(test,start=c(1980,1),freq=12) }
	if(method=="corrected")  { 
	 test=read.table("e:/OneDrive/climate/data/steig/aws_Steigcorrdata.dat",sep="\t",skip=1)
	 names(test)=scan("temp.dat",n=ncol(test),what="")
	 Data$aws=ts(test,start=c(1980,1),freq=12) }
	}
	

      if(download_method=="online"){
 	download.file(file.path(loc,"Info.tab"),"Info.tab",mode="wb")
	load("Info.tab")
	download.file(file.path(loc,"recon_aws.tab"),"recon_aws.tab",mode="wb")
	load("recon_aws.tab")
	dimnames(recon_aws)[[2]]=Info$aws_grid$id
	download.file(file.path(loc,"Data.tab"),"Data.tab",mode="wb")
	load("Data.tab")
	if (method=="original") {
	 download.file(file.path(loc,"aws_Steigorigdata.dat"),"temp.dat")
	 test=read.table("temp.dat",sep="\t",skip=1)
	 names(test)=scan("temp.dat",n=ncol(test),what="")
	 Data$aws=ts(test,start=c(1980,1),freq=12) }
	if (method=="corrected") {
	 download.file(file.path(loc,"aws_Steigcorrdata.dat"),"temp.dat")
	 test=read.table("temp.dat",sep="\t",skip=1)
	 names(test)=scan("temp.dat",n=ncol(test),what="")
	 Data$aws=ts(test,start=c(1980,1),freq=12) }
        }   

   #edit surface 
	surf=Data$surface
	surf=window(surf,start=1957,end=c(2006,12)) #ends in 2006 per Steig
	surf <- surf[,!(dimnames(surf)[[2]]=="68994")] #deletes column 26 - Marion
	surf <- surf[,!(dimnames(surf)[[2]]=="68906")] #deletes column 17 - Gough
	if(mario=="nodup") surf <- surf[,!(dimnames(surf)[[2]]=="89662")] #deletes duplicate Mario/Terra Nova
	ny=ncol(surf) #42
	Data$surface=surf
	Info$surface=Info$surface[!is.na(match(Info$surface$id,dimnames(surf)[[2]] )),]
	sanom=as.matrix(surf)
	if(method=="january") {
		for (i in 1:ny) sanom[,i]=anom(surf[,i],reference=1957:2006) } else {
		for (i in 1:ny) sanom[,i]=anom(surf[,i],reference=1957:2003) } 
            #Steig probably used different reference period and this can be experimented with later

  #AWS reverse engineered (rather than READER scrape)
	if( (method=="original")|(method=="january")) {
	  anoms =window(recon_aws, start = 1980);dim(anoms) #324 63
	  dat_aws = window(Data$aws[,colnames(recon_aws)],start=1980,end=c(2006,12));
	  dim(dat_aws) #324 63
	  anoms[is.na(dat_aws)] =NA
	  nx=ncol(anoms)  } #63

	if( method=="corrected") {
	  dat_aws = window(Data$aws[,colnames(recon_aws)],start=1980,end=c(2006,12));
	  nx=ncol(dat_aws)
	  anoms =as.matrix(dat_aws)
	  for (i in 1:nx) anoms[,i]=anom(dat_aws[,i],reference=1980:2003) } 
   
   # combine anomalies
	anomalies=list()
	anomalies$aws=anoms;anomalies$surface=sanom
	make.anomalies=list(Info=Info,Data=Data,anomalies=anomalies,recon_aws=recon_aws)
	make.anomalies
	}