0

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

karthikeyan matheswaran
  • 5
  • 2

1 Answers1

0

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)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63