I have 11 NDVI raster data, representing single summer. I want to create raster by day (Julian day or name of raster) that observed maximum value.

In other words, I want to know which day (which raster) has maximum value in each cells of NDVI data and I want to get the result by raster file.

Actually, the day of maximum NDVI value should be in each grid instead of maximum value.

  • Would you be able to edit you Question to focus it on whichever is your preference to use, i.e. R or ArcGIS, please? If that does not get you your first choice Answer then you can always research/ask about the other separately. – PolyGeo Jun 16 '14 at 9:27
  • @PolyGeo It would be preferable to provide a generic answer to a generic question like this one rather than having a proliferation of threads covering the same question in different software environments. – whuber Jun 16 '14 at 19:33
  • @PolyGeo If you need exact one I prefer R. ArcGIS would be ok, if this is impossible in R. – Vandka Jun 18 '14 at 8:28

If you want to use R instead of ArcGIS, I would stick to the suggestion by @radouxju using which.max() from the raster package (which returns an object of class 'RasterLayer'). Like that, you get information about the layer holding the maximum value per pixel, and you can use those values as indices to assign the corresponding Julian day to each cell (provided that you previously extract the Julian day from the raster filenames...). Note that I worked with enhanced vegetation index (EVI) instead of NDVI, but the procedure remains the same.

# Required packages

# List files of interest, e.g. MODIS Aqua EVI files for 2013
fls.evi <- list.files("myd13q1/processed/", 
                      pattern = "^CRP_MYD13Q1.*2013.*EVI.tif$", 
                      full.names = TRUE, recursive = TRUE)

# Extract julian day from available filenames
jdn.evi <- as.numeric(substr(basename(fls.evi), 18, 20))

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

# Per pixel, identify layer holding the maximum EVI value
rst.evi.max <- which.max(rst.evi)

# Replace no. of layer with corresponding Julian day
rst.evi.max[] <- jdn.evi[getValues(rst.evi.max)]
  • It was great functions. I have replaced my data instead of yours. But there is another problem message displayed. In the following line. > rst.evi.max <- which.max(rst.evi) Error in which.max(rst.evi) : not yet implemented for large objects My data is actually very huge as like below 3153, 8597, 27106341, 11 (nrow, ncol, ncell, nlayers). Is there any solution for that? – Vandka Jun 17 '14 at 10:19

With ArcGIS, you can use the spatial analyst tool called highestPosition . It will return the number of the raster where you have the maximum value.

in R, you can use which.max()

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.