Steig's Antarctic Peninsula Pac-mann

WUWT readers may recall a couple of weeks ago that I suggested that the weather stations with different climatic influences of the Antarctic peninsula, which might very well merit its own separate climate designation from the Antarctic mainland, was heavily weighting the Steig et al results ( Nature, Jan 22, 2009).  Essentially that weighting “gobbled up” the trends on the mainland, such as the trend at the south pole station which shows a long term cooling.

Jeff Id took that advice and did an analysis which I have reposted by invitation below. But, I just couldn’t help notice that this graph below looks a lot like Jeff’s results .

Above: Peninsula Pac-mann gobbles up the trend. See Figure 8 in Jeff’s analysis.

Antarctic Warming – The Final Straw

Guest posted by Jeff Id of the Air Vent

This is the first post I’ve done which gets to the heart of where the trends in Steig et al. came from. Steve M did a post on TTLS reconstruction TTLS in a Steig Context which makes the point that despite the PCA and truncation the result of RegEM is still a linear recombination of station data. This post is the result of a back calculation of station weights to determine which stations were weighted and by how much to create the final trend of Steig et al.

Before I succeeded in this calculation yesterday, I tried it once before some time ago and it didn’t work. There were a couple of errors which prevented me from getting a solution and I was too lazy to fix it. The Climate Audit post pushed me to try again and this time I got it right. I think you’ll find the result a bit telling.

The satellite reconstruction from Steig et al is based on two halves. The pre-1982 half is entirely surface station data, the post 1982 data is satellite based data. The satellite half is easily replicated from the satellite data while the surface station half is simply a linear weight and sum of the surface stations. If the surface station temperature is SST, and the weights are c the net result of all this complex math prior to 1982 looks like this

T output = (C1 * SST1) + (C2 * SST2) ……. (Cn * SSTn)

That’s it!

So in order to calculate the C’s involved in this equation we can back solve a series of linear equations having the form above. There are 42 SST’s in the reconstruction and 1 Satellite trend. Since the satellite is not used pre-1982 we can ignore that for determining the pre-1982 portion of the reconstruction. So we have 42 SST’s but not all of those have any data before 1982. After removing the stations which don’t have any pre-1982 data only 34 remain. These 34 are the only ones mathematically incorporated in the reconstruction and are shown in Figure 1.

Figure 1 – Location of 34 Stations Used in Reconstruction

It’s odd that Steig et al included the extra stations at all. I’m not sure if they understood what they were doing when they included stations which had no data in the pre-1982 timeframe. I need to run RegEM without them to see for sure but they may affect the weightings of the other 34 stations but IMO it isn’t likely to be helpful.

The code to perform the reconstruction and sort the correct 34 stations out is as follows:

#perform RyanO SteveM RegEM reconstruction

clip=form.steig

dat=window(calc.anom(all.avhrr), start=c(1982), end=c(2006, 12))

base=window(parse.gnd.data(all.gnd, form=clip), start=1957, end=c(2006, 12))

base=calc.anom(base)

pcs=get.PCs(dat,3)

dd=ts.union(base,pcs[[1]])

reg3=regem.pttls(dd,maxiter=50,tol=.005,regpar=3,method=”default”, startmethod=”zero”, p.info=”Unspecified Matrix”)

dim(reg3[[35]]$X) #600 45

#extract surfacestations and PC’s

regemSST=ts(reg3[[35]]$X[,1:42],start=1957,deltat=1/12)

dim(regemSST) #600 42

regemPC=ts(reg3[[35]]$X[,43:45],start=1957,deltat=1/12)

#calculate full reconstruction

recon = regemPC %*%t(pcs[[2]])

recontr=ts(rowMeans(recon),start=1957,deltat=1/12)

coef(lsfit(time(recontr),recontr))[2] #0.01190505

##find stations which have data pre 1982

mask=colSums(!is.na(base[1:300,]))!=0

sum(mask) #34

sstd=base[1:300,mask]

reconSST=base[,mask]## these stations are actually used in recon

nam=as.character(all.idx[1:42,1][mask])

lats=all.idx[[2]][1:42][mask]

lons=all.idx[[3]][1:42][mask]

After the 34 stations are sorted the task is to set up a matrix which has the form of the equation above.

c1 * SST1(x) + C2 * SST2 (x) …… = output(x)

Where x is the value of each surface station and RegEM output on that particular date. Since we have 34 unknowns we need 34 independent equations to solve. All the SST data has values infilled for all dates from 1957 – 2007 but the infilled values are combinations of the non-infilled values. This makes the matrix singular and indeterminate (unsolvable). Our task then is to find 1 row (date) for each station for which the station has have at least 1 unique measured value. To do this I used the raw data and looked for independent months which contain at least 1 value for each row. (this is where I got lazy last time)

##backsolve regem weights

