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.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
#!/usr/bin/env python

# 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::

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

for example::

    $ python write_ndvi.py MER_RR__1P_test.N1 my_ndvi.raw

'''

from __future__ import pront_function

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
    product = epr.open(argv[1])

    # 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('read "%s" data' % band1_name)
    band1.read_raster(offset_x, offset_y, raster1)

    logging.info('read "%s" data' % band2_name)
    band2.read_raster(offset_x, offset_y, raster2)

    # Open the output file
    logging.info('write ndvi to "%s"' % argv[2])
    out_stream = open(argv[2], 'wb')

    # 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')

    out_stream.close()
    product.close()


if __name__ == '__main__':
    main()

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

The name of the product is in the first argument of the program call. 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:

    product = epr.open(argv[1])

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

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:

    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)

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

                                             subsampling_x, subsampling_y)

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

    logging.info('read "%s" data' % band1_name)
    band1.read_raster(offset_x, offset_y, raster1)

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

    band2.read_raster(offset_x, offset_y, raster2)

    # Open the output file

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}

    out_stream = open(argv[2], 'wb')

    # 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

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.

            else:

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)

Previous topic

Exporting bitmasks

Next topic

API Reference

This Page