rfunctionrasterquantileterra

Reclassifying a large rasterstack in Terra based on another raster


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: enter image description here and then take the 90th quantile: enter image description here

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))

Solution

  • 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.