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

[Thiago V. dos Santos]

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