Sargasso Sea dO18

Oct 14, 2007 – 10:10 PM

Recently Keigwin’s Sargasso Sea dO18 temperature reconstructions have been mentioned in the climate public eye. Keigwin’s reconstruction famously has a warm MWP as shown below. This reconstruction uses modern dO18 measurements at Station S to calibrate two cores (with modern Station S values specifically shown in the graphic below):

I observed last year that the Juckesian version of the Moberg reconstruction had a very different appearance using Sargasso Sea instead of Arabian Sea G Bulloides and using Indigirka instead of Yamal – a point that I re-iterated recently in light of the final Juckes paper. I discussed Juckes and the Sargasso Sea last year and in my Review Comments at CPD. We’ve discussed the Juckesian pretext for withholding the results with Indigirka. His pretext for excluding the Sargasso Sea SST reconstruction from dO18 is that the data ends in 1925 (a point was recently picked up at realclimate).

The Sargasso Sea series used by MSH2005 is sometimes mistakenly presented as having a 1975 data point representing the 1950–2000 mean, but the actual end point is 1925 and so this series is also omitted.

Today I want to discuss exactly what is and isn’t in the Keigwin Sargasso Sea data, amending some first thoughts based on clarifications provided below by several readers. Keigwin (1996) (see here) archived a data set with SST, salinity and dO!8 values from 1955 to 1995 here . Keigwin 1996 says:

As expected, the effects of decreased SST and increased salinity during the late 1960s combined to increase the dO18 value of equilibrium calcite. Over the full 42-year series, linear regressions between the SST, salinity and dO18 values show that temperature accounts for about two-thirds of the isotopic signal ( r2 = 0.61), whereas salinity accounts for one-third (r2 = 0.29). Thus, by comparison to sea surface changes during the past several decades, it is reasonable to interpret the foraminiferal isotopic data mostly in terms of SST change.

This statement, combined with the data set, looks like a calibration set. The above discussion – in which variance is allocated to factors – looks very much like a conventional statistical analysis. As I observed in my first cut at this:

The statistical model from the two-factor model is so uncannily accurate (r2= 0.9988) that I’m wondering whether some of the data here has been forced somehow. Actual and fitted dO18 values are shown below: you have to look very closely to see any discrepancy. The dashed line is not simply a replot of the data: you can notice some very small differences. I’ve never seen any proxy model that is remotely this accurate – so it is puzzling. Some time, I’ll try to check whether (say) the salinity measurements are forced and that’s why the model is so accurate.

My surmise that the model was forced was correct, but I missed what was being forced. The explanation is latent in the article – it’s just that what was being done was something that I hadn’t expected. Keigwin forced the dO18 values from SST and salinity measurements at Station S and then did a regression using the forced data set.

In order to gauge the influence of the annual variability of SST and salinity on the oxygen isotope ratios that might be recorded by surface-dwelling foraminifera at Station “S”, I calculated the 6% value of calcite precipitated in oxygen isotopic equilibrium with seawater (Fig. 3) (26).

What exactly does this calculation do – I guess that it somewhat quantifies the expected contributions of salinity and SST to dO18 variation, given the measured variability of SST and dO18. Here’s a way of representing the calculation which clarifies things for me: the black points show reconstructed SST from core sediments against the dO18 of the core sediments (in a linear relationship according to the reconstruction transfer model), while the red points show two 20-year averages of reconstructed (“modeled”) dO18 from measured SST and salinity using the Deuser/Keigwin methodology. Keigwin observed that core dO18 values (which integrated decadal periods) were outside the range of modern (estimated) values from Station S using current transfer models. His plot shown above amalgamates the modern and core data along these lines.

sargas22.gif

While I was looking at this data, I became intrigued with some of the properties of the original calibration of dO18 to SST in the Station S data set, which illustrate rather nicely some of the issues involved in reconstructions of any kind. The dO18 values here are forced and that’s why the model is 100% accurate. But let’s suppose that we actually did have a 100% accurate model and see what the effect of not knowing salinity does to the error structure of SST reconstructions. These are just thoughts right now.

Two-Factor dO18 Model
Although dO18 is used to reconstruct SST, the actual physical relationship goes from SST and salinity to dO18. As an experiment, I tried a simple regression of dO18 values against SST and salinity data as follows:

fm=lm(dO18~SST+Sal_psu,data=sarg.mod)

sargas16.gif

First, let’s look at the dO18 model (note that we’re fitting dO18 here rather than SST) if we just use SST for the reconstruction, as shown below. You can readily see that the model using only SST is less successful and that the residuals from the model leaving out a factor have positive autocorrelation (and a failed Durbin-Watson statistic.) This is a rather nice demonstration of how the Durbin Watson statistic gives evidence that a salient factor has been left out of the model. In econometrics, this is evidence that the model is mis-specified – a statistical method that Juckes et al for example simply repudiated without providing any statistical authority. Juckes’ say-so was high enough authority for editor Goosse, I guess. In this case, there is a fair amount of autocorrelation in the salinity data and this carries forward into the residuals from the incompletely specified model.

sargas18.gif

The above plot illustrates the situation for the modeling of dO18. However in climate studies, the procedure of interest is the reverse: the use of dO18 as a “proxy” to “reconstruct” SST. A common strategy is the use of “inverse regression” of SST against dO18 -even though the physical relationship is the opposite. I use the term “inverse regression” here as it’s used in calibration statistics off the Island. There’s been much made of inverse regression in climate science (von Storch, Juckes), where the term “inverse regression” does not coincide with “inverse regression” as used in conventional statistics, but to describe a variant of partial least squares used on the Island in various guises- a topic that I’ve discussed here. If the calibration is done from the regression of dO18 against SST, this is called “classic calibration” in calibration studies. For the discussion here, don’t worry about the terminology in the various studies, but about what is being done. Here I’m regressing SST against dO18 as shown below:

fm.inv=lm(SST~dO18,data=sarg.mod)

The figure below shows the SST reconstruction (top) and residuals (bottom) using inverse regression based only on SST – a defective model to be sure, but a model containing some information. I’ve also shown the reconstruction using variance matching. In a univariate case, this simply dilates the reconstruction a little, in which the price for matching the variance is the increase in the average size of the residuals. As you can readily see, the autocorrelation in the reconstruction residuals is intimately related to the defective model (in which the autocorrelation of the salinity is imposed onto the residuals.) You can also see that the variance-matching reconstruction and the inverse regression reconstruction are intimately related, with the one just being a dilation of the other.

sargas19.gif

Thus, estimating the confidence interval of the reconstruction depends on the size of the impact of unreconstructed salinity. Because the mis-specified factor has considerable autocorrelation, the confidence intervals cannot be calculated in an elementary way. A common recipe in climate studies is to estimate the reduction in number of degrees of freedom. While the approach is not absurd, the devil is always in the details. I think that the salient question is: how many measurements are needed in an average so that the probability of the average error being within \epsilon of 0 corresponds to a “standard” value using i.i.d. normal errors. Some time, I’ll re-examine the typical recipes on the Island.

It should also be a little sobering to see how large the errors are in this case when one factor is unavailable for a proxy where, in this case, due to forcing, the model is specified to 100% accuracy. The standard deviation of the residuals was 0.2, while the standard deviation of the SST was 0.29 – this is in a perfect model with unknown salinity. If the dO18 model explained (say) 50% of the variance in the model, then these errors would have to be added in. I’ll do the calculations some time to try to figure out where the model “breaks even”.

By Stephen McIntyre | Posted in General | Tagged , | Comments (41)

USHCN #3

Oct 14, 2007 – 4:40 PM
By Stephen McIntyre | Posted in Surface Record | Comments (285)

Unthreaded #22

Oct 14, 2007 – 4:35 PM
By Stephen McIntyre | Posted in General | Comments (559)

Almagre Vistas

Oct 13, 2007 – 8:54 PM

I hope that many of you have visited Pete Holzmann’s photo gallery here. If you click on each picture on this page, you get a separate “gallery”. Today, I’m going to post up a quick overview together with a few maps. Continue reading →

By Stephen McIntyre | Posted in Almagre | Comments (85)

A Little Secret

Oct 12, 2007 – 5:44 AM

Don’t you think that someone on the Team might have been a little curious as to what bristlecone ring widths have done during the past 25 years? For this, we have the classic excuse of Michael Mann and the Team for not updating bristlecone and proxy records is that it’s not practical within the limited climate budgets:

While paleoclimatologists are attempting to update many important proxy records to the present, this is a costly, and labor-intensive activity, often requiring expensive field campaigns that involve traveling with heavy equipment to difficult-to-reach locations (such as high-elevation or remote polar sites). For historical reasons, many of the important records were obtained in the 1970s and 1980s and have yet to be updated.

From the first moment that I got involved with paleoclimate, it seemed obvious to me (as it is to anyone not on the Team) that, if the classic “proxies” are any good and not merely opportunistic correlations, that there is an ideal opportunity to perform out-of-sample testing of the canonical Team reconstructions by bringing the proxies up-to-date. I wrote an Op Ed in February 2005 for the National Post entitled “Bring the Proxies Up to Date”, where I expressed the view that this was really the first order of business in Team world. While the addition of new proxies is also important and nice, this is not the same thing as out-of-sample testing of the proxies used in MBH99, Crowley and Lowery etc – especially the bristlecones.

I’ve continued to satirize this failure pointing out that several of Graybill’s classic bristlecone sites were easily accessible from UCAR world headquarters in Boulder and that no heroic expedition was required to update, for example, the Graybill sites to the west of Colorado Springs.

To get to these sites from UCAR headquarters in Boulder, a scientist would not merely have to go 15 miles SW of Colorado Springs and go at least several miles along a road where they would have to be on guard for hikers and beware of scenic views, they would, in addition, have to go all the way from Boulder to Colorado Springs. While lattes would doubtless be available to UCAR scientists in Colorado Springs, special arrangements would be required for latte service at Frosty Park, though perhaps a local outfitting company would be equal to the challenge. Clearly updating these proxies is only for the brave of heart and would require a massive expansion of present paleoclimate budgets. No wonder paleoclimate scientists have been unable to update these records since Graybill’s heroic expedition in 1983.

Pete Holzmann (Mr Pete), who lives in Colorado Springs, agreed with this satire and this led to what I’ll call the Starbucks Hypothesis: could a climate scientist have a Starbucks in the morning, collect tree rings through the day and still be home for dinner?

To make a long story short, last summer, when my wife and I visited my sister in Colorado Springs and I thought that it would be rather fun to test the Starbucks Hypothesis and I gave a bit of a teaser report in late July, promising some further reports in a few weeks, but I got distracted by the Hansen stuff. At the time, I mentioned that, together with CA reader Pete Holzmann and his wife Leslie, we visited some bristlecones in the Mt Almagre area west of Colorado Springs.

