0

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?

  • Decompose the problem into manageable sub-tasks and implement those. You want a single layer where each pixel value is the layer number with the max value of that pixel? The raster package has many ways to do that. At most brute you can just loop over seq_len(ncell(r[[1]])) and use extract with base::which.max to get the layer for each pixel. – mdsumner Jun 17 '14 at 13:18
2

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

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?

  • Are you running a 64bit instance of R? By default the installer installs both 34 and 64 bit. Be sure that you are launching the R x64 version. The calc function is memory safe in that it looks at your available RAM and adjust the size of the processing block accordingly. However, I do not believe that it looks at the R version being run and just assumes that R can allocate appropriate available RAM. Since there is a memory ceiling in the 32bit version of R this would cause the error you are seeing. – Jeffrey Evans Jun 18 '14 at 2:35
  • My computer has 32bit and R also. – Vandka Jun 18 '14 at 6:27
  • Please edit/modify your original question to expand and do not post an answer! Are you saying that your OS is 32bit? If you are working with spatial data in R a 64bit OS is a must. If you are in a 32bit OS then you may want to play with the rater packages default setting for blockSize. I edited my answer to demonstrate a custom read/write function using blockSize(). – Jeffrey Evans Jun 18 '14 at 17:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.