I'm trying to do the following:
I have a stack of monthly temperature rasters, and from this I calculate the 90th quantile of the stack.
Then, I want to have a count for how often each raster cell is equal to or above this quantile. This way I'll end up with a raster with summed values for the frequency with which the 90th quantile is exceeded. I'm aware this may seem unnecessary since I'm taking the 90th quantile, but it will make it easier to then repeat the process for a stack of projected (future) rasters, so I can look at the change in frequencies of exceeding the threshold between past and present.
So I start with: and then take the 90th quantile:
I then try to look at each layer of hist_temp_stack and want, for each cell in the raster, the count where the layer exceeds the value for T90:
hist_extreme_temp <- app(hist_temp_stack, fun = function(x) ifelse(x >= T90[], 1, 0))
hist_extreme_temp_count <- sum(hist_extreme_temp)
resulting in
There were 50 or more warnings (use warnings() to see the first 50)
> warnings()
Warning messages:
1: In x >= T90[] : longer object length is not a multiple of shorter object length
And after all those warnings...
> hist_extreme_temp_count
class : SpatRaster
dimensions : 180, 360, 1 (nrow, ncol, nlyr)
resolution : 1, 1 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +no_defs
source : spat_362c1f613e88_13868.tif
name : sum
min value : 4
max value : 62977
But this max value does not make any sense. The max value should not be above 492 (the number of layers in hist_temp_stack), as I'm trying to sum, for each cell, the number of layers in hist_temp_stack that exceed the 90th quantile (T90).
To avoid the warnings related to hist_temp_stack being longer than T90, I've also done the following, but this results in a memory error saying I'd need 8TB of memory...
T90 <- rep(T90, nlyr(hist_temp_stack))
hist_extreme_temp <- app(hist_temp_stack, fun = function(x) ifelse(x >= T90[], 1, 0))
Example data
library(terra)
r <- rast(nrow=5, ncol=5, nlyr=100, xmin=0, xmax=1, ymin=0, ymax=1, crs="local")
r <- init(r, runif)
Get the 90th quantile
q <- quantile(r, .9)
q
#class : SpatRaster
#dimensions : 5, 5, 1 (nrow, ncol, nlyr)
#resolution : 0.2, 0.2 (x, y)
#extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter)
#source(s) : memory
#name : q0.9
#min value : 0.8204033
#max value : 0.9477986
Sum the number of cells with a value >= the quantile
sum(r >= q)
#class : SpatRaster
#dimensions : 5, 5, 1 (nrow, ncol, nlyr)
#resolution : 0.2, 0.2 (x, y)
#extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
#coord. ref. : Cartesian (Meter)
#source(s) : memory
#name : sum
#min value : 10
#max value : 10
The way you were going about it is not valid. The function applied to app should not access the data from another SpatRaster. For something like that you could instead use lapp
with a SpatRasterDataset (one dataset is the temperature data, the other the quantile); but in this case that is unnecessary complex given the above.