Creating a sea surface temperature animation in R

We’re going to use R and the NOAA optimal interpolation sea surface temperature data set to create an HTML/Javascript animation of global sea surface temperatures, with the continents overlaid on top.

The data

We’ll be using the NOAA OI-SST v2 data. You can grab the data in NetCDF format from that previous link. You’ll want to grab both the SST data and the land mask data that describes the placement of the continents on the same coordinate system.

The data is 1 by 1 degree gridded data of the estimated sea surface temperature at a given point in time. The numeric data are:

  • latitude in degrees north (positive) or south (negative)
  • longitude in degrees east from 0ยบ, which means that 80W corresponds to 360-80 = 280E
  • date as the number of days since 1 Jan 1800

They provide weekly and monthly data; I chose weekly for myself (sst.wkmean.1990-present.nc).

Setting up R

Make sure you have the ncdf and animation libraries installed.

The code

The following R code will load the files from the working directory and save the animation to a temporary directory. See the ani.options documentation for information on changing the defaults.

# Load the data-handling libraries
library(ncdf)
library(chron)

# Load SST and landmask data
sst.oisst <- open.ncdf("sst.wkmean.1990-present.nc")
landmask  <- open.ncdf("lsmask.nc")

x <- get.var.ncdf(sst.oisst, "lon") # longitude for the x axis
y <- get.var.ncdf(sst.oisst, "lat") # latitude for the y axis

# z holds the SST values, accessed as z[x_index, y_index, t_index]
# Note that these are the indexes in the x, y, and t vectors that
# correspond to the location of interest. If you want 55.5W 39.5S at
# the week centred around 9 Dec 1990, the appropriate call is
# z[205, 130, 50]
z <- get.var.ncdf(sst.oisst, "sst")

# The dates are in number of days since 1 Jan 1800
t <- get.var.ncdf(sst.oisst, "time") # t holds the time points
t <- chron(t, origin=c(month=1,day=1,year=1800))

# Grab the land mask
mask <- get.var.ncdf(landmask, "mask")


# Create the animation
library(animation)

# Start recording the animation
# Basically, each time you clear the device buffer, it saves that as a frame in the
# animation.
ani.start()

# For every frame (that is, every time point in t)...
for (i in 1:length(t)) {
  # ...draw a contour map of the SST data
  # filled.contour expects x and y to be monotonically increasing, so we must
  # reverse y. In order for the z data to be drawn correctly, we also need to 
  # reverse is on its y axis.
  filled.contour(x, rev(y), z[,ncol(z):1,i],
                 color=function(x)rev(heat.colors(x)),
                 asp=1,
                 levels=pretty(c(-2,36), 20),
                 plot.title=t[i],
                 plot.axes = {
                   # This runs after generating the contour plot, using the contour plot's axes
                   # for its coordinate system. We use this to overlay other objects onto the plot

                   # Add the land mask as an 'image', with a contour plot to smooth the edge
                   # Don't forget to reverse the y and mask vectors appropriately, and
                   # make sure to 'add' them to the current buffer instead of flushing
                   # the buffer first.
                   image(x, rev(y), mask[,ncol(mask):1], add=TRUE, col=c("black", "transparent"))
                   contour(x, rev(y), mask[,ncol(mask):1], add=TRUE, drawlabels=FALSE, nlevels=2)

                   # Add the date to the bottom of the frame
                   text(180, -105, t[i], pos=1)
                 }
                 )
}

# Stop recording the animation
ani.stop()

Results

SST Animation

The resulting animation can be seen here—give it some time to load, it’s about 65MB of PNGs. Once it’s loaded, you can click the ‘>’ to start it playing, and you’ll probably want to click ‘+’ a few times to speed it up.