Sunday, April 4, 2010

More GHCN results.

I've put up version v1.1 of my least squares GHCN code. It now makes it possible to easily analyse subsets using any of the information on the inventory file. I've used it to calculate a number of much discussed examples. The code now runs very fast. On my fairly old PC, the biggest case (global) takes less than 100 secs. So although each special case is calculated pretty much from scratch, there's not much penalty there. The code is here, and also on the document repository.

New code features.

The program has now been split into two. The first is a preprocessor which reads v2.mean, does the merging of duplicates, and prints a condensed file called v2mean.txt. Unlike ver 1, it doesn't require prior editing of v2,mean. The idea is that once you have merged the duplicates, you don't need to do it again un til you want to deal with a new version of v2,mean.

The second part is somewhat like V1, but much more clearly laid out. It's even shorter (with the preprocessing gone) - about 170 lines. There's a lot of infrastructure with printed tables, plots with trends etc.

The selection process works thus. The 18 parameters from the .inv file have been given names, like tv.$country, listed in colnames. You declare a list of names that you'll want to use and have appear on graphs etc: 
title = c("Globe","NH","SH","ConUS","Airport","Rural","Urban"); 
choices = 1:7 #c(4,5,6,7);

In the second line, you declare the ones you want to run.
The logic of the choices is expressed in the list:
if (tc == "NH")   i01 = i01[cellnums>1296]; 
if (tc == "SH")   i01 = i01[cellnums<1297]; 
if (tc == "Airport")  i01 = i01[tv$airpt == "A"]; 
if (tc == "Rural")   i01 = i01[tv$urban == "A"]; 
if (tc == "Urban")   i01 = i01[tv$urban != "A"]; 
That's my current list, but you can change it.

Some results

Firstly, a table of trends, in deg C/decade.
No. Stations..728060561224239027124568

No real surprises. The main outlier is the SH, but remember, this is land stations only. The classifications of airport, rural and urban make little difference. The urban/rural differentiation is by night lighting, and I've included small towns with urban.

Remember that this is based on raw GHCN data, with no corrections made by me.

Update - there was a bug in V 1.1, which had a minor effect on the results, but which could cause a crash. It didn't affect global calcs - only subregions. Here are the recomputed trends:
No. Stations728060561224239027124568
I'll fix the posted code in all locations, calling it v1.11, and the results file.

I've posted the results file on the document repository.

Some plots

NH 1900-2009

SH 1900-2009

NH 1978-2009

SH 1978-2009

Rural 1900-2009

Urban 1900-2009

Rural 1978-2009

Urban 1978-2009

Airport 1900

Airport 1978

1 comment:

  1. Nice work Nick... something else you may want to do: add uncertainties to your series and trend fits. Part of doing the trend fit is remembering that your uncertainty generally is increasing as your number of stations drop (especially before 1950), i.e., so really one also needs to weight by year, if you fit to an interval before that.

    Of course, this is already done in CRUTEM3, so if you haven't seen it, here is a plot including 95% uncertainty bounds