Wednesday, March 30, 2011

What's it like to run out of data?

In my previous post, I described a controversy which had been created about Keith Briffa's apparent decision to start a plot of paleo reconstructed NH temperature at 1550, when some earlier data was available. Plotting the data earlier than 1550 showed increasing variations, very likely due to a rapid diminution in the supporting data. I showed plots of the number of sites available, and indeed they did reduce rapidly between 1550 and 1400, when there were only 8 left.

There are at least two reasons why the reduction in data might lead to spurious oscillations. One is simple averaging - fewer sites means greater variance in the mean. The other is geographical spread. I might be able to say something later about spread, but for this post I'm looking at the averaging effect.

I simulated with some artificial data on the actual site histories (from Briffa's 1998 paper). I took a 40-year sinusoid and added AR1 noise (with a 5-year time constant). And I plotted the results below the jump.

Update - I should add that another effect, maybe as important, is the loss of replication within sites. This analysis only covers reduction in sites.

Keith Briffa and the Renaissance

While I was busy with the Antarctica analysis a blogfuss started at CA on Briffa's 1999 (!) Science presentation of various paleo reconstructions, particularly relating to whether a particular graph should have started in 1400 or 1550. A file had been discovered which showed data down to 1400, and if you plot it, it goes into oscillations in the years before 1550. Since it is clear that this is in a period of rapidly diminishing data, and very likely caused by that, I thought that would die fairly quickly, but no, as these things go, it was promoted to a grand ethical violation, megaphoned at WUWT, and taken up at the Air Vent, where it was seen as "unbelievable fraud".

The evidence is supposed to be this graph, taken from this paper: SEEING THE WOOD FROM THE TREES. KR Briff and TJ Osborn, Science, 1999

The pink bits have been helpfully added by Steve McIntyre, based on numbers he found in an XL file which had been apparently inadvertently attached to another paper's data archived at NOAA. It's undocumented there.

Well, it seemed clear to me that the available data is just getting low as we go back beyond 1550, and the wild swings are just the result of the growing noise, as you'd expect. And I haven't found anyone who seems to seriously think they reflect any kind of reality. So Briffa sensibly stopped at 1550 to avoid misleading the public. But no, it's apparently academic misconduct, another Yamal even. [ed.. But, but, wasn't that about actually proceeding with plots where data was running low? O, never mind, ethical violation etc].

So I belatedly saddled up and sallied forth to dispute at CA and tAV. I don't want to rehash that here, nor transfer the debate, which I'll continue there. In the course of it, there was argument about how much the data was actually reducing, and Steve McI pointed me to the archived list of sites at UEA. There isn't one for the Science paper but there is a list for the Nature paper that it references:
Briffa KR et al. (1998) Influence of volcanic eruptions on Northern Hemisphere summer temperature over the past 600 years. Nature 393, 450-455.
and a list for a subsequent paper
Briffa KR et al. (2001) Low-frequency temperature variations from a northern tree-ring-density network. J. Geophysical Research 106, 2929-2941.

From these, I could extract a graph for the rate at which the number of sites was reducing in the relevant period. I also made lists of the actual sites concerned.

Thursday, March 24, 2011

Trends in Antarctica

In a comment, Eric Steig suggested direct estimation of temperature trends in West Antarctica(WA) (as opposed to estination via fitted EOFs). For trend estimation, the methods of S09 and O10 have good and bad points. A plus is that they do estimate for a prescribed area. Against that is that it isn't obvious, with variable sea ice, what area to prescribe. A minus is the indirectness. Temperatures are used to fit basis functions (EOFs), and the trend is obtained from the fitted function. Whether through a small number of EOFs, or via regularisation (kgnd), this interpolated step restricts degrees of freedom.

TempLS is designed for the direct estimation of trends via least squares model fitting, so I thought I would try that direct approach. I deferred the project while I looked for a suitable area weighting scheme, and that is now available.
Background to this post can be found in earlier discussions here, here, here, here, here, and here

The discussion is prompted by two papers:
S09, a 2009 Nature paper by Steig et al and RO10 (or O10), a 2010 J Climate paper by O'Donnell et al. There has been heated blog controversy about these papers - this post is just one entry point. Here is a Real Climate post following S09's original appearance, and here is Eric Steig commenting on O10 (this led to other heated postings).
There are earlier discussions of area weighting by RO10 authors at tAV, eg here and here.

The results below broadly confirm the earlier spatial analyses, and follow the regional breakdown of S09 and O10, though often with higher trends. The continent ground station trend, 1958 to 2009 and not area weighted, was 0.175 ± 0.03 °C/decade. This is consistent with my earlier EOF-based result, and higher than S09 or O10. Area weighting brought it down to 0.114 ± 0.035 °C/decade, similar to S09 (with AVNRR). Adding equally weighted AVHRR pushed it up to 0.137 ± 0.037 °C/decade. I think these results add a degree of confirmation to those earlier analyses, and to my EOF-based analysis, and tend to confirm a positive warming trend, especially outside East Antarctica.

