Sunday, April 18, 2010

V1.4 with maps, conjugate gradients

I'm posting version 1.4 of TempLS. It has some new features. One is computational - it uses a conjugate gradient solver instead of the direct solver. This is much faster - at least five times, and more than doubles overall speed. The solution step goes from being the big time consumer to taking just a few seconds for even the biggest problems.

The next is a mapping facility. To use it, you'll need to install the "maps" package. Installing a package in R is easy - just type (with R running)
and answer queries about mirror sites etc. With linux, this is better done as root.

If you don't want to do this, there's a variable at the top of the GlobTempLSv1.4.r called UseMaps - set it to F and everything but maps will be done.

It provides for each selection that you run a world map with the stations that qualified marked as dots. It's a useful check.

V1.4 also provides plot of number of stations in your sample, by year. This helps to keep from overinterpreting when numbers get small.

V1.4 also corrects an error with the ocean modified data, which had the years mislabelled by a year - it said 1881-2009 when it was actually showing 1880-2008. That made a minor difference to results. Plots look not much different now, but slightly better.

The selection function chooser() is simplified using a R switch statement. Now you just add lines to the list like:
"SST" = tv$urban=" ", # the last comma is needed
This works here because for sea "stations" many fields like urban are blank. As usual, the code is on the repository here, and at the document store.

A summary of this series of posts is here.

An SST example

Here's a map of the artificial SST stations that are used in the global analyses, and a plot of station numbers over the years. Note the WW1 and WW2 dips

A curiosity is that the Great Lakes are counted as SST - I didn't know that. So lets focus on that. Detroit's lat/lon is about (42,-83), and I'm guessing that a circle of 1200 km radius would pretty much cover them. So I'll try a new criterion:
"GtLakes" = tv$urban=" " & getdist(42, -83),

Here's what the metadata shows:
So there were almost no "stations" until 1960, and after that about 6. The map looks a bit odd because it doesn't plot with the1200 km radius requested, but only as much as needed to cover the stations, which are at the centres of their 5x5 grid cells. This actually makes them appear to be on land, but locating them at the centre of their cells is only a convention. So despite appearances, the full area of the  Lakes are included. Adjusting the plot and trend request accordingly:
So the Great Lakes appear to be warming in recent years, at about the global rate, although the patchy data does not justify an emphatic statement.


  1. very nice work Nick. I finally had a break through on the duplicate stations with the help of guys on the R help list. There are some nice matrix tools coming in a new package MatrixStats, but for now I didnt use it. i think I may be able to use some of Steve's code to replicate what GIss do in combining duplicates. For now I've implemented two approaches:

    for a set of duplicates ( say d1, d2, d3 ) I have two options:

    1. Simple average: (d1,d2,d3)
    2. Midpoint of the range (max(d1,d2,d3)+min(d1,d2,d3))/2.

    Glad to see you added the Map stuff. If your interested I can shoot you the kml code I got from peter oneill, It needs to be cleaned up
    but if you feed it a file of locations and station names it creates
    a google map with pins. Let me know if you want it.

  2. Thanks, Steven - yes, I'd be glad to look at the KML file. I've never really done anything with xml. I thought of adding names to the map, but it's hard to avoid a lot of clutter. Also cities often have several sites.

  3. Steven, I've added a contact link under resources on the right.

  4. Nick,

    Nice work! I ran my own model with SST data, and our results are nearly identical. However, there does seem to be some oddities with the SST data we are using, at least compared to Hadley SST. Specifically, it seems to be one year off when you plot them side by side. Working to make sure it wasn't an error on my part...

  5. Zeke,
    Yes, the data were out by a year. I had left a comment back on the original thread. I've left a new zip file on the doc repository

  6. Nick,

    Can you update the output-2010-04-13-22-00-landsea-sst-nhlandsea-shlandsea.txt file with a run using the corrected SST data? Finishing up a comparison chart for global temp reconstructions, and I don't want to use a messed up version of your model.

  7. Zeke,
    Yes, I'll do that. It will be an hour or so before I can get on to it.

  8. Zeke,
    I've uploaded the new file. The old one is called oldLandSea...
    I've been puzzling over an oddity with the Had files I've been comparing with, which came from WFT. None seem to have a zero average over 1961-1990. The GISS file does have a zero average over 1951-80.

  9. Thanks Nick!

    I just rebaseline everything before plotting to be safe.