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

Question

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?

Answer 1

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])
    scale_factor = float(data.GetMetadata()['scale_factor'])
    array = data.ReadAsArray() * scale_factor
    return array

def get_geoinfo(filepath, subdataset):
    sourceds = gdal.Open(filepath, gdal.GA_ReadOnly)
    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()