[GIS] Calculating median of all pixels of all layers of raster stack in ArcGIS for Desktop

Question

I have a stack of 15 layers and I need to calculate the median of all the pixels of all the layers.

How can i do so in ArcGIS?

I tried zonal statistics but it only calculates the median among pixels in one layer.

I also tried Cell Statistics but then it calculates the median of every single pixel trough time.

I need to consider all the values of all the pixels of all the layers together and calculate the median of that distribution.

Answer 1

The following code should work. With a few lines more you could also calculate the percentiles you need.

import arcpy
import numpy as np

arrays = []    
arcpy.env.workspace = r"C:\..."  # The directory containing your rasters

for raster in arcpy.ListRasters():
    arr = arcpy.RasterToNumPyArray(in_raster=raster)
    arrays.append(arr)

stack = np.stack(arrays)   # 3d numpy array
median = np.median(stack)  # returns one single value