1

I want to find the maximum/minimum value from a serials images at each pixel location. These images have the same resolution and projection. such as Z(m)= max(Z1,Z2,Z3....) where Z indicates a raster file.

I know ArcMap raster calculator tool contains a Function 'con'. Are there any rapid methods in Python or GDAL?

2
  • do you want to find the maximum value - from all images - at each pixel location? if so, this question is probably related
    – Steven Kay
    Dec 16, 2016 at 22:51
  • yes, that;s what I mean.@Steven Kay Dec 17, 2016 at 14:09

2 Answers 2

4

If you a have a series of 2D rasters with the same shape, you can stack them to a 3D array:

import numpy as np
Z = np.stack((Z1, Z2, Z3, ...))

Then use min, max or similar along the first axis (i.e. axis=0):

Zmin = Z.min(0)
Zmax = Z.max(0)

Which should result in a 2D array the same shape as one of the input arrays.


For very large files, there are a few ideas. If your 2D rasters are larger than your available RAM, you will need to do multiple passes with window subsets (nXOff, nYOff, nXSize, nYSize) on each raster, which you can do with GDAL. Each window stack would find the minimum, and write the result.

However, if you can hold a few 2D rasters in your RAM at once, you could read them one-by-one and accumulate the result:

shape = (10000, 10000)
Zmin = np.inf * np.ones(shape, 'f')  # Float32
Zmax = -np.inf * np.ones(shape, 'f')
for file in file_list:
    ds = gdal.Open(file)
    ar = ds.ReadAsArray()
    ds = None  # close
    assert ar.shape == shape
    sel = ar < Zmin
    Zmin[sel] = ar[sel]
    sel = ar > Zmax
    Zmax[sel] = ar[sel]
3
  • ok,thanks. firstly need to stack all images to a 3D array in this method.maybe the 3D file is too large to operate. shall we read 2D images pixel by pixel meanwhile write the max pixel and write? Dec 19, 2016 at 1:52
  • @chaobinzhang I've added a few ideas on how to handle large datasets.
    – Mike T
    Dec 19, 2016 at 22:36
  • ,Thanks very much. I think np.stack() is faster than second method. But for np.stack(), how to tranverse a folder and put all fiiles in this Function to stack? Dec 27, 2016 at 14:36
2

The numpy.maximum function creates a new array from two arrays, with the maximum value of the two at each index. It can be applied to 3+ arrays with functools.reduce(). The rasters will need to have the same pixel size and shape, and you can read them as numpy arrays using gdal.ReadAsArray

import numpy as np
# Python 3x only
from functools import reduce

max_arr = reduce(np.maximum, [array_1, array_2, array_3, ...]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.