enter image description hereI am working with MODIS Global Primary Production (GPP) rasters and I am trying to pick the max GPP pixel value in a series of 15 dates for a study area and create a new raster layer. MODIS GPP rasters contain values greater than 30,000 that are either NoData, unclassified or other categories that do not correspond to GPP values. I am actually looking for the max value between 0 and 30,000.

gdal_calc.py will only take 2 rasters as input and while SAGA accepts more than two rasters as input, it does not let specify a range for values to be taken into account.

Any suggestions how this can be done?

1 Answer 1


You can accomplish this using gdal and numpy. First you have to read all your rasters as numpy arrays and then stack them to create a single 3D array. Then you can mask this array to exclude values outside of the 0 - 30,000 range and then you can get the maximum value for each pixel. Finally, you have to fill the remaining masked values (pixels where the values for the 15 dates were outside the range) with a NoData value and create the output GeoTIFF.

import glob
import os

import gdal
import numpy as np

# change working directory

# get rasters' file names
fnames = glob.glob('*.tif')

# read general properties of the first raster (assuming all the rasters share these properties)
ds = gdal.Open(fnames[0], 0)
gt = ds.GetGeoTransform()
sr = ds.GetProjection()
xsize = ds.RasterXSize
ysize = ds.RasterYSize
nd = ds.GetRasterBand(1).GetNoDataValue()
del ds

# read each raster and create a 3D array
arrays = []
for fn in fnames:
   ds = gdal.Open(fn, 0)
   arr = ds.ReadAsArray()  # 2D array (rows by columns)
   del ds
arr = np.stack(arrays)  # 3D array (date by rows by columns)

# mask the array to exclude values outside a defined range
mask = (arr < 0) & (arr > 30000)
arr = np.ma.array(arr, mask=mask)

# get maximum value for each pixel
arr = np.max(arr, axis=0)

# fill the array with the NoData value
arr = arr.filled(nd)

# create the output raster
out_fn = 'max_values.tif'
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(out_fn, xsize, ysize, 1, gdal.GDT_Float32)  # you might want to change the pixel type
out_band = out_ds.GetRasterBand(1)

# save and close output file
del out_ds, out_band

Note: I am assuming all your rasters are in the same directory and that no other raster is located in that folder. The code above will take all the rasters in the same directory and get the max value. Furthermore, make sure all of your rasters share the same dimensions (i.e. number of rows and columns).

  • Thank you I am attempting to run this through an anaconda prompt and this is the error I am receiving... arr = np.stack(arrays) # 3D array (date by rows by columns) File "<stdin>", line 6 arr = np.stack(arrays) # 3D array (date by rows by columns) ^ SyntaxError: invalid syntax
    – Bill Sands
    Jun 26, 2019 at 18:17
  • Are you running line by line by copying and pasting into anaconda prompt? Create a .py file, run that file and tell me if you get the same error. Jun 26, 2019 at 19:18
  • thanks It runs but the output file doesn't look correct. my range of values is 32762-32766 with most of the area having no value. I don't know how to post screenshots but I will try.
    – Bill Sands
    Jun 26, 2019 at 19:45
  • Oh, my bad. I mixed up the signs in the mask. I edited the code and it should be working now. Jun 26, 2019 at 20:09
  • perfect thank you!!
    – Bill Sands
    Jun 26, 2019 at 20:33

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.