#COLLATION FUNCTIONS reader.info=function(url){ fred=readLines(url) fred=fred[substr(fred,1,13)==""] 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)=="") fred[temp]=substr(fred[temp],1,n[temp]-5) fred=strsplit(fred,"") 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,"
") info$name=sapply(test, function(x) x[1]) test=sapply(test, function(x) x[2]) test=gsub("","",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 }