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?