A classic GIS question you might ask is: "What is the maximum or minimum cell value from a set of multiple overlapping rasters?". When your raster data is in the form of two 2-dimensional NumPy arrays, this is very straightforward as you can just use the numpy.maximum or numpy.minimum functions to get a new array with the minimum and maximum values. However, if you have more than two arrays you need another solution.
To start, lets create the arrays shown in the picture at the top of this post:
import numpy as np
# Create some arrays to stack
a = np.array([[0,1,2],[3,4,5],[6,7,8],[4,6,7]])
b = np.array([[5,4,8],[6,3,9],[1,5,9],[0,0,0]])
c = np.array([[5,1,9],[5,9,7],[3,1,5],[3,9,1]])
The next step is to use the numpy.dstack function to stack the arrays depth wise - which in a GIS context is to say on top of each other:
# Stack the arrays
abc = np.dstack((a,b,c))
It is at this point that things get a little more complicated. What we can do is use the numpy.argmax or numpy.argmin functions to pull out not the minimum and maximum values themselves, but rather an index to the location of the minimum and maximum values within the array stack. The axis argument is set to 2 so that the index is for minimum and maximum values going depth wise into the stack (the first two axes would relate to looking at the stack in terms of top to bottom or left to right):
# Get the index values for the minimum and maximum values
maxIndex = np.argmax(abc, axis=2)
minIndex = np.argmin(abc, axis=2)
Having got the index values for the minimum and maximum values in terms of their depth into the array stack, we just need to create two other arrays that can act as the matching row and column index positions:
# Create column and row position arrays
nRow, nCol = np.shape(a)
col, row = np.meshgrid(range(nCol), range(nRow))
And then its a simple process of indexing out the maximum and minimum values within the array stack to create new arrays containing the minimum and maximum values as shown in the picture at the top of this post:
# Index out the maximum and minimum values from the stacked array based on row
# and column position and the maximum value
maxValue = abc[row, col, maxIndex]
minValue = abc[row, col, minIndex]
Putting the whole process together we get:
import numpy as np
# Create some arrays to stack
a = np.array([[0,1,2],[3,4,5],[6,7,8],[4,6,7]])
b = np.array([[5,4,8],[6,3,9],[1,5,9],[0,0,0]])
c = np.array([[5,1,9],[5,9,7],[3,1,5],[3,9,1]])
# Stack the arrays
abc = np.dstack((a,b,c))
# Get the index values for the minimum and maximum values
maxIndex = np.argmax(abc, axis=2)
minIndex = np.argmin(abc, axis=2)
# Create column and row position arrays
nRow, nCol = np.shape(a)
col, row = np.meshgrid(range(nCol), range(nRow))
# Index out the maximum and minimum values from the stacked array based on row
# and column position and the maximum value
maxValue = abc[row, col, maxIndex]
minValue = abc[row, col, minIndex]