But I have a little secret which I’ll share with you as long as you promise not to tell anyone: our objective was to locate the precise site sampled by Graybill. Not just that. Prior to the trip, I obtained a permit from the U.S. Forest Service to take dendrochronological samples from bristlecones on Mount Almagre and we did more than look at pretty views; we obtained up-to-date bristlecone samples. I only went up Almagre on the first day. Our permit lasted a month and Pete and Leslie spent two more days on Almagre, finally locating and sampling tagged Graybill trees on the third day.

Altogether (and primarily through the efforts of Pete and Leslie), our project collected 64 cores from 36 different trees at 5 different locations on Mount Almagre. 17 Graybill trees were identified, of which 8 were re-sampled. All the cores are currently at a dendrochronological laboratory, where sample preparation and scanning steps have been completed. Cross-dating is now taking place. For the most part, we tried to sample non-stripbark trees in keeping with NRC recommendations, but some stripbarks were sampled to reconcile with Graybill. Of the tagged Graybill trees, the tag numbers do not reconcile with the archive identification numbers and, in the absence of any concordance, reconciliation may prove more difficult than one would think.

We will archive at WDCP detailed information on the location of all samples (current spreadsheet is here) , which has already been sent to the U.S. Forest Service. Photographs of each tree are shown gallery here. Here’s a fun presentation that Pete prepared of our Day 1 itinerary. Here is a Google Earth tour. If you run it and when Google Earth comes up, go to Tools “ Play Tour , you’ll have some fun.

Some expenses have been incurred for this expedition. Leaving aside travel expenses (which were vacation expenses that I was going to incur anyway), the jeep got a bad scratch on the first day and cost about $500 to repair plus some more repairs from Days 2 and 3; it’s going to cost a few thousand at the dendrochronology laboratory for the sample prep, scanning and cross-dating as this has been done on a contract basis. I’ve submitted an abstract to Rob Wilson’s divergence session at AGU and would like to present these results (and to cover Pete’s expenses if he can come). On the basis that we submit a data paper to a journal, publication expenses will be another $800-1000 or so (academic authors pay the journals to publish). This has been a Climate Audit project so I’d like readers who contribute to the top jar to think about a special contribution for the bristlecone sampling. Maybe Martin Juckes, James Hansen and Michael Mann will contribute as well – I’m sure that they are all anxious for the results.

I’ll add some more information later in the day. Right now I’m off to visit the dendro lab and see how things are coming along. In 2002, Malcolm Hughes sampled bristlecones at Sheep Mountain and nothing has been reported or archived from this study. In 2003, Lonnie Thompson sampled ice cores at Bona-Churchill and we’ve heard nothing about it. One might guess that 20th century dO18 levels were not high as, at the nearby site of Mount Logan, 20th century dO18 levels were lower than earlier levels, attributed to regional changes in circulation rather than temperature.

I’ve obviously been very critical of what appears to be opportunistic reporting of results. With my experience in mining speculations, I fully understand how much temptation that there is to delay reporting of “bad” results in the hope that later drill holes in the program will salvage things. But you don’t have any choice in the matter – you’re obliged to report the results. Plus investors are smart enough to now that delayed results are virtually never good results.

Right now I have no idea what the sampling will show – maybe it will show a tremendous response by the bristlecones in the past 20 years – perhaps due to CO2, nitrate or phosphate fertilization, perhaps due to temperature increases. Maybe they won’t go up and we’ll hear more about the divergence problem. I don’t expect these particular measurements to settle anything. But jeez, doncha think that someone would have tried to find out?

Anyway I promise one thing: the measurements are going to be made public as soon as I get them. Just like a mining project. No waiting for 5 or 10 or 25 years like certain people. No losing the data like other people. Whatever they show. As soon as I get the cross-dated measurement data, we will immediately send it to the World Data Center for Paleoclimatology (which I expect to take place within a few weeks.) I hope that this will set an example to the trade as to the type of turn-time which is practical.

I’m visiting the dendro lab today to say hello and I wanted this to be on the record before my visit. I’m not sure how far along they are, but I think that they’ve finished sample prep and scanning and have started cross-dating. I don’t expect any results, but they may be far along enough so that I’ll have an impression of what the result growth will have been. So I wanted to be on the record on the planned schedule prior to my knowing anything about the results, just in case I get an impression today of what the recent growth has been.

UPDATE 2 pm. OK, I’m back from the dendro lab in Guelph. They are further along than I expected. The longest core is 883 years (Tree 30A). This had a Graybill tag 84-55, but if you go to the archived measurement data
 and look for ALM55 (which would presumably be the match), there are no corresponding measurements; there is obviously a sequence ALM01, ….; there is an ALM53 and an ALM60. Is there an alter ego somewhere or is it missing? Right now we don’t know.

After sanding the core is scanned. The measurement of ring widths is semi-automatic. For these bristlecones, earlywood and latewood was easy to distinguish. Using a magnified version of the scan, each ring is picked out (with the computer recording the pick). The computer then yields back the measurements. I’ve posted up a couple of print-screens showing the most recent widths for 30A and the widths in the mid-19th century.
.

Below is a print screen showing the 30A ring widths from 1124 to 2007. I’ll post up a re-plot at some point with a legible x-axis. For orientation in the absence of a scale, the upspike on the left is 1174; there are low values from 1353 through the early 1400s; there is a 1690 spike; 1865 and 1880 are upspikes; 1941 is a small upspike.

