Calculate anomalies on time-series rasters
Thanks for the input, Michael and Loic. It looks like remote's "anomalize" might work, but the overlay approach was more convenient to use based on the shape of my data. Greetings, -- Thiago V. dos Santos PhD student Land and Atmospheric Science University of Minnesota
On Sunday, February 12, 2017 5:24 PM, Lo?c Dutrieux <loic.dutrieux at conabio.gob.mx> wrote:
On 12/02/2017 03:14, Michael Sumner wrote:
I believe the "remote" package has functions for doing exactly this. HTH On Sun, Feb 12, 2017, 17:42 Thiago V. dos Santos via R-sig-Geo < r-sig-geo at r-project.org> wrote:
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?
If you're confident about the order of the layers in your rasterStacks
you could use overlay. See the note about recycling in the function help.
# lenght(x) needs to be a multiple of length(y), so that it can be recycled
# e.g. c(10, 10, 10, 10, 10, 10) - c(5, 7)
fun <- function(x, y) {
x - y
}
annomalies <- overlay(x = s91.2010, y = s61.90.mon, fun = fun)
However, it looks like your climatologies are in alphabetic and not
chronological order, so you'll have to double check that.
Cheers,
Lo?c
Thank you very much in advance, -- Thiago V. dos Santos PhD candidate Land and Atmospheric Science University of Minnesota
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo