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.



The sinusoid had unit amplitude. On each site I convolved a unit uniform variate (default rnorm) with the function 3*exp(-0.2*i). In choosing these numbers I was trying to best illustrate the effect. I did this 383 times, assigning to each sequence the starting point of one of the Briffa sites (in sequence). The first run produced this plot:

The site starting years are shown in red, and the true sinusoid in green. The common signal, the sinusoid, is coming through fairly well above 1600, but gets ragged earlier on from there. But the raggedness is random - doing it again gives

and again

So that is what is bad about drawing plots when the data runs low. You produce spurious deviations which are not reproducible in artificial data, and are misleading. That's why Briffa didn't do it, as most scientists wouldn't.

Update. RuhRoh in comments has suggested that I should have shown smoothed noise comparable to that in the Briffa/Steve plots. I think that is a good idea. I originally used AR1 noise as a well-recognised kind of red noise, but I'll now use the same filter that Steve did. The formula is
 e1=truncated.gauss.weights(50)*20
The 20 is the factor that regulates the size of the noise relative to the sinusoid. Here are the corresponding 3 plots (random repeats):







Appendix.
Here is my R code. yy is the vector of starting years (minus 1400) (see prev post) and u1 is the ordered version (I could have used u1-1400 instead  of yy). Here is the code:
# Add AR1 noise to a sinusoid
i1=1:550; h=pi/40; s1=sin(i1*h);
 e1=3*exp(-i1*0.2);
 for (k in 1:3){   # Do it 3 times
   M=z=i1*0;   N=383;

   for(i in 1:N){   # Looping over sites
    z1=rnorm(600);
    z2=filter(z1,e1,method="conv",sides=1,circular=T)
# z2 is the AR1 noise
    j1=max(yy[i],1):550;
    K=M[j1]=M[j1]+1;
    z[j1]=z[j1]+(z2[j1]-z[j1])/K;  # Averaging
  }
  jpeg(paste(c("sinew",k,".jpg"),collapse=""));
  plot(1400+i1,z+s1,type="l");    # Output
  lines(u1,1:383/100-2,col="red");    # Site numbers
  lines(c(1400,1950),c(-2,-2),col="red");
  lines(1400+i1,s1,col="green");   #  Sinusoid
  dev.off();
}



























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.



The reason for plotting the JGR paper as well is that Steve M pushed the claim that the JGR paper had plotted down to 1400, while the Science paper stopped at 1550. They referenced different data sets, so I wanted to check if they perhaps had more data in the later paper. As you'll see, the answer is no, but they did claim in the second paper to have an improved method, which might well allow them to go further back in time.

Anyway, here is the plot of number of sites vs time, over all time.

As you can see, the number of sites is dropping rapidly before 1600, and is down to about 40 near 1550. Here is the expanded region between 1400 and 1600:

As you can see, the rate of decrease is quite sharp near 1550. There's no absolute rule on where you have to say that a plot has to be stopped. The noise rises relative to the signal in a continuous way, and I don't curently know how to quantify whether 40 sites is likely to be sufficient. But neither do the critics. What is clear is that the observed rapid changes observed in McIntyre's graph are closely associated with the steep reduction in data. In those circumstances, I would be very uncomfortable about presenting them as real. And I don't think referees would let me.