According to the Team hypothesis of a positive linear relationship between temperature and ring widths, the warm 1990s and 2000s should have yielded the widest ring widths in history. What do you think? This is only one tree, but my quick impression was that recent growth was not elevated. So this looks like a Divergence “Problem”. If CO2 or other fertilization has been a factor, then I hate to think what the growth would have been without the fertilization. Remember the NAS panel saying that the Divergence Problem only affected high-latitude sites? Maybe they should have done some testing before they opined on this.

BTW while I’m critical of how the Team uses dendro information, I think that it is well worth supporting the collection of dendro information, even if it’s meaning is not clear right now. It has the advantage of being well-dated – and when you see the problems with dating ocean sediments, it’s nice to have some records that are well dated. There’s a role for it; so please – no posts dumping on dendrochronology. The dendrochronologists who’ve been doing this work (and who I will credit in due course) are excellent people.
almagr8.jpg

Tree 30A Ring widths. 2007 on right. i’ll replot this some time soon. For orientation in the absence of a scale, the start is 1124; the upspike on the left is 1174; there are low values from 1353 through the early 1400s; there is a 1690 spike; 1865 and 1880 are upspikes; 1941 is a small upspike.

By Stephen McIntyre | Posted in Almagre, bristlecones, Post-1980 Proxies, Proxies, Replication | Comments (377)

Juckes and the West Greenland Stack

Oct 10, 2007 – 6:32 PM

Last fall, I discussed information sources on West Greenland ice core series, noting that the West Greenland version attributed by Juckes to Jones et al 1998 was a version that I’d not seen before. While I was looking at the proxy decisions in Juckes et al, I noticed the following intriguing rationalization:

The Greenland stack data used by MBH1999 is a composite of data analysed by Fisher et al. (1996), but the precise nature of the composite is not described by Fisher et al. (1996). … The Crete ice core series [of Jones et al 1998] is preferred to the “stack” series used by MBH1999 because it is better documented.

So one asks: exactly what is the “better documentation” preferred by Juckes and the Euro Team? Or is this “better documentation” a figment of their imaginations?

Jones et al 1998 and MBH98-99 both cite Fisher et al 1996 as a reference. As it happens, Fisher et al 1996 does not appear to describe either series (it only illustrates rather short versions, although it describes Fisher’s stacking methodology.) So any documentation advantage in this reference is a draw.

Fisher (2002) Table 1 reports the use of a version identified as DELNORM6.CWG (with that exact nomenclature) – “CWG” standing for Central West Greenland.: the file itself was ASCII. A couple of years ago, Fisher sent me a CD with a variety of series on it. Experimenting a little, I determined that the series DELNORM6.CWG, included on the CD, matched the MBH99 version 100%. The DELNORM6.CWG file listed its components for each step e.g. the following components in one of the steps. So any person who received the DELNORM6.CWG version from Fisher would have received a relatively adequate recipe for how it was constructed. In Mann’s MBH98 data archive, the email transmitting this series to Mann (from Frank Keimig of the University of Massachusetts) is preserved for some reason.

A CT85A-1Y.CRT 1983 AD 362 212 144
B CT84B-1Y.CRT 1982 267 211 50
D CT84D-1Y.CRT 1982 217 211 0
CRETE CT74-1YS.CRT 1973 922 202 714
MILCT MC73-1Y.MIL 1966 791 195 590
GISP2 GISP1YR.GSP 1985 718 214 498 DEUTERIUM
GRIP SIG#1YR1.GRIP 1979 917 208 703
GRIP SIG#2YR1.GRP 1979 917 208 703
GRIP SIG#3YR1.GRP 1979 917 208 703

I’ve attempted to locate the Jones et al version within the many variations in Fisher’s files, but so far have been unable to locate an exact reference. There are multiple series from Crete, including many of the series in the Fisher stack described above (Crete, A, B, C, D, H). The Jones version is not the same as any of the potential components. The Jones version goes from 1000 to 1983; the only version going back to 1000 (actually to 553) ends in 1973; there are a number of potential components ending in 1983, but all start after 1000. So the Jones version appears with a high degree of probability to also be a CWG stack – except in this case, unlike DELNORM6.CWG, we don’t know the precise components of the stack, although the Jones version is obviously closely related to these components.

So back to the question: is the “better documentation” of the Jones version a figment of the Euro Team’s imagination – the full team including Briffa, Osborn, Moberg, Allen, Hegerl, Weber, Esper. It sure seems that way to me. Perhaps Juckes can identify some such documentation, but he didn’t do so in the article itself. I guess Goosse and the referees didn’t ask him about it.

In addition, given Juckes’ professed concern about ensuring that all series versions have properly received permissions from their originators, it’s interesting to contemplate what permission (if any) Juckes has for archiving the Jones et al version of the Fisher data (or for that matter the MBH version.) I am unaware of any public archiving of either version by Fisher himself.

I presume that Frank Keimig got the DELNORM6.CWG version from Fisher either directly or indirectly. In the email at Mann’s FTP site, Keimig mentions re-formatting the series – in which process the documentation in the Fisher version was removed. Mann got this version from Keimig; Juckes presumably got it from Mann, either directly or from Mann’s archive. Did anyone in this chain ever get permission from Fisher to archive the data publicly? I doubt it. I don’t suppose that Fisher has any objections to people using his data provided they acknowledge him – why would he? But within Moberg’s definition, Fisher himself, to my knowledge, never publicly archived the version now disseminated by Juckes.

