Thursday, August 26, 2010

TempLS Version 2 Release

Here is the long promised Version 2. I posted some trailers months ago:
There are lots of new capabilities, and it's taken a while to get them into a framework suitable for general use. But I think it's there now.

The headline feature is the new spatial capability. This is a logical extenstion of the least squares formulation. But the new organisation scheme may also be attractive. Instead of tinkering with the R script, as formerly, now you just create a user file for each job. The file should be fairly brief, giving just changes to a default user file. I'll list examples.

Code, plus a user guide, is in a TempLSv2.zip file on the doc repository. The files are:
  • TempLSv2.r - the main algorithm file
  • TempLSfnsv2.r - auxiliary functions
  • TempLSuserv2.r - the default user file (acts as a template too)
  • Some sample user files with output directories
  • TempLSprepv2.r - the data preprocessor
  • TempLSprepinputv2.r - a user file for the data preprocessor
  • TempLSv2.pdf - a user guide.
I've put in bold the files that are called by name to run TempLS.


How to run TempLS

  • You'll need R installed, preferably with the maps package.
  • You'll need the four core R files listed above, plus the two user R scripts. These can be quite brief.
  • Run the preprocessor TempLSprepv2.r. You can do this either in R interactive mode, with the command source("TempLSprepv2.r"). Or you can use the command line
    R CMD BATCH TempLSprepv2.r

    The input file will list the station data file (say v2.mean), the inventory, and the SST data (eg HADSST2), and the name you want to give the output files, which are a .data file and a .inv file. Once you have preprocessed the data files, you don't need to do it again.
  • Set jobname="myjobname" and run TempLSv2.r as for ht epreprocessor. If you want to use batch mode, just make a 2-line master script which does what you would do in interactive, and BATCH that.
  • Look for your graphic and text output in a directory myjobname under the directory in which you entered R.

Output

For model 1, for each run there are plots of station numbers and temperature plots (with trend) over the intervals you specify. There's also a map of stations. There's a text output file with the temperature data, and more info on trends, for the intervals you specify.

For models 2 and 3, the output is a shaded plot of trend or temperature, and a plot of station numbers over time.

Where to find input data:

GHCN v2 (v2.mean.Z and v2.temperature.inv)
HADSST2 (HadSST2_SST_1850on.txt.gz (3.9Mb))
GSOD (Ron Broberg - scroll down to Data section)

Example user files

Here is the file for the Nepal post. It runs two cases - one for the country Nepal, and one for the part of the Himalaya in India, Pakistan and Nepal. The two extra cases of high altitude are also given, but there is too little data to run them.

{   # run 1
inputstem="GSOD";
num_runs=2;
name="Nepal"; station_mask=expression(tv$country==217);
model = 1; oCG=F;
yrB=1980;
yr0 = yrB-1; yr1 = 2009;
#model 1
plotyrs = c(yrB,yr1);
trendyrs = c(yrB,yr1,1990,yr1);
}
{   # run 2
name="Himalaya";
station_mask=expression({a=tv$country;(a==217 | a==219 | a==207) & tv$alt>1000 & 14*tv$lat+8*tv$lon>1046});
}
{   # run 3
name="HiHimal";
station_mask=expression({a=tv$country;(a==217 | a==219 | a==207) & tv$alt>2000 & 14*tv$lat+8*tv$lon>1046});
}
{   # run 4
name="HiNepal";station_mask=expression(tv$country==217 & tv$alt>2000);
}


Global Trend and Year map.
Here is the code for the global trend map
Note that it won't look as in that post. This version uses default  contour levels. It was quite fiddly getting a match for the GISS version. The second run plots the 2008 year average, as in this post. Again, the levels don't match.

{
name="Trend12";
inputstem="../TempLS/GHCN";
station_mask=expression(tv$country<9999);
num_runs=2; model=2; num_orthogs=12; isglobe=T;
yr0=1900; yr1=2009; yrx=2008;
type_orthogs=1; # Spherical Harmonics
col_lev=expression({ #Varying the shading colors (optional)
shadcol <- colorRampPalette(c("#ddd4ee","white","#f4de99",
"#edbd6e","#e27d5a","#db6a67","#854746", "#d12e5e", "#8c181d"));
}
{
name="Globyr12";
model=3; num_orthogs=12;
yr0=1900; yr1=2009; yrx=2008;
}


Spatial US year map.

This is the US map for 2008, as in this post.

Again the levels don't match. The year 2008 is not a good choice, since most of the US stations had been moved out of GHCN. Before 2005 is better. However, it's chosen because of the NOAA comparison plot.

{
name="USyr2008";
inputstem="../TempLS/GHCN";
station_mask=expression(tv$country==425 & getdist(c(39,-100))<2700);
num_runs=1; model=3; num_orthogs=6;
yr0=1900; yr1=2009; yrx=2008;
type_orthogs=0;# Trig Functions;
}

0 comments:

Post a Comment