#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
}
|
---|