[R-sig-Geo] Calculate anomalies on time-series rasters

[Loïc Dutrieux]

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
>>