Home > Uncategorized > SST with Raster. Complete

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.

About these ads
Categories: Uncategorized
  1. p
    August 30, 2010 at 3:35 PM | #1

    Is this possible to change color scale?

    • Steven Mosher
      August 30, 2010 at 8:15 PM | #2

      I’m sure it is, but havent looked into it yet. Graphics is the last thing I tend to do.

  2. p
    August 31, 2010 at 11:04 AM | #3


    I’ve changed drawGlobalbrick function:

    jet.colors <-colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

    and change:

  3. steven Mosher
    August 31, 2010 at 12:53 PM | #4

    ok, lemme test and revert.. works.

    I would like to figure a way to get a consistent scale month to month

  4. steven Mosher
    August 31, 2010 at 1:02 PM | #5

    zlim I’ll play a bit

  5. p
    August 31, 2010 at 1:18 PM | #6

    And i have last question – temp. scale is different for each year. Is possible to change it and set one scale (from -6 to +6 for example)?

  6. p
    August 31, 2010 at 1:37 PM | #7

    oops i forgot to refresh page :)

  7. p
    August 31, 2010 at 1:49 PM | #8

    Sorry for this spam, but yes – zlim works.

    simply add zlim=c(-6,6) param to plot function.


  8. steven Mosher
    August 31, 2010 at 2:09 PM | #9

    P I added -10 to 10, lemme post up some stuff and you tell me. New post.

    Also need to see how to supress labels on the map.. you’ll see

  1. February 25, 2011 at 8:40 PM | #1

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s


Get every new post delivered to your Inbox.

Join 31 other followers

%d bloggers like this: