Antarctic Temperature RegEm Forensics

Thanks to Molan Labe who was quite helpful in identifying what I believe is the correct coefficient into regparm for the purposes of recreating the RegEm calculations. After some considerable time I am still not able to recreate the results of Steig’s paper perfectly, but as you can see below, I am fairly close.

First let’s look at the mean trend of the reconstructed AWS station data by Steig. All data used was from the Climate Audit directory (thanks to Steve McIntyre) so the updated reconstructions and data are not involved.

The slope is written in the bottom right. Next is my attempt at the reconstruction.

Overall the shapes are visually quite similar so I subtracted them.

The step function is quite visible. The data pre-1980 is completely imputed/extrapolated while after 1980 has at least some content from one automatic weather station or other. Next are some of the differences in individual series from the data Eric Steig archived for Nature and my reproductions.

The flat line sections are sections where the actual data exists in each series, the other data was infilled by RegEm.  The step in the difference is quite visible here as well.  I did take the time to confirm that offsetting the surface station data had no effect on the output reconstruction whatsoever so an anomaly in the offset is not the issue.   I will try some other settings in RegEm but I doubt that this will resolve the difference.  All in all, I feel like I’m pretty close now.  Perhaps someone has a suggestion which could resolve the final differences.

Here is the code I used to make the above reconstructions in R.  If you cut and paste it you will need to fix the wordpress quotes.

download.file(“http://www.climateaudit.org/data/steig/Data.tab”,”Data.tab”,mode=”wb”)
load(“Data.tab”)
download.file(“http://www.climateaudit.org/data/steig/Info.tab”,”Info.tab”,mode=”wb”)
load(“Info.tab”)
download.file(“http://www.climateaudit.org/data/steig/recon_aws.tab”,”recon_aws.tab”,mode=”wb”)
load(“recon_aws.tab”)
dimnames(recon_aws)[[2]]=Info$aws_grid$id

#########backcalc anomalies
anoms =window(recon_aws, start = 1980)
dat_aws = window(Data$aws[,colnames(recon_aws)],start=1980,end=c(2006,12))
anoms[is.na(dat_aws)] =NA
#write.csv(anoms,file=”c:/agw/antarctic paper/bc anoms1.csv”)

#########calculate surfacestatoin anoms
surf=Data$surface
surf=window(surf,start=1957,end=c(2006,12))
a=as.numeric(surf)
dim(a)=c(600,44)
surfanomaly=array(NA,dim=c(600,44))
for(i in 1:44)
{
b=a[,i]
dim(b)=c(12,50)
an=rowMeans(b,na.rm=TRUE)

surfanomaly[,i]=a[,i]-an
}
#write.csv(surfanomaly,file=”c:/agw/antarctic paper/surface temp anoms.csv”)
sanom=ts(surfanomaly,start=1957,deltat=1/12)
##### combine anomalies

anomalies=ts.union(anoms,sanom)

write.csv(anomalies,file=”c:/agw/antarctic paper/combined anomalies r4 offset.csv”)
######## after matlab RegEm run this #################
recon_id=read.csv(file=”c:/agw/antarctic paper/recon redo R4 offset e=3 rp=3.csv”,header=FALSE)

recon_id_aws=ts(recon_id[1:63],start=1957,deltat=1/12)

#######plot differences

for(i in 1:63)
{
plot(recon_aws[,i]-recon_id_aws[,i],main=paste(“RegEm Steig – IdnStation # = “,colnames(recon_aws)[i]),xlab=”Year”,ylab=”Anomaly C”)
savePlot(paste(“C:/agw/antarctic paper/id recon difference”,colnames(recon_aws)[i],”.jpg”),type=”jpg”)
}

#######total trend calcs
idmean=ts(rowMeans(recon_id_aws),start=1957,deltat=1/12)
steigmean=ts(rowMeans(recon_aws),start=1957,deltat=1/12)

plot(idmean,main=”Antarctic AWS Recon TrendnId eigs=3 regpar=3″,xlab=”Year”,ylab=”Anomaly C”)
abline(lsfit(time(idmean),idmean,NULL), col=”red”)
text(1962,-6,paste(“Slp =”,round(  coef(lsfit(time(idmean),idmean,NULL))[2],3)  ),col=”red”)
savePlot(“C:/agw/antarctic paper/id recon trend A.jpg”,type=”jpg”)

plot(steigmean,main=”Antarctic AWS Recon TrendnSteig eigs=3 regpar=3″,xlab=”Year”,ylab=”Anomaly C”)
abline(lsfit(time(steigmean),steigmean,NULL), col=”red”)
text(1962,-6,paste(“Slp =”,round(  coef(lsfit(time(steigmean),steigmean,NULL))[2],3)  ),col=”red”)
savePlot(“C:/agw/antarctic paper/steig recon trend A.jpg”,type=”jpg”)

plot(steigmean-idmean,main=”Antarctic AWS Recon Trend DifferencenSteig – Id eigs=3 regpar=3″,xlab=”Year”,ylab=”Anomaly C”)
savePlot(“C:/agw/antarctic paper/steig minus id recon trend A.jpg”,type=”jpg”)

######recenter anomaly calculations based on 1980 and redo analysis

sanomw=window(sanom,start=1980)
means=colMeans(sanomw,na.rm=TRUE)

for(i in 1:44)
{
if(!is.na(means[i])){    sanom[,i]=sanom[,i]-means[i]}
}

######

anomalies=ts.union(anoms,sanom)
write.csv(anomalies,file=”c:/agw/antarctic paper/combined anomalies r2.csv”)

Since most of the work was done in R, the matlab code was quite simple.

A = dlmread(‘c:/agw/antarctic paper/combined anomalies r3 offset processed.csv’, ‘,’)

OPTIONS.regress = ‘ttls’;
OPTIONS.neigs=3
OPTIONS.regpar=3
A=regem(A,OPTIONS)

dlmwrite(‘c:/agw/antarctic paper/recon redo R4 offset e=3 rp=3.csv’,A)