How to calculate the Min & Max value for a pixel from multiple Raster Images?

I have downloaded a Bunch of MODIS MOD13Q1.006 Images for my Study Area, covering a period of 4 years (same scene across time).

Now I want to calculate the Min & Max NDVI value for each Pixel, so that I can then calculate the VCI Value with future data.

Using GDAL, how do I calculate the Min & Max Value for each Pixel?

Using gdal in python, you could open the first file and use those values as your "base-values". Then looping over the rest of your files, comparing the base-values with the just opened values and each time selecting the maxes (or minuses). Something like this:

``````# Create a list of files, have a look at glob.glob() to list many files.
files = [r"D:\some_path\MOD11A1.A2005001.h28v07.005.2007347211605.hdf",
r"D:\some_path\MOD11A1.A2005002.h28v07.005.2007348210508.hdf",
r"D:\some_path\MOD11A1.A2005003.h28v07.005.2007348221535.hdf",
r"D:\some_path\MOD11A1.A2005004.h28v07.005.2007349081734.hdf",
]

# Select the first file.
filepath = files[0]

# Get some general info from that file.
ndv, xsize, ysize, geot, projection = get_geoinfo(files[0], 0)

# Open the file twice to get your "base-values".
max_values = open_hdf_file(files[0], 0).astype(np.float)
min_values = open_hdf_file(files[0], 0).astype(np.float)

# Set the NDV pixels to NaN values.
max_values[max_values == ndv] = np.nan
min_values[min_values == ndv] = np.nan

# Loop over the rest of your files.
for filepath in files[1:]:

next_array = open_hdf_file(filepath, 0).astype(np.float)
next_array[next_array == ndv] = np.nan

max_values = np.nanmax(np.array([max_values, next_array]), axis = 0)
min_values = np.nanmin(np.array([min_values, next_array]), axis = 0)

# Save the created arrays as, for example, geotiffs.
out_filepath_max = r"C:\some_path\max_values.tif"
out_filepath_min = r"C:\some_path\min_values.tif"

geo_info = (ndv, xsize, ysize, geot, projection)

create_geotif(out_filepath_max, max_values, geo_info)
create_geotif(out_filepath_min, min_values, geo_info)
``````

Using these function:

``````import gdal, osr
import numpy as np

def open_hdf_file(filepath, subdataset):
hdf_file = gdal.Open(filepath)
data = gdal.Open(hdf_file.GetSubDatasets()[subdataset][0])
return array

def get_geoinfo(filepath, subdataset):
sourceds = gdal.Open(sourceds.GetSubDatasets()[subdataset][0])
ndv = sourceds.GetRasterBand(1).GetNoDataValue()
xsize = sourceds.RasterXSize
ysize = sourceds.RasterYSize
geot = sourceds.GetGeoTransform()
projection = osr.SpatialReference()
projection.ImportFromWkt(sourceds.GetProjectionRef())
return ndv, xsize, ysize, geot, projection

def create_geotif(out_filepath, array, geo_info):
datatypes = {gdal.GetDataTypeName(i).lower() : i for i in range(1, 12)}
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(out_filepath, geo_info[1], geo_info[2], 1, datatypes[array.dtype.name])
array[np.isnan(array)] = geo_info[0]
dataset.GetRasterBand(1).SetNoDataValue(geo_info[0])
dataset.SetGeoTransform(geo_info[3])
dataset.SetProjection(geo_info[4].ExportToWkt())
dataset.GetRasterBand(1).WriteArray(array)
dataset = None
``````

Note that:

• This script would only work with the images from the same scene.

• The second input in the functions get_geoinfo and open_hdf_file (all zero in the example) is used to choose which subdataset to use from the hdf-files. Have a look at the code below to make sure you select the right subdataset:

hdf_file = gdal.Open(filepath)
hdf_file.GetSubDatasets()