I have monthly rasters for multiple years from which I would like to compute seasonal sum based on specific months (October to February).
library (terra)
#create rasters
r1 <- rast(nrows=50, ncols=50)
rr <- lapply(1:36, function(i) setValues(r1, runif(ncell(r1), min = 0, max = 1000)))
#create raster stack
stk <-rast(rr)
#create date vector
dte <- seq(as.Date("2015-1-1"), as.Date("2017-12-31"), by="month")
#apply date vector to raster names
names(stk) <- dte
Using the months assigned in the raster names, I would like to create a new raster stack by computing season sum for five months every year starting from October to February.
There are examples of extracting specific month from the stack Subsetting a Rasterstack by month or computing average for trisemester with seasons going across the years R How to calculate the seasonal average value for each year using stackApply?. In this case, the indices are of equal length.
**indices <-**
seasonal_sum <- stackApply(stk, indices, sum)
Wondering how to get indices that are unequal across years before applying stackapply (five months - October to February and seven months - March to September). Appreciate your help if there better ways to do this
You can use terra::tapp (the equivalent to raster::stackApply)
library (terra)
set.seed(1)
stk <- rast(nrows=50, ncols=50, nlyr=36)
values(stk) <- runif(size(stk), min = 0, max = 1000)
Assign the time (it would be possible to use the names as well, but this is simpler)
time(stk) <- seq(as.Date("2015-1-1"), as.Date("2017-12-31"), by="month")
I would first remove the months that are not of interest.
dte <- time(stk)
m <- as.numeric(format(dte, "%m"))
i <- (m < 3 | m > 9)
x <- stk[[i]]
Now sum the values by season (Oct - Feb). To do so, I shift January and February to the preceding year.
d <- time(x)
y <- as.numeric(format(d, "%Y"))
m <- as.numeric(format(d, "%m"))
y[m<3] <- y[m<3]-1
z <- tapp(x, y, sum)
z
#class : SpatRaster
#dimensions : 50, 50, 4 (nrow, ncol, nlyr)
#resolution : 7.2, 3.6 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84
#source : memory
#names : X2014, X2015, X2016, X2017
#min values : 5.9095, 600.6971, 582.8348, 128.1672
#max values : 1983.775, 4843.461, 4422.809, 2896.181
Below is an alternative route that would be less efficient if you only care about the Oct-Feb season.
d <- time(stk)
m <- as.numeric(format(dte, "%m"))
y <- as.numeric(format(d, "%Y"))
y[m<3] <- y[m<3] - 1
m <- ifelse(m < 3 | m > 9, "w", "s")
ym <- paste0(y, m)
zz <- tapp(stk, ym, sum)
subset(zz, grep("w", names(zz)))
#class : SpatRaster
#dimensions : 50, 50, 4 (nrow, ncol, nlyr)
#resolution : 7.2, 3.6 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84
#source : memory
#names : X2014w, X2015w, X2016w, X2017w
#min values : 5.9095, 600.6971, 582.8348, 128.1672
#max values : 1983.775, 4843.461, 4422.809, 2896.181
For others looking at this question. For the simpler case where you just want to summarize values by calendar year you can do (with terra >= 1.6.0)
r <- tapp(stk, "years", sum)