Something that needed attention is the regional tagging of the data I am using, which comes from the data Ryan put on line, and earlier from S09. There are four regions, Peninsula, WA, Ross Sea, EA. But the satellite tags are bounded inconsistently with the ground stations, so I retagged them as shown in this map. Ground stations are shown with large dots, and the subset of AVHRR stations (described here) are small dots, both colored by region.

That done, I have computed trends for these regions and the continent using the OLS methods of TempLS V1. As with some earlier posts, I have mixed AVHRR and ground stations, with various relative weightings (including ground only).

Friday, March 18, 2011

Area weighting and a 60 stations global temperature.

I've been experimenting more with area-based weighting. I found strict Voronoi tessellation was slow for global coverage. The reason is that one really has to loop over triangles, and in R that is slow. So I investigated a variant in which a triangular mesh is used, but to get the area associated with the node, the mid-point of each side is connected to the median instead of the circumcentre. That is equivalent to apportioning the triangle area equally to its three nodes. It is also equivalent to using linear finite element basis functions.

It's much quicker, because most calculation can be done on whole vectors of nodes. I found that a 9800 node global calc took nearly a minute per month to make Voronoi meshing, compared with about 2 min for a complete lat/lon grid weighting calc. That slowness isn't prohibitive, but the modified method brings it back to a few seconds per month.

So as a test I decided to revisit the just 60 stations calculation. Looking back, I was actually surprised that that calc worked as well as it did, because the stations were not evenly distributed, and there was no effectively discriminating weighting for area. The usual grid-based weighting was used, but because the sparse stations were generally one to a cell, the only differentiation was based on grid cell area, which there was just a function of latitude.

In that earlier calc I wanted to stick to a simple objective criterion for station selection. This was that the stations be rural, have at least 90 years of data, and have been still reporting in 2009. That cut the number of stations down to 61.

Now I can use the weighting in the selection process too. It's still objective - the actual temperatures are not consulted. What I do is to take a larger set (I relaxed the record length to 50 years) and mesh that. There is a big range of weights, so I pick the top two thirds (approx), the stations covering the largest areas. Then I remesh, and select again, and then repeat (stages 170, 120, 90, 60). This gets down to 60 stations, but more evenly distributed. And they are now properly weighted by area, with no missing areas.

Wednesday, March 9, 2011

Area weighted trends for Antarctica

In an earlier thread I discussed the use of a least-squares program TempLS for calculating regional trends in Antarctica from mixed ground and satellite (AVHRR radiometer) data. It was in the context of recent papers by Eric Steig and colleagues (S09) and Ryan O'Donnell and colleagues(O10), where various statistical methods were used to achieve these objectives. In that analysis the stations were not (down-)weighted by their area density, which has the potential to excessively weight the Antarctic Peninsula warming, and also West Antarctica. I don't think S09 or O10 did this either, and it isn't obvious how to do it. But I've been developing Voronoi tessellation methods, and now I have been able to get them to work for this analysis. I'll also look at some residuals and explained sum of squares (R2).

Monday, March 7, 2011

New Blogger - Isaac Held

h/t Eli.

Isaac's blog is here.  Introduction here.  It looks very promising, and I wish him well. He's on the blogroll here, with RSS connection, so you'll see when there is a new post.

While I'm talking about bloggers, I've been finding the weekly index of news compiled by H.E Taylor, and hosted by Coby Beck, very useful. I made an index of posts, which I'll put below the jump. It goes into the future - H.E Taylor doesn't foretell, but if Coby keeps up his systematic numbering system, the posts should be there when the date arrives. There is a copy of the index on the document store.

Sunday, March 6, 2011

Voronoi weighting of temperature stations

Last post, we had a Voronoi weighting, but restricted by being able to only to convex shapes, and with fine scale weighting variations that seemed arbitrary, and likely to increase variance. Those have now been fixed, and I think area weighting can go ahead.

The revised code is here (in, and is more structured.

Friday, March 4, 2011

Spatial weighting and Voronoi tessellation.

I've been spending far too much time lately on triangular meshes and Voronoi tessellation. The idea for combining temperature measurements is to weight stations according to the area they could be said to "represent". This is important in Antarctica, for example, where the distribution is very uneven, and the Peninsula is likely to dominate just because of the number of stations.

For the global indices, this is handled indirectly through gridding. A cell average is calculated, and cells combined, usually weighted by their easily calculated area.. There's a problem, usually unresolved, with empty cells. And this is one reason why I'm attracted to Voronoi tessellation. Every cell has a reading, provided you are prepared to create a new tessellation every month.

I experimented with a simpler variable cell scheme last year. But there were still empty cells.

A Voronoi tessellation just divides the space up according to the nearest nodes. The patch around node A consists of points nearer to A than any other node. There are still some practical problems in getting it into TempLS, because the tessellation would have to be done for every month. So it has to be reasonably efficient, but more importantly, done without any intervention to fix things. In this post I'll describe my experience so far.