Raster algebra on multiple bands in R Growing degree days calculation

3
votes

So I am building a 30-year average growing degree days raster from the daily min and max temperature for North America. The input data is composed of 60 rasters, 30 tmax, and 30 tmin . Each of those rasters consists of 365 bands (one for each day). The calculation for GDD is fairly simple it is ((tmax-tmin))/2)-10. The problem is I can't seem to find a way for all the bands and raster to be calculated together.

My initial plan was to stack all tmin rasters and all tmax rasters and subtract them together. But in that case, only the first band (day 1) gets computed. Is there any ways to do the operation without having to extract all bands one by one.

library(raster)

tmax <- stack("C:/gddtmax")
tmin <- stack("C:/gddtmin")
average <- ((tmax-tmin)/2)-10
1
rstackrastermodelingr-raster
hi @PhilippeTremblay, I am trying to calculate something similar. Please, is here your base temperature 10? I.e. then you calculate the sum of the positive temperature differences? if is is sum over 100 days of positive difference, can you call it temperature sum (100 dd)? – maycca

1 Answers

2
votes

Your example works fine for me.

library(raster)
f <- system.file("external/rlogo.grd", package="raster")
tmax <- brick(f) / 10
tmin <- tmax - 10

tmax and tmin each have three layers

avg <- ((tmax-tmin)/2)-10

And avg also has three layers

avg
#class       : RasterBrick 
#dimensions  : 77, 101, 7777, 3  (nrow, ncol, ncell, nlayers)
#resolution  : 1, 1  (x, y)
#extent      : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#data source : in memory
#names       : red, green, blue 
#min values  :  -5,    -5,   -5 
#max values  :  -5,    -5,   -5

EDIT It looks like the problem was due to a memory management bug (fixed in version 2.7-3). Here are two work-arounds.

Use overlay (which also allows you to specify an output filename)

avg <- overlay(tmin, tmax, fun=function(x, y){ ((y-x)/2)-10 })

Or use rasterOptions

rasterOptions(todisk=TRUE)
avg <- ((tmax-tmin)/2)-10
rasterOptions(todisk=FALSE)