# Determining Highest Position in QGIS for raster stacks with 10+ images?

I want to follow up on this question from last year

Emulating ArcGIS Highest Position tool in QGIS?

Here the question is if QGIS has a tool like the "highest position" tool in ArcGIS, where, in a layer stack, the tool will return the position or index value of the raster with the highest value for each cell. For example, if the third raster in a 10 raster stack has the highest value for a particular cell, than that cell will have the value 3.

The solution in the above question, rather elegantly, uses the Raster calculator. However, if you're trying to find the highest position for all cells in a layer stack with 10 or more images, this particular solution becomes difficult. Is there any other solution to this in QGIS that can be used for raster stacks with 10+ images?

I found a solution using R, and wanted to report that, but would like the question to stay open, because I think it is important to see if there is a solution to this using QGIS.

Anyway, through some web-searching, I found this email thread on an R site, asking more or less the same question as I did.

https://stat.ethz.ch/pipermail/r-sig-geo/2013-October/019572.html

It looks like the author of the R raster package (among others) offered insights here. So, anyway, using the Raster package on R (and Rgdal), the following code, which I have appropriated directly from the link above, worked for my purposes:

The following code, using "stackApply", produced a max value composite of the layer stack, with the number of "1's" in the "indices" corresponding to the number of layers in the stack.

max_NDVI <- stackApply(raster_data, indices=c(1,1,1,1,1),max)

The following code produced an index of max values, such that the pixel value = the number of the layer in which the max value occurred for that pixel in the layer stack. You have to create a function using "which.max".

miwhichmax <- function(x,na.rm=TRUE) which.max(x)

index_max <- stackApply(raster_data,indices=c(1,1,1,1,1),miwhichmax)

As I said, this suited my purposes, but I'm not sure how good this approach would be for someone who has a truly large layer stack, for example if you're trying to index daily MODIS data.