Calculate anomalies on time-series rasters

Dear all,

I have a netcdf file with monthly temperatures values covering the period
of January 1961 to December 2010:

library(raster)
# Create date sequence
idx = seq(as.Date("1961/1/1"), as.Date("2010/12/31"), by = "month")

# Create raster stack and assign dates
r = raster(ncol=5, nrow=5)
s = stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)))))
s = setZ(s, idx)

First, let's consider only the period from 1961-1990 and take the global
monthly means (aka climatologies)

# Split 1961-1990 period and take climatology
s61.90 = subset(s, which(getZ(s)>=as.Date('1961-01-01') &
getZ(s)<=as.Date('1990-12-31')))
s61.90.mon = zApply(s61.90, by=months, mean, name=month.abb[])

s61.90.mon is a raster with 12 layers, one for each monthly mean
temperature.

Now what I need to do is calculate the monthly 'anomaly' for the remainder
period (Jan 1991 to Dev 2010), i.e. the difference of a single month of the
period from the corresponding month in the climatological mean.

# Here I split the remainder of the time series (1991 to 2010)
s91.2010 = subset(s, which(getZ(s)>=as.Date('1991-01-01') &
getZ(s)<=as.Date('2010-12-31')))

In other words, I need to subtract 'jan' in the raster s.61.90.mon from
every 'jan' in the raster s91.2010, and so on for every other month.


I was wondering what is the best way to do that using raster (and possibly
zoo) functions?

Thank you very much in advance, -- Thiago V. dos Santos

PhD candidate
Land and Atmospheric Science
University of Minnesota