I believe that the Jones version is an earlier version of the DELNORM6.CWG series (or some close variation.) Jones may have got his version directly from Fisher who’s very cooperative (or he could have got it from Ed Cook or someone like that indirectly). Presumably Juckes got the Jones version from Briffa or Jones. (BTW Jones refused to provide the versions that he used to me when I requested.) Moberg expressed concern about the possibility of Sidorova issuing a new version and that he did not want to preempt any adjustments by Sidorova to her series. It seems inescapable that the DELNOM6.CWG version (which is used in Fisher 2002) is a version that Fisher prefers to the version used in Jones et al 1998. Yet Juckes seems to have no compunction about not only archiving the Fisher version used in Jones et al 1998 version – a version that never previously had been posted on the internet digitally – and in using a version seemingly not preferred by the originator of the data, all using the pretext of “better documentation”, though no evidence of such “better documentation” was provided anywhere.

References:
David A. Fisher, 2002. High-resolution multiproxy climatic records from ice cores, tree-rings, corals and documentary sources using eigenvector techniques and maps: assessment of recovered signal and errors, The Holocene 12,4 (2002) pp. 401-419
Fisher, D. A., Koerner, R. M., Kuivinen, K., Clausen, H. B., Johnsen, S. J., Steffensen, J. P., Gundestrup, N., and Hammer, C. U.: Inter-comparison of ice core and precipitation records from sites in Canada and Greenland over the last 3500 years and over the last few centuries in detail using EOF techniques, NATO ASI Ser, Ser. I, vol. 41, edited by: Jones, P. D., Bradley, R. S., and Jouzel, J., pp. 297–328, Springer-Verlag, New York, 1996.

By Stephen McIntyre | Posted in Juckes et al 2006 | Tagged , , | Comments (44)

Juckes and "Restricted" Data

Oct 7, 2007 – 8:47 AM