Here is a list of the first 61 sites from the 1998 Nature paper. I won't list the corresponding JGR sites, as it is very similar.
"ID""Name""Abbr""Genus""Lon""Lat""Start""End"
"276"" 862A "" Polar-Ural (historisch) "" POU_LA "" LASI "65.6366.879141990
"253"" 928A "" Mangazeja (hist. + rez.) ""MANGAZPC "" PCOB "82.366.6812461969
"252"" 928D "" Mangazeja (historisch) ""MANGAZLA "" LASI "82.366.6813071698
"301"" 913A "" Zhaschiviersk ""ZHASHIST "" LASI "142.6267.4513111708
"351"" 1041A "" Majakit river/village ""MAJAKILA "" LADA "151.9561.2813381994
"278"" 930A "" Seimchan-river ""SEIMCHLA "" LADA "151.7263.5213621991
"215"" 745A "" Sylvan Pass bei Cody "" SYLVAN "" PCEN "-110.1344.3713881983
"264"" 1017A "" Nuleger river ""NULEGELA "" LADA "127.4371.2213911990
"174"" 691A "" Jasper "" SUNW "" PCEN "-11752.2514001983
"383"" 9999Z "" Tornetrask "" TORNXX "" PISY "19.568.314001980
"206"" 731A "" Medicine Bow Peak "" MEDI "" PCEN "-107.741.314011983
"372"" 992A "" Qamdo "" QAMDOFI "" PCBA "96.9531.0814061994
"337"" 1014A "" Bilibina,Mali Anuj river ""BILIBILA "" LADA "167.6767.4714321991
"382"" 1060A ""Tschokurdach,Ochotingna r ""TSCHOKLA "" LADA "148.0570.2814341990
"48"" 630A "" Sierra da Crispo "" CRISPO "" PILE "16.2339.914411980
"332"" 1032A "" Andryuschkine ""ANDRYULA "" LADA "145.7769.2814491991
"268"" 899A "" Olenjok ""OLENJOLA "" LAGM "112.5368.5214501990
"214"" 686A "" Snow Bowl, San.Fr. Peak "" SNOW "" PCEN "-110.235.4314531983
"377"" 1063B "" Adycha river (flach) ""ADYBOALA "" LADA "137.2565.814811991
"304"" 865A "" Olenok-River ""ZWANZGLA "" LAGM "119.168.5214821990
"316"" 540A "" Lofoten, Loedingen "" LOFOTN "" PISY "16.0368.4814851978
"333"" 955A "" Balshoia Anui ""BALANULA "" LADA "165.4266.2214921991
"198"" 692A "" Highland Fire Outlook "" HIGHL "" PCEN "-112.5345.7514961983
"209"" 733A "" Powder River Pass "" POWDER "" PCEN "-107.0544.1514961983
"352"" 981A "" Omoloya river ""OMOLOYLA "" LADA "132.9870.9514961991
"380"" 1070A "" Balygichan-river ""BALYRILA "" LADA "154.5862.1314971994
"188"" 762A "" Barlow Pass, am Mt.Hood "" BARLOW "" PSME "-121.6545.3215041983
"200"" 746A "" Granite Pass, Hunt Mtn. "" HUNT "" PCEN "-107.8744.7815081983
"101"" 689A ""Hidden Peak, Wasatch Mtn. "" ALTA "" PCEN "-111.6340.5715111983
"184"" 702A "" Yosemite Park, E Eingang "" YOSE "" PICO "-119.2537.815131983
"87"" 629A "" Col de Sorba, Mt.Renoso "" SORBA "" PINI "9.242.0715181980
"202"" 760A "" Lassen National-Park "" LASSEN "" TSME "-121.5240.4515251983
"196"" 723A "" Galena Pass, Sawtooth NF "" GALENA "" PCEN "-114.7243.8715301983
"208"" 756A "" Pike Peaks "" PIKES "" PCEN "-105.0339.3315301983
"341"" 974A "" Cherskij ""CHERSKLA "" LADA "163.0568.815381991
"220"" 961A "" Ayakli river ""AYAKLILA "" LAGM "97.5369.5315401990
"192"" 755A "" Electric Lake "" ELEC "" PCEN "-111.3339.5815421983
"345"" 1047A "" Julietta north ""JULINOLA "" LADA "153.9761.1715471994
"204"" 732A "" Lone Lake "" LONE "" PICO "-118.4737.1715481983
"336"" 973A "" Batagay,Chandon-river ""BATAGALA "" LADA "138.1770.2515501991
"221"" 924A "" Ayandina-River ""AYANDYLA "" LADA "143.1768.4215531991
"125"" 704A "" Denali National-Park "" DENALI "" PCGL "-149.5863.6715541983
"210"" 700A "" Sierra Blanca, (Ruidoso) "" RUID "" PSME "-105.6733.3315541983
"339"" 956A "" Bulun river ""BULUNLAD "" LADA "154.9365.115541991
"353"" 958A "" Rossocha river ""ROSSOCLA "" LADA "149.4365.1515551991
"187"" 726A "" Baldy Peak "" BALDY "" PSME "-109.5533.9715561983
"348"" 969A "" Srednie-Kolymsk ""KOLYMSLA "" LADA "153.767.2515561991
"381"" 1061A "" Kubaka right bank middle ""KUBA2LAG "" LAGM "159.9563.6715601994
"240"" 864A "" Kotuykan-River ""KOTUYALA "" LAGM "104.2570.5815631990
"191"" 748A "" Crater Lake, NE-Medford "" CRATER "" TSME "-122.1742.9715641983
"246"" 926A "" Kuonamka-River (trocken) ""KUONLATR "" LAGM "112.8269.9315641990
"302"" 906A "" Zhigansk ""ZHIGANPI "" PISY "122.3366.5215641991
"354"" 970A "" Sartan river ""SARTANLA "" LADA "132.9364.9315641991
"190"" 683A "" Cottonwood Pass "" COTTON "" PCEN "-107.5838.6715651982
"356"" 980A "" Tirekhtjakh river ""TIREKHLA "" LADA "137.4767.6215651991
"203"" 695A "" Mt. Lemon "" LEMON "" PSME "-110.7832.4515681983
"234"" 927A "" Khandiga-River ""KHANDILA "" LADA "137.7562.4715681990
"238"" 914A "" Khotugn-Uladan-Tukulan ""KHOTUGPI "" PISY "125.863.3815681991
"228"" 1011A "" Balschaya Kamenka river ""KAMENKLA "" LAGM "93.9771.3215691990
"243"" 871A "" Kulyumbe River ""KULYUMLA "" LASI "89.0768.0515741990


















