0

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?

2

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()

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.