Many climateaudit readers will remember Mann’s “CENSORED” directory, in which Mann calculated principal components on a network that excluded bristlecone pines (which needless to say didn’t have a HS shape. Now Juckes et al introduces us to a new type of climate data: “restricted” data. The Team has introduced a novel data classification system – PG and R. Juckes et al say that the Indigirka series is R-rated and so it can’t be used in their reconstruction. Yes, R-rated tree ring data. Data so salacious that you have to keep it under lock and key.

Is it only under-18s that are not allowed to see R-rated tree ring data? Can we show it here if Kristen Byrnes promises not to look?
Or are all climateaudit readers prohibited? Is this a bit like pornography that is only available to priests? You think that I am juck-ing? Here are their exact words from the Euro Hockey Team for excluding the Indigirka series:

The Indigirka series used by MSH2005 is not used here because it is not available for unrestricted use.

I wonder what went through the minds of editor Goosse when he read that this was a “restricted use” proxy? Did Goosse ask what the restrictions were? Or referee Gerd Bürger or the other two anonymous referees? Since they don’t appear to have asked or weren’t bothered by the answer, let’s ask the question here. And, by the way, for readers who may be offended by salacious tree ring data, proceed at your own risk as the R-rated data is shown graphically in the post continuation. Yes, you too can see what the climate priests keep in their secret cupboard. If you are not over 18, please do not continue without parental guidance. Continue reading →

By Stephen McIntyre | Posted in Juckes et al 2006, Moberg [2005] | Tagged , | Comments (243)

Gridding from CRN1-2

Oct 4, 2007 – 1:26 PM

In order to assess the 2006 anomaly in John V’s 2006 calculation reported previously, I’ve calculated a gridded anomaly for the U.S. using standard spatial smoothing software in R. This enables the calculation to be done in a few lines as explained below. I’ll also show results which reconcile most of the results to date. I’ve shown the details of calculation step so that there will be no confusion. I also want to show how simple it is to produce gridded data from station data. It only takes a couple of lines of code.

CRN Ratings
First I downloaded the current ratings from Anthony Watts from http://www.surfacestations.org/downloads/USHCN_stationlist.xls and saved it as a tab-separated ASCII file. I first deleted the comment columns and changed the apostrophes to _. This is important for R programs. I then added these ratings to my USHCN info file details as follows:

url=”d:/climate/data/station/watts”
#downloaded from http://www.surfacestations.org/downloads/USHCN_stationlist.xls
#some columns deletedl also replace ‘ with _
loc=file.path(url, “stationlist.txt”)
watts=read.table(loc,sep=”\t”,fill=TRUE,nrow=1221,header=TRUE)
details$watts=NA
test=match(details$id,watts[,”USHCN.ID”]);temp=!is.na(test); sum(temp) #1221
details$watts[temp]=watts[ test[temp],5]
details$watts[(details$watts==0)&!is.na(details$watts)]=NA

I then loaded my collation of USHCN station data and extracted the stations with CRN ratings of 1 or 2 as follows:

load(“d:/climate/data/station/ushcn/ushcn.collation.tab”)
temp=(details$watts< =2)&!is.na(details$watts);sum(temp)
index_good=(1:1221)[temp]
good=ushcn.collation[[2]][,index_good] #2544 53
M=ncol(good)

Anomaly Calculation
I then converted this to an anomaly basis and calculated annual averages. I did a separate calculation for 2006 using available months only. The functions anom and ts.annavg are available in scripts/gridcell/collation.functions.txt.

good.anom=good
for(i in 1:M) good.anom[,i]=anom(good[,i])
good.ann=ts( array(NA,dim=c(nrow(good)/12,M)),start=tsp(good)[1])
for(i in 1:M){ good.ann[,i]=ts.annavg(good.anom[,i])}
temp=(c(time(good.ann))>=1890)
good.ann=ts(good.ann[temp,],start=1890)
time0=c(time(good.ann)) #1890:2006
for(i in 1:M) good.ann[N,i]=mean(good.anom[(nrow(good.anom)-11):nrow(good.anom),i],na.rm=T)

Gridding
Gridding is a type of spatial smoothing for which there are a variety of standard functions and no need for people to ponder over triangles or hexagons (which don’t matter much) or to write new code. If one wants to vary smoothing parameters (in this perspective, GISS smooths much more than CRU), then one should do so in the context of standard software where there is working experience. In this case, I used the function interp from the R-package akima (which I had previously installed.) I used three “masks” to screen out values which otherwise would get interpolated into the ocean or into Canada (tho not much turns on this.) Because each gridpoint represents an equal area, I simply calculated the mean of the values.

library(akima)
N=length(time0)
annual=rep(NA,N)
for (j in 1:N){
X=data.frame(details[index_good,c(“long”,”lat”)],good.ann[j,]) #locations of “good” stations
names(X)=c(“long”,”lat”,”anom”)
temp=!is.na(X$anom)
gridded < – interp(X$long[temp], X$lat[temp], X$anom[temp])
gridded$z=gulfmask(gridded$z)
gridded$z=atlmask(gridded$z)
gridded$z=canadamask(gridded$z)
annual[j]=mean(gridded$z,na.rm=T)
}
annual.tobs=ts(annual,start=1890)

The mask functions were not particularly elegant. The akima package with this data yields a 40×40 grid and I simply set values at inappropriate gridpoints to NA as follows.

gulfmask=function(Z){
Z[21:31,1:2]=NA;Z[21:30,3]=NA;Z[22:30,4]=NA;
Z[23:30,5]=NA;Z[c(23:24,27:30),6]=NA; Z[27:28,7]=NA; Z}

atlmask=function(Z){
Z[40,c(1:26,29:31)]=NA; Z[38:39,1:26]=NA; Z[37,1:23]=NA;
Z[36,1:16]=NA;Z[35,1:15]=NA; Z[34,1:13]=NA;Z[33,1:11]=NA;
Z}

canadamask=function(Z){
Z[40,37:40]=NA; Z[37:39,35:40]=NA; Z[36,34:40]=NA;
Z[34:35,32:40]=NA;Z[32:33,29:40]=NA; Z[31,c(28:29,33:40)]=NA;
Z[30,36:40]=NA;Z}

Here’s an illustration of the points that are masked out by these masks. I didn’t bother with masks in other places because the algorithm didn’t interpolate over there, but it could be esily done.

usgrid52.gif

Here’s an example of the type of annual map that is generated by the smoothing function. This was generated by doing a single step and plotting from the object gridded produced in the above routine as follows (using the R package fields). You’ll note that this smoothing function does not calculate gridded values at the edges which are an extrapolation from the stations. (This map also shows the stations used in the calculation.)

par(mar=c(3,3,2,1))
breaks0=seq(-3,3,.25)
image.plot (gridded$x,gridded$y,gridded$z, col=tim.colors(length(breaks0)-1),breaks=breaks0,ylab=””)
contour(gridded, add=TRUE,ylab=””)
points (X$long[temp],X$lat[temp], pch = 19,cex=.7)
US(add=TRUE,lwd=2)
title(main=paste(“USA “, time0[j]))

usgrid53.gif

So it’s pretty easy to do a gridding calculation using standard software. Hansen’s program has very ad hoc statistical methods in which stations are adjusted from other stations and gridcells calculated. At the end of the day, the best approach to this aspect of his program will probably be to see what smoothing parameter best emulates his method and thereby restore it to a known statistical methodology.

COMPARISONS

First here is a comparison of my calculation using CRN1-2 data to NOAA results from 1890-2007. (I re-iterate for the N-th time that analyses on U.S. data are not representative of the ROW.) NOAA has values to Aug 2007 and I’ve estimated the last 4 months as the same as the corresponding 2006 months for now. NOAA has run warmer by about 0.2 deg for the past decade and was about 0.1 deg C cooler in the 1930s.

usgrid54.gif

Second here is a comparison of CRN1-2 data to GISS pre-Y2K results. As you see the GISS with the Y2K error was running warm by about the same as NOAA – which may be one reason why they didn’t spot the error. GISS runs a bit warmer than NOAA in the early and mid portions of the graph due to their greater efforts to do an urban adjustment.

usgrid55.gif

Third, here is a comparison of CRN1-2 data to GISS post-Y2K results (but not including the September change from SHAP to FILNET, which seems much too opportunistic.) Relative to CRN1-2, GISS tends to warm the cool periods (1960s, 1900s) and to cool the warm periods (1930s), reducing the size of the mid-20th century cooling. For the U.S., as noted before, by using the trends from “unlit” stations to establish trends – in a context where there are a lot of decent rural stations with long records – the GISS methodology for USHCN stations seems more appropriate than the NOAA (and probably the CRU methods). Unfortunately, this method is not used in the ROW, where there is nothing equivalent to the USHCN.

usgrid56.gif

Finally, here is a reconciliation to John V’s annual results sent ot me yesterday. His 2006 results run much colder than the ones that I calculate- which is presumably an artifact of the 3 months stub period that he used in his calculation. So the anomaly with his 2006 results can be traced definitely to this (which shows up a “cool” difference in the graph below.)

usgrid57.gif

So what does my calculation using CRN1-2 data look like: see below, which includes a 2007 estimate (using the 2006 difference between the CRN1-2 and NOAA to estimate 2007 CRN1-2 from the YTD NOAA results.) In this calculation, the warmest year remains 1934 by a margin of nearly 0.2 (enough to be “statistically significant” according to things that GISS has said, then 1921, then 1998, then 2006 in fourth place. Values since 1998 have been at uninterrupted high levels. Looking at the graph, it’s definitely noticeable – but is it the sort of change that’s inconsistent with the sort of stochastic behavior in the early part of the series – hard to tell here. If you cut the series off in 1997, you certainly wouldn’t have noticed anything odd about the 1990s. You can sort of see in retrospect why MBH98-99 were so excited about how warm 1998 and in such a rush to include it in their article. The main verification issues right now is why the ROW history is so different from the US – is it climatic or nonclimatic? And can one tell?

usgrid58.gif

By Stephen McIntyre | Posted in Scripts, Surface Record | Tagged , , , , | Comments (108)

2006 and CRN1-2

Oct 3, 2007 – 11:13 PM

John V has posted some graphics recently arguing that CRN1-2 yielded pretty much the same results as major temperature indices and, in some sense, vindicated these results. As Gavin Schmidt has pointed out, the U.S. is only 2% of the world’s surface; and as I’ve observed on many occasions, the statistical methodologies and data quality in the U.S. are different from the ROW, as is the history.

I’ve been focussing more on ROW issues recently rather than U.S. issues, but thought that it would be interesting to drop in briefly on the U.S. issues. Since there seems to be an interest in gridding, I’ll probably do a version, using the methods that I used to contour the trend maps a while ago. I also thought that it would be interesting to verify John V’s CRN1-2 estimate to NOAA and NASA versions myself.

I mention 2006 in the header because there’s something interesting about 2006 which hasn’t been mentioned in John V’s posts: 2006 values for both NOAA and NASA are as much as 0.8 deg C higher than CRN1-2 values.
Continue reading →

By Stephen McIntyre | Posted in Hansen, Surface Record | Tagged , , , | Comments (81)

A First Look at the CRU Station List

Oct 3, 2007 – 5:50 PM

On Sept 28, 2006, Willis Eschenbach sent an FOI for CRU station data. A year later, after many letters, we still do not have the CRU data as used, but do have a list of stations used, a list which is slightly shorter (4138) than the 4349 stations said to have been used in Brohan et al 2006, their most recent publication. It looks as though they sanitized the list somewhat – in Brohan et al 2006, they said that they removed 55 duplicates. I guess that they’ve identified 156 more duplicates (4 times as many as reported in Brohan et al.) But perhaps there’s another reason. In my opinion, they should have delivered a list of 4349 stations – I’ve asked for this from Phil Jones.

Secondly, the list has not been delivered in working order: the data is supposed to be mostly (“98%”) derived form GHCN, but the identification numbers do not tie in precisely to GHCN numbers. The CRU ID numbers are 6 digits versus 11 digits at GHCN. Many of the CRU numbers tie in to GHCN numbers as follows: the GHCN number is in the form CCCWWWWWDDD where CCC is the country code, WWWWW is the WMO number and DDD is the station – for example, nearby sites (but different) sites can have the same WMO number. GHCN DDD identifiers seldom get out of single digits. CRU identifications of WWWWWD tie to GHCN numbers for 1782 sites. As seen below, many sites can be identified with GHCN sites, but not without a further concordance. CRU says that they have a look-up table, but failed to disclose it. I’ve requested it.

Thirdly, and somewhat unbelievably, the CRU identifications are not unique. In one case, there are 6 stations with an identical ID number. Perhaps there is some still undisclosed list and they’ve delivered it in non-working form for some reason of their own. IF there is no proper list, then I have no idea how they can define a look-up table that functions for non-unique ID numbers. For my own attempt at a concordance, I’ve added a duplicate number for each group of sites with overlapping ID numbers so that the new ID number is unique. (I wasted a considerable effort before I figured out that they had non-unique ID numbers – imagine.)

After doing this, in order to make a concordance, for each CRU site unmatched in a first pass, I then selected GHCN stations that were within 1 degree latitude and within 1 degree longitude and had the first 6 letters identical. If there was only one, I declared a match and assigned the GHCN number. This didn’t match as many sites as could be matched, but reduced the unmatched sites to about 354 sites.

I then made an ASCII tab-separated table in which I wrote down the CRU station plus lat, long and altitude and GHCN stations within a degree plus the same information and manually inspected the stations. I probably could have figured another semi-automatic method of reducing it further, but I also wanted to inspect the matches. In many cases, there was a fairly obvious match with the previous method failing due to multiple candidates or spelling variations. In this way, I added 175 matches, getting up to 3959 matches (a little under 96%), leaving 179 unmatched.

Here are ASCII files listing all CRU stations together with proposed GHCN identifications (ID, name, lat, long shown from GHCN) and the unmatched list. All Unmatched These are ASCII tab-separated and can be opened in Excel or read in R.


Unmatched Stations

The distribution of unmatched stations is really very strange – and I add here, that, for each country specified below, I’ve double checked manually against the GHCN inventory to confirm that there was at least one unmatched CRU station from that country. In total, I identified 29 countries where there were CRU stations that were not present at GHCN, including surprisingly stations from Canada, Australia and even the U.S. Here is a list of countries with at least one station that is unmatched in GHCN:

Argentina: CRU had quite a few stations not at GHCN.
Australia: nearly all match, but two CRU stations didn’t match GHCN – Maryborough and Brisbane Airport. Why these?
Austria – quite a few stations not at GHCN
Bolivia – a couple didn’t match
Brazil – a couple didn’t match
Canada – quite a few stations not at GHCN. I noticed a duplicate GHCN for Parry Sound, which is near Toronto and which occurs in two alter egos in GHCN.
Chile – a couple didn’t match. A couple were called “UNKNOWN” in the CRU list. Perhaps they are connected to the UCAR “Bogus Stations”.
China – quite a few stations not at GHCN
Denmark – a couple didn’t match
Dominica – a couple didn’t match
Germany – one didn’t match
Finland – one (Kuopio) didn’t match
Greenland – possibly a couple didn’t match
Guinea – one didn’t match
Iran – one didn’t match
Ireland – one may not match (Phoenix Park)
Israel – a few don’t match
Italy – a couple don’t match
Kyrgyz republic – one doesn’t match
Netherlands – a couple don’t match
Norway – a few don’t match
Oceania — a few don’t match
Peru – one doesn’t match
Sweden – a few don’t match
Syria – several don’t match
Taiwan – quite a few don’t match
UK – a couple don’t match (Kirkwall, Wick)
USA – about 25 don’t match e.g. Moroni, Lahontan
Russia – a couple may not match

IT is quite weird to see these oddball stations crop at CRU. I’m sure we’ll quickly track down where Moroni and Lahontan and their ilk come from, but it doesn’t seem to be GHCN.

Prior Excuses

With these results in mind, let’s review the history of CRU excuses as to why they should not be required to disclose information under the FOI Act – and it’s taken slightly over a year and many letters and appeals to even get this station list. Their original refusal CRU stated that the data was already located at GHCN as follows:

Datasets named ds564.0 and ds570.0 can be found at The Climate & Global Dynamics Division (CGD) page of the Earth and Sun Systems Laboratory (ESSL) at the National Center for Atmospheric Research (NCAR) site at: http://www.cgd.ucar.edu/cas/tn404/ Between them, these two datasets have the data which the UEA Climate Research Unit (CRU) uses to derive the HadCRUT3 analysis. The latter, NCAR site holds the raw station data (including temperature, but other variables as well). The GHCN would give their set of station data (with adjustments for all the numerous problems). They both have a lot more data than the CRU have (in simple station number counts), but the extra are almost entirely within the USA. We have sent all our data to GHCN, so they do, in fact, possess all our data.

In accordance with S. 17 of the Freedom of Information Act 2000 this letter acts as a Refusal Notice, and the reasons for exemption are as stated below

In response to a further request trying to pin them down, they stated that “more than 98%” of CRU data and the remaining 2% was collected under confidentiality agreements.

Our estimate is that more than 98% of the CRU data are on these sites. The remaining 2% of data that is not in the websites consists of data CRU has collected from National Met Services (NMSs) in many countries of the world. In gaining access to these NMS data, we have signed agreements with many NMSs not to pass on the raw station data, but the NMSs concerned are happy for us to use the data in our gridding, and these station data are included in our gridded products, which are available from the CRU web site. These NMS-supplied data may only form a very small percentage of the database, but we have to respect their wishes and therefore this information would be exempt from disclosure under FOIA pursuant to s.41. The World Meteorological Organization has a list of all NMSs.

Obviously, none of this justified not providing a list of stations, but that has taken another 6 months. In connection with the supposed confidentiality agreements, as reported previously, Doug Keenan asked for the countries with which there were confidentiality agreements that restricted access and was told:

Dear Doug,
I have done some searching in files – all from the period 1990-1998. This is the time when we were in contact with a number of NMSs. We have also got datasets from fellow scientists and other institutes around the world. All supplied data (eventually and sometimes at cost), but we were asked not to pass on the raw data to third parties, but we could use the data to develop products (our datasets) and use the data in scientific papers. It is likely that some of the NMSs and Institutes have changed their policies now – and that the people we were corresponding with (all by regular mail or fax) are no longer there or are in different sections. The lists below don’t refer to all the stations within these countries, nor to all periods, but to some of the data for some of the time.
The NMSs
Germany, Bahrain, Oman, Algeria, Japan, Slovakia and Syria

Scientists/Institutes (data for these countries)
Mali, India, Pakistan, Poland, Indonesia, Democratic Republic of the Congo (was Zaire), Sudan and some Caribbean Islands.

These are the only ones I can find evidence for. I’m sure there were a few others during the 1980s, but we have moved buildings twice since 1980.

Not sure how you will use this data.
Phil Jones

Above I summarized the countries for which there are stations that are not matched at GHCN. Remarkably, these include virtually none of the countries where Jones said that they had received data subject to confidentiality agreements – so that the confidentiality agreement excuse cannot apply for any of these countries. And for each of the countries for which Jones said that there was a confidentiality agreement (Bahrain, Oman, Algeria, Japan, Slovakia, Mali, India, Pakistan, Poland, Indonesia, Zaire and Sudan), I was able to cross-identify all CRU stations with GHCN identifications so that the confidentiality excuse didn’t affect anything.

At this point, the only unmatched stations which would appear to be covered by a reported “confidentiality agreement” are about 6 stations in Syria (about half at GHCN) and one German station (Wahnsdorff). Otherwise there is no valid excuse for not disclosing this station data. Of course it is possible that Jones has confidentiality agreements with Canada and the Australia, but was embarrassed to report them and thus omitted them in the above list. We’ll see.

It is disappointing that the pretexts for not providing the data previously have turned out to be untrue. However, it should be possible to now develop a reasonable concordance for the CRU stations to GHCN where applicable and to identify provenances for the oddball stations to make a concordance up to a very small number of stations – at which analysis can begin.

By Stephen McIntyre | Posted in CRU | Tagged | Comments (35)