Briffa et al Nature 1998: ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/n_hem_temp/nhemtemp_data.txt

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).



Continent ground station trends

These are difficulties with the sparse samples in some tof the regions, most notably with West Antarctica. So the whole continent is simplest. I'll start with unweighted (by area):

Stations reporting in each year

A map of ground stations colored by region

The trend is quite high at 0.175°C/decade. Now I'll add weighting by area, as described here. The map shows for a specific month, Dec 1987,  the triangular mesh created, the subdivision used to weight each station, and the circles show, by area, the size of the weights. I have not smoothed, nor tried to cut out the Weddell Sea.




The trend is lower, at 0.114°C/decade. This is in line with earlier observations that area weighting reduces the trend, probably by boosting the effect of the less rapidly warming East Antarctica stations, which occupy the largest region.

Adding AVHRR readings


As with previous posts, I've introduced the satellite AVHRR data, which comes as grid data, as artificial stations. I've selected only 1 in 4, so there are 1392 of them. A factor is needed to determine the relative weighting for combining them with ground stations. I now use a range from 0 (all satellite) to 1 (all ground), with 0.5 indicating that total weights for satellite and ground are equal. That is for all years, so that before 1980, there is ground data only, and the factor doesn't matter. Because of the way it is totalled, this means sat data gets a higher weighting (at 0.5) after 1980. This is a small effect.

So here are results from satellite data only. Like O10, I'm using Eric Steig's archived AVHRR data, which goes up to 2006 only.

As expected, the trend is quite high. One of the reasons for interest in S09 was that it looked into the discrepancy between the ground trend (lowish) and AVHRR (high).

So now looking at the combination of satellite and ground, remembering that before 1980 it is ground only:

So the trend is lower - closer to ground than AVHRR, reflecting the shorter period of AVHRR.

Antarctic Peninsula trends

So now I look at regional trends. I'll start with the Peninsula, which does not have problems of years without data. Here are ground stations only:

AS expected, the trend is high. Now if the AVHRR stations are included, with equal (total) weighting:

East Antarctica trends

Again firstly the ground stations:


As reflected in work in previous posts (and S09, O10) ground stations in East Antarctica have been warming very slowly. The inclusion of AVHRR results at nominally equal weighting increases this by a fairly small amount:




West Antarctica and Ross Sea

I found a specific difficulty with WA, which I think Ryan and Eric have mentioned. In about 1980, about the time AVHRR satellite information is available, there was also a change to measurement (AWS). Only one station in WA was reporting before then, Byrd, and there was a gap of four years with no readings at all. Byrd changed to AWS, and this is listed as a new station. So there is no way of bridging the gap without making some assumption about how the old Byrd relates to Byrd AWS. And even that would be slender information.

A first option for WA is to combine it with the Ross Sea stations, where there is a fuzzy boundary anyway. There are just enough stations pre-1980 to get a result:

Again, as expected, the trend is high. Adding AVHRR stations does not change this.



West Antarctica trends

The other option for WA is to look at results after 1980 only. That gives:

The trend is very high indeed, but so is the uncertainty, and the stations are few.

It seems the AVHRR data stabilizes the average when it is available, so just checking those years:


This did reduce the trend slightly, and improve the error range.








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.



Here is the plot of the initial selection and mesh, from 4 perspectives. The circles show (by area) the weights attached:

If it looks ragged around the rim, that's because I plotted only triangles where all nodes could be seen - otherwise there are end-of-world effects. The faint lines are the underlying mesh, and the dark lines show the area assigned to each station in the weighting.

Here are the stages of concentration, seen from N America.

And here is the final stage, from other perspectives:

and a lat/lon plot:



For comparison, here is the original 61-station selection:


Finally, here is the time series of numbers of these stations reporting in each year:


Results

Firstly. for comparison, here were the results of the original study compared with CRUTEM3 (land only. The observation was that the same general trend was followed, with of course more noise in the 61-station sample.

Now here is the new plot compared with CRUTEM3 and with GISS Land/Ocean (all to base period 1961-90): (Update - corrected following Anon comment below - the graph I originally showed is here)



It's a lot less noisy. It is also rather closer, in both trend and yearly, to GISS Land/Ocean rather than to CRUTEM3. This is interesting, because although it has 60 land stations only, they are weighted for total global coverage, and with islands heavily represented. This suggests that the difference between the two data sets may reflect weighting rather than the kind of measurement.

Next

I think the smoothing achieved with refinement by weighting worked rather well. I'll look at rationalising larger subsets that way.

I also want to look at schemes for avoiding monthly meshing. Creating these rather larger cells around a few stations means that the stations could be weighted not by individual areas but by their share in larger areas. Since several stations would inhabit such cells, frequent remeshing would be unnecessary - only required when there were no stations at all in a cell.