# Get maximum pixel value within numeric range

I 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?

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
os.chdir(r'path/to/your/rasters')

# 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)
arrays.append(arr)
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)

# 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_ds.SetGeoTransform(gt)
out_ds.SetProjection(sr)
out_band = out_ds.GetRasterBand(1)
out_band.SetNoDataValue(nd)
out_band.WriteArray(arr)

# save and close output file
out_band.FlushCache()
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 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. 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!! Jun 26, 2019 at 20:33