SST with Raster. Complete

Update: new zip, correcting bug found by Steve  McIntyre:

if(!file.exists(HadSST2ncdf)) downloadHADSST2()

if(!file.exists(HadSST2ncdf)) downloadHadSST2()

issue pending with another line as well. Checking raster versions. I’ve also, added some code into “downloadHadSST2″ that corrects for the “NA” problem with HadSST. (currently commented out). There is an issue with “ncdf” handling CF standards, which has been addressed in raster (need to test) until the ncdf folks address it.  So, you can try uncommenting that code to do the fix manually,like I used to do.

UPDATE: try out the box.net on the right hand sidebar to get the two scripts.

UPDATE: NA issue with the maps.. the little red dots on Decade maps are actually NAs. Tracking that down.

UPDATE: Windows and Linux users should switch to 1.4.10. Available on CRAN, MAC is coming. (or get it from Forge)

Two lines will change as applyStack() will now take na.rm as a parameter. less filling and tastes great!

Annual <- stackApply(SST,indices=yeardex,fun=function(x)mean(x ,na.rm=TRUE))

Annual <- stackApply(SST,indices=yeardex,fun=mean, na.rm=TRUE)

Decade <- stackApply(Annual,indices=decadeIndex,fun=function(x)mean(x ,na.rm=TRUE))

Decade <- stackApply(Annual,indices=decadeIndex,fun=mean ,na.rm=TRUE)

Update: Admiring my own stuff, I found this:

decadeIndex <- rep(1:16, by=10)

decadeIndex <- rep(1:16, each=10). Needs testing the first is not correct

Headed out for bit, test on MAC when I return. Latest version is in the box to the right.

Analysis of SST processing with R’s ‘raster’ package is now complete.  The process has gone really smoothly primarily because of some help by the key package developer. My goal was to complete SST processing using only raster objects and only raster functions. When you download the zip file and unzip it, you will see the following:

Those are the only two files you should have in your directory. Make this folder your working directory and make sure that you have all the packages installed. This will not work properly unless you have the lastest version of raster from  R Forge. 1.4.9. The other pacakges you need should be clear from the library load statements in the R Script: SSTVerify.R. With all the packages installed, just execute the script, SSTVerify.R

That script will download and unzip all the files you need. It will download the SST datafile ( which takes a while ), a land/ocean mask, and the latest results from Hadley that we will validate against. When the program completes you will se your directory populated like so:

While the program runs you can watch the progress as files fill up.  I’ll go over the files in order as they are written ( the ascci files are self explanatory. Arrg,doubled the date columns.. easy fix). The first calculation we do is a simple unweighted average. This is compared to the figures from Hadley.Next we weight the temperatures by the area of grid cells. small difference is primarily from a difference in how they calculate a global average. There average is the mean of the NH and SH. Moshtemp calculates a global average. Next I take the unweighted temps and generate Annual area maps for every year in the series ( and yes my colors suck ). Then I collate by decade.Finally I test applying a weight for  the actual area  of water in a cell. This would ONLY be applied when doing a complete land/ocean analysis. There we would have temperature from land stations and from Ocean data sets in the cells that are coastal. So, don’t make too much of the difference between these series. What is shows is that coverage is increasing in cells that are coastal.

To show you how powerful raster is, We can calculate the monthly SST like this in 3 lines

SST <- brick(HadSST2ncdf,varname=’sst’)

Weights <- area(SST,na.rm=T,weight=T)

Monthly <- cellStats(SST*Weights,sum)

and Monthly will contain the mean monthly SST for 1850 to 2010/7

OR in 2

SST <- brick(HadSST2ncdf,varname=’sst’)

Monthly <- cellStats(SST*area(SST,na.rm=T,weight=T),sum)

What this means is that when the land is rewritten we will be able to combine a land series and SST series in a few simple lines.