The Best Laid Plans

My job for the Antarctic paper so far has been simple. Present an area weighted reconstruction of the surface stations as a simple sanity check of the serious math reconstruction Ryan and Nic have been doing and I’ve been reading along on by email. So if you’ve been visiting once in a while you know I’ve presented several area weighted reconstruction’s. The algorithm I prefer so far is one which infill’s each satellite gridcell by the closest station which has data for that year I’ve created several ‘dozen’ versions, the latest one worth publishing was shown Some Other Area Weighted Reconstrucions and others were presented at Area Weighted TPCA Check .

Since the latest reconstructions utilize the noisy automatic weather stations which weren’t employed until after 1982 there are a lot of new series starting in the second half of the data. This creates some problems. First, we need to recognize that the Manned stations are more consistent and of a better quality than AWS stations which can become buried in snow for extended periods of time. A second issue is also important, each anomaly is centered about it’s mean. Average = 0. So if you have two thermometers which measure exactly the same values always and one starts 10 years after the other. Both have a mean of all their data of zero so if any trend exists thermometer 2 which measures exactly the same temperature as thermometer 1 ends up having a slightly different anomaly.  When the two are averaged this creates a sudden step in the data.

Currently all presented plots do not correct for this factor. Since trends are fairly low, it isn’t a terrible method to use and we can expect that whatever corrected versions I come up with we’ll have a similar result. I’ve had several ideas of how to fix this annoying effect, some more stupid than others and none of which were worth posting on.

One particularly misguided method I used was to start at 1957 – the beginning year of our reconstruction – and offset existing Jan 1957 temps by zero. Then each time a new series came in I took the mean of all the other stations and applied that as the offset for the new data. I got a trend of 0.12 C/Decade which is out of whack completely and matches Steig et al. We know it’s out of whack by checking the surface station versions of the same algorithm which typically come to 0.04 or 0.05C/Decade. What the algorithm did by taking mean was to ignore the fact that the peninsula and Ross ice shelf are heavily oversampled in relation to the rest of the continent.

Ryan O made a great point that the offset in this manner assumes that when a new station is introduced and the mean is used, the temperature record for the new station prior to its existence is assumed to be the same as the continent. It’s an obviously incorrect assumption. In this post, I tried another method for correction. Since any offset assumes the trend of the new station prior to its own existence the next best option would be to choose the next closest station with existing data to compute the offset!

Well the algorithm I used was designed to calculated the offset between the new station and next closest station by using an average of the closest data points. I took the matrix of all ground station data on the sat grid 63 manned and aws stations and ran them through a filtering algorithm.  This matrix FYo (filtered Yo)was stored separately and used only for calculating offsets.

squarefilter=function(Yo=Yo,windowsize=3)

{

FYo=Yo #allocate memory

for(j in time(Yo))

{

FYo[time(Yo)==j,]= colMeans(window(Yo,start=j-windowsize,end=j+windowsize),na.rm=TRUE)

}

FYo

}

FYo=squarefilter(Yo,2)

The last line calls the function with a half width of two years. The filter does no endpoint projection and simply averages the data inside the specified window. Seemed like a good idea, however the data is so noisy that I got some unusual results. The function to go through the series wasn’t terribly complicated but isn’t exactly simple either.

for (i in 1:600)

{

mask= is.na(offset)

newmask=((!is.na(Yo[i,]))-!mask)>0

if (sum(newmask)>0)

{#new value started

for(j in which(newmask==TRUE))

{

#get next closest distance with offsets

dist =circledist(surfcgrid[j,1],surfcgrid[j,2],surfcgrid[,1],surfcgrid[,2],)

hh=order(dist[!mask],decreasing=FALSE) #order all stations which have offsets

namask=is.na(Yo[i,hh]) #find stations with no data values

hh=hh[!namask][1] #remove stations with no values

stationindexb=which(dist==dist[!mask][hh])[1]

offset[j]=FYo[i,stationindexb]+offset[stationindexb] – FYo[i,j] #calculate offset as a difference between

#filtered anomalies

if(j==4)

{

print(stationindexb)

print(time(Yo)[i])

}

}

}

}

Formatting by WordPress.

The function moves through all surface stations, as new stations are introduced it calculates an offset based on the filtered existing data from the next closest station. It took quite a while to write and debug so I expected interesting results and wasn’t terribly dissapointed to that end. The following plots are from a randomly chosen gridcell #130 of 5509 total. Each color represents a different surface station. Each satellite gridcell is infilled with the closest available surface station data so the introduction of new stations over time results in a complex combination of surface measurements for each point. Actual stations used are enumerated in the legend.

These plots are of  1 of 5509 reconstructed satellite gridcells.

Figure 1 - One year square filter offset

The light blue surface station #24 is consistently below station #64 in the trend. Of course this is completely unnatural and means that my offsets are too noisy. I confirmed that by trying a variety of different filter levels which demonstrated the offsets were better and worse depending on the filter level. The problem is that in order to get values that look reasonable, the filter level had to get ridiculously long to reduce any high frequency noise.

Below is a plot of he same values using a 4 year filter windows.

Figure 2 - Four year offset filter

The four year window results in much better ‘looking’ values but the trend of the reconstruction is dependent on the filter level. For the four year trend, this method calculated — 0.083C/decade. I believe this trend to be about 60%-80% higher than actual from the manned surface station only reconstructions. The same plot without any offsets is in Figure 3. This is still my preferred area weighted reconstruction despite the fact it has offset issues described in the beginning.

Figure 3, Station offset =0

The four year filtered area weighted plot looks like this:

Figure 4. Four year filter area weighte plot. Trend = 0.083C/Dec

Compared Figure 4 to the non-offset version which has a trend of 0.051 C/Decade in Figure 5.

Figure 5, Antarctic temperature trend no offset

Well, all I can say is another 10 hours spent on a great idea which didn’t pan out. For those of you who are perennially on the wrong side of every issue, several trends I found were in the opposite direction from reality – cooling. The slightly higher trend, of the offset reconstruction is perfectly acceptable and that was the one I presented, however the method just didn’t work.

So with all this, I’ve got a new idea. My next version will take all the available points between this station and the next closest station and regress them to determine the offset. This will use the maximum available information to align the two stations. It’s Sunday though and I’ve spent four hours this morning and yesterday on the Antarctic trend and I’m going to take a break.  BTW, We’re not even allowed to visit anymore due to the possibility of environmental damage (what a pile of bull).  The Antarctic is huge in area and impossible to damage through a bit of tourism.  Not that I wanted to go anyway but it would have been fun for a short time.