##find unique rows which have 1 value for each station

index=array(0,dim=34)

for(i in 1:34)

{

j=1

while( (is.na(sstd[j,i]) == TRUE) | (sum(index==j)!=0) )

{

j=j+1

}

index[i]=j

}

##use index rows to backsolve RegEM: Index =

# [1] 65 1 109 2 135 3 52 26 4 25 5 6 7 8 9 171 10 165 148 11

# [21] 12 13 74 292 50 73 14 240 280 275 15 16 27 17

##setup square matrix a from infilled data

The value index listed in the code above is the row (month) number from jan 1957 = 1 forward for which at least 1 value was measured. You can see the first station on the list has a value of month 65 for the starting value, the second has a value of 1 which means the second station has data for the first month. The fourth station has a value in the first month but we can’t reuse the same value or the matrix would be singular so it found the next open value at month 2. The algorithm continues in this fashon through the 34 stations.

After these values are gathered we can set up the matrix and solve the following equation for c.

a * c = b

I like simple. The code looks like this.

##setup square matrix a from infilled data

a = regemSST[index,mask]

dim(a) #34,34

b=recontr[index]

c= solve(a,b)

aa=regemSST[,mask]

#a%*%c

m = aa %*% c

m=ts(m,start=1957,deltat=1/12)

The matrix aa is multiplied times weights c to create the surface station temperature reconstruction m. Here is the replicated trend by RegEM we’ve seen before, thanks primarily to Ryan O and SteveM code.

Figure 2

For the first time we can see the Steig et al reconstruction as determined by the surface station temperatures only.

Figure 3

I was a bit shocked the trend was still so high. After all we know the area weighted surface station trend sits at about 0.04 C/Decade.

Just to make it clear, Figure 4 is the difference between the above plots.

Figure 4

The pre-1982 data is a perfect match up to rounding error the post 1982 difference is the satellite data difference which I have to point out boosts the final recon trend a bit higher than the weighted surface stations. The surface stations and weights “C’s” required to recreate the pre – 1982 Steig reconstruction are in Table 1.

Now we get to the fun part. Surface station weights for this reconstruction are shown in Figure 5. The graph is color coded the same as Figure 1 by region. I’ve moved Byrd from Ross Ice shelf to West Antarctica which is the only change from Ryan O’s color coding in his posts.

Figure 5

You can see the dominant number of (black) surface stations located in the peninsula. The Y axis is normalized to 1 equals 1/34 of the total contribution for 1 in 34 stations. This area is of course known to have high warming trend, however 4 stations have strong negative net weights – an oddity I mentioned in my earlier work on this paper explaining RegEM ignores trend in favor of high frequency correlation. It is of course nonsensical to flip temperature data upside down when averaging but that is exactly what Steig et al does. This alone should call into question the paper’s result.

This isn’t the end of the story however, in Figure 6 I multiplied the individual (infilled by RegEM) station trends times their weights and created another bar plot

Figure 6

Ok, at this point my eyes are widening. Figure 6 represents the contribution of each stations trend to the positive total output trend. Negative values here are acceptable if they come from negative trend, so the 4 black bars and one near zero blue which were negative in Figure 5 are incorrect, and the ones which changed sign for Figure 6 are a result of a truly negative trend in temperature.

Figure 7 is a Pie chart showing the station weights for each region- same as Figure 5 – different colors.

Figure 7

It’s telling in Figure 7 that station weights for the tiny peninsula region were not contained well spatially in that the sum of the weights adds up to an area equal to the entire East Antarctica. A correct reconstruction would contain this information to a section of the pie reasonably equivalent to the geographic area of coverage.

And finally the graph we’ve all been looking for since this all started, the contribution of each region to the total reconstruction trend.

Figure 8

There it is, we can now say conclusively that the positive trend in the Antarctic reconstruction comes primarily from the well known peninsula warming trend.

If we recall Figure 3 is the actual Steig et al reconstruction using both pre and post 1982 surface station data only and yet the trend is nearly the same as the final RegEM. This trend is quite different from simple methods of determining station weights using methods such as these.

Maximum Triviality Reconstruction

Closest Station Antarctic Reconstruction

My final check was to add up the area contribution to trends as a check. These values created Figure 8. The four values in order are from Peninsula, West Antarctica, Ross Ice Shelf, East Antarctica in degrees C/Decade:

0.0709 + 0.0115 + 0.0028 + 0.0134 = 0.0987

This was in fact an exact match (7 figures) of the trend in Figure 3 above. Demonstrating the correctness of the last equation in Steve McIntyre’s CA post linked above.

It will be interesting to see how well RyanO’s latest holds up to the same analysis – don’t expect any favoritism around here .

NOTE: be sure to read Steig et al falsified as well. Real Climate has post direct from Dr. Steig on why they don’t want to discuss issues like this one. – Anthony