NDVI computation

This tutorial shows how to use PyEPR to open a MERIS L1B product, compute the Normalized Difference Vegetation Index (NDVI) and store it into a flat binary file.

The example code (examples/write_ndvi.py) is a direct translation of the C sample program write_ndvi.c bundled with the EPR API distribution.

The program is invoked as follows:

$ python write_ndvi.py <envisat-oroduct> <output-file>

The code have been kept very simple and it consists in a single function (main()) that also performs a minimal command line arguments handling.

#!/usr/bin/env python3

# This program is a direct translation of the sample program
# "write_ndvi.c" bundled with the EPR-API distribution.
#
# Source code of the C program is available at:
# https://github.com/bcdev/epr-api/blob/master/src/examples/write_ndvi.c


"""Example for using the epr-api

Demonstrates how to open a MERIS L1b product and calculate the NDVI.

This example does not demonstrate how to write good and safe code.
It is reduced to the essentials for working with the epr-api.

Calling sequence::

    $ python3 write_ndvi.py <envisat-product> <output-file>

for example::

    $ python3 write_ndvi.py MER_RR__1P_test.N1 my_ndvi.raw

"""

import sys
import struct
import logging

import epr


def main(*argv):
    if not argv:
        argv = sys.argv

    if len(argv) != 3:
        print("Usage: write_ndvi <envisat-product> <output-file>")
        print("  where envisat-product is the input filename")
        print("  and output-file is the output filename.")
        print("Example: MER_RR__1P_TEST.N1 my_ndvi.raw")
        print
        sys.exit(1)

    # Open the product
    with epr.open(argv[1]) as product:

        # The NDVI shall be calculated using bands 6 and 8.
        band1_name = "radiance_6"
        band2_name = "radiance_10"

        band1 = product.get_band(band1_name)
        band2 = product.get_band(band2_name)

        # Allocate memory for the rasters
        width = product.get_scene_width()
        height = product.get_scene_height()
        subsampling_x = 1
        subsampling_y = 1
        raster1 = band1.create_compatible_raster(width, height,
                                                 subsampling_x, subsampling_y)
        raster2 = band2.create_compatible_raster(width, height,
                                                 subsampling_x, subsampling_y)

        # Read the radiance into the raster.
        offset_x = 0
        offset_y = 0

        logging.info(f"read {band1_name!r} data")
        band1.read_raster(offset_x, offset_y, raster1)

        logging.info(f"read {band2_name!r} data")
        band2.read_raster(offset_x, offset_y, raster2)

        # Open the output file
        logging.info(f"write ndvi to {argv[2]!r}")
        with open(argv[2], "wb") as out_stream:

            # Loop over all pixel and calculate the NDVI.
            #
            # @NOTE: looping over data matrices is not the best soluton.
            #        It is done here just for demostrative purposes
            for j in range(height):
                for i in range(width):
                    rad1 = raster1.get_pixel(i, j)
                    rad2 = raster2.get_pixel(i, j)
                    if (rad1 + rad2) != 0.0:
                        ndvi = (rad2 - rad1) / (rad2 + rad1)
                    else:
                        ndvi = -1.0
                    out_stream.write(struct.pack("f", ndvi))
            logging.info("ndvi was written success")


if __name__ == "__main__":
    main()

The ENVISAT epr.Product is opened using the epr.open() function.

    # Open the product
    with epr.open(argv[1]) as product:

As usual in modern python programs the with statement has been used to ensure that the product is automatically closed as soon as the program exits the block. Of course it is possible to use a simple assignment form:

product = open(argv[1])

but in this case the user should take care of manually call:

product.close()

when appropriate.

The name of the product is in the first argument passed to the program. In order to keep the code simple no check is performed to ensure that the product is a valid L1B product.

The NDVI is calculated using bands 6 and 8 (the names of these bands are “radiance_6” and “radiance_10”). epr.Band objects are retrieved using the epr.Product.get_band() method:

        # The NDVI shall be calculated using bands 6 and 8.
        band1_name = "radiance_6"
        band2_name = "radiance_10"

        band1 = product.get_band(band1_name)
        band2 = product.get_band(band2_name)

band1 and band2 are used to read the calibrated radiances into the epr.Raster objects that allow to access data matrices with the radiance values.

Before reading data into the epr.Raster objects they have to be instantiated specifying their size and data type in order to allow the library to allocate the correct amount of memory.

For sake of simplicity epr.Raster object are created with the same size of the whole product (with no sub-sampling) using the epr.Band.create_compatible_raster() method of the epr.Band class:

        # Allocate memory for the rasters
        width = product.get_scene_width()
        height = product.get_scene_height()
        subsampling_x = 1
        subsampling_y = 1
        raster1 = band1.create_compatible_raster(width, height,
                                                 subsampling_x, subsampling_y)
        raster2 = band2.create_compatible_raster(width, height,
                                                 subsampling_x, subsampling_y)

Then data are actually loaded into memory using the epr.Band.read_raster() method. Since epr.Raster objects have been defined to match the whole product, offset parameters are set to zero (data are read starting from specified offset):

        # Read the radiance into the raster.
        offset_x = 0
        offset_y = 0

        logging.info(f"read {band1_name!r} data")
        band1.read_raster(offset_x, offset_y, raster1)

        logging.info(f"read {band2_name!r} data")
        band2.read_raster(offset_x, offset_y, raster2)

Note

in this simplified example it is assumed that there is enough system memory to hold the two epr.Raster objects.

After opening (in binary mode) the stream for the output

        # Open the output file
        logging.info(f"write ndvi to {argv[2]!r}")
        with open(argv[2], "wb") as out_stream:

the program simply loops over all pixel and calculate the NDVI with the following formula:

NDVI = \frac{radiance_{10} - radiance_8}{radiance_{10} + radiance_8}

            # Loop over all pixel and calculate the NDVI.
            #
            # @NOTE: looping over data matrices is not the best soluton.
            #        It is done here just for demostrative purposes
            for j in range(height):
                for i in range(width):
                    rad1 = raster1.get_pixel(i, j)
                    rad2 = raster2.get_pixel(i, j)
                    if (rad1 + rad2) != 0.0:
                        ndvi = (rad2 - rad1) / (rad2 + rad1)
                    else:
                        ndvi = -1.0
                    out_stream.write(struct.pack("f", ndvi))
            logging.info("ndvi was written success")

This part of the code tries to mimic closely the original C code (write_ndvi.c)

out_stream = fopen(argv[2], "wb");
for (j = 0; j < height; ++j) {
    for (i = 0; i < width; ++i) {
        rad1 = epr_get_pixel_as_float(raster1, i, j);
        rad2 = epr_get_pixel_as_float(raster2, i, j);
        if ((rad1 + rad2) != 0.0) {
            ndvi = (rad2 - rad1) / (rad2 + rad1);
        } else {
            ndvi = -1.0;
        }
        status = fwrite( & ndvi, sizeof(float), 1, out_stream);
    }
}
epr_log_message(e_log_info, "ndvi was written success");

and uses the epr.Raster.get_pixel() method to access pixel values and perform computation.

The Python struct.pack() function together with file.write() is used to write the NDVI of the pixel n the file in binary format.

                    out_stream.write(struct.pack("f", ndvi))

Note

the entire solution is quite not “pythonic”. As an alternative implementation it could be used the numpy.ndarray interface of epr.Raster objects available via the epr.Raster.data property. The NDVI index is computed on all pixels altogether using vectorized expressions:

# Initialize the entire matrix to -1
ndvi = numpy.zeros((height, width), 'float32') - 1

aux = raster2.data + raster1.data

# indexes of pixel with non null denominator
idx = numpy.where(aux != 0)

# actual NDVI computation
ndvi[idx] = (raster2.data[idx] - raster1.data[idx]) / aux[idx]

Finally data can be saved to file simply using the numpy.ndarray.tofile() method:

ndvi.tofile(out_stream)