Error in which.max() : not yet implemented for large objects in R

Question

I am going to create raster by Julian day that observed maximum value of NDVI within a single year. Actually, I want to know which julian day (or raster) has max value in each cells and get result by a single layer of raster. I am working with huge raster stack (3153, 8597, 27106341, 11 (nrow, ncol, ncell, nlayers)). I have found following script to do that, but it reported error message as like "Error in which.max() : not yet implemented for large objects".

Is there any solution for this kind of error?

# Required packages
library(raster)
library(rgdal)

# List files of interest, e.g. MODIS Terra NDVI files for 2000
fls.ndvi <- list.files("D:\\NDVI_MODIS_500\\Mod_ndvi_500_real_16", 
                      pattern = "mod13a1.a2000.*.*.img$", 
                      full.names = TRUE, recursive = TRUE)
fls.ndvi

# Extract julian day from available filenames
jdn.ndvi <- as.numeric(substr(basename(fls.ndvi), 14, 16))
jdn.ndvi

# Import files into R
rst.ndvi <- stack(fls.ndvi)
rst.ndvi

# Per pixel, identify layer holding the maximum EVI value
rst.ndvi.max <- which.max(rst.ndvi)   
### Error in which.max(rst.ndvi) : not yet implemented for large objects 

# Replace no. of layer with corresponding Julian day
rst.ndvi.max[] <- jdn.ndvi[getValues(rst.ndvi.max)]

What should I do to avoid above errors?

Answer 1

This is a very good example of why we request a repeatable code example. Without seeing how you are using which.max() we can only speculate on what the issue is.

I am assuming that you are trying to call the raster wrapper function by passing a brick object directly to which.max(). Rather than using the raster wrapper function, I use the base which.max() function in conjunction with calc() or overlay().

Here is an example of using which.max and calc() to return the layer index of the maximum brick/stack raster brick/stack cell value. To keep everything memory safe you will likely need to write the resulting raster to disk using the "filename" argument in calc().

require(raster)

# Create example data
r <- raster(ncols=100, nrows=100)
  r[] <- runif(ncell(r)) 
    rb <- brick(r)      
      for(i in 2:11) {
        r[] <- runif(ncell(r))
          rb <- addLayer(rb, r) 
      }

# Use calc to return max raster layer index   
max.index <- calc(rb, function(x) {which.max(x)})
  plot(max.index)   

If you need to explicitly control the size of the read/write chunks you can modify the default blockSize and write custom code for read/write. Here is an example expanding on the code above.

# CONTROL READ/WRITE BLOCKSIZE 
( tr <-  blockSize(rb) )
  s <- writeStart(rb[[1]], filename="test.img", overwrite=TRUE)
    for (i in 1:tr$n) {
      v <- getValuesBlock(rb, row=tr$row[i], nrows=tr$nrows)     
    writeValues(s, apply(v, MARGIN=1, FUN=which.max), tr$row[i])
    }
  s <- writeStop(s) 

Answer 2

I have replaced following function that has error

# Per pixel, identify layer holding the maximum EVI value
rst.ndvi.max <- which.max(rst.ndvi)   
### Error in which.max(rst.ndvi) : not yet implemented for large objects 

by

### Cause of this error I have used following functions.
max.index <- calc(rst.ndvi, function(x) {which.max(x)})
### 1. Error: cannot allocate vector of size 569.3 Mb  

Still I have errors as an above. What should I do to avoid these errors?