Chris Holden (ceholden@gmail.com) - https://github.com/ceholden

Chapter 2: Your first remote sensing vegetation index

Introduction

Now that we can read our data into the computer, let's calculate some vegetation indices.

The Normalized Difference Vegetation Index (NDVI) is so ubiquitous that it even has a Wikipedia entry. If you're here for learning how to do remote sensing image processing using GDAL and Python, I suspect you don't need any introduction to this section. If you need a refresher, please visit the Wikipedia URL for NDVI.

This chapter will be very straightfoward. We've already seen how we can read our imagery into a NumPy array -- this chapter will simply extend these tools by showing how to do simple calculations on NumPy objects.

Let's bring up our previous code for opening our image and reading in the data:



In [18]:

    
# Import the Python 3 print function
from __future__ import print_function

# Import the "gdal" and "gdal_array" submodules from within the "osgeo" module
from osgeo import gdal
from osgeo import gdal_array

# Import the NumPy module
import numpy as np

# Open a GDAL dataset
dataset = gdal.Open('../example/LE70220491999322EDC01_stack.gtif', gdal.GA_ReadOnly)

# Allocate our array using the first band's datatype
image_datatype = dataset.GetRasterBand(1).DataType

image = np.zeros((dataset.RasterYSize, dataset.RasterXSize, dataset.RasterCount),
                 dtype=gdal_array.GDALTypeCodeToNumericTypeCode(image_datatype))

# Loop over all bands in dataset
for b in xrange(dataset.RasterCount):
    # Remember, GDAL index is on 1, but Python is on 0 -- so we add 1 for our GDAL calls
    band = dataset.GetRasterBand(b + 1)
    
    # Read in the band's data into the third dimension of our array
    image[:, :, b] = band.ReadAsArray()
    

print('Red band mean: {r}'.format(r=image[:, :, 2].mean()))
print('NIR band mean: {nir}'.format(nir=image[:, :, 3].mean()))









    



Red band mean: 589.379808
NIR band mean: 3442.297712

Even from simple mean statistics over the entire image we can see the contrast between the red and the near-infared (NIR) bands.

NDVI

To calculate NDVI, we can simply use standard arithmetic operators in Python because these operations in NumPy are vectorized. Just like MATLAB, R, and other higher level languages, NEVER loop over a NumPy array unless you can avoid it.



In [27]:

    
b_red = 2
b_nir = 3

ndvi = (image[:, :, b_nir] - image[:, :, b_red]) / (image[:, :, b_red] + image[:, :, b_nir])

print(ndvi)
print(ndvi.max())









    



[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 ..., 
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]]
0

Oooops?

From looking at the printout of our NDVI result, we can see that all of the values are equal to zero.

NumPy is stupid then, right?

Well, we did tell NumPy that our input arrays were to be of an integer datatype. As mentioned in the documentation for division in NumPy, an integer array divided by another integer will result in an integer. As far as I'm concerned, this behavior is more respectful of my design decisions because I have very good control of when and where datatypes are converted. Just something to remember.

To fix this, let's cast one as a float:



In [35]:

    
ndvi = (image[:, :, b_nir] - image[:, :, b_red]) / \
        (image[:, :, b_nir] + image[:, :, b_red]).astype(np.float64)

print('NDVI matrix: ')
print(ndvi)

print('\nMax NDVI: {m}'.format(m=ndvi.max()))
print('Mean NDVI: {m}'.format(m=ndvi.mean()))
print('Median NDVI: {m}'.format(m=np.median(ndvi)))
print('Min NDVI: {m}'.format(m=ndvi.min()))









    



NDVI matrix: 
[[ 0.71390828  0.71079741  0.69352291 ...,  0.79392185  0.81408451
   0.79165379]
 [ 0.68064263  0.6787194   0.6643924  ...,  0.81387182  0.79880597
   0.77389811]
 [ 0.66904762  0.67268446  0.66332892 ...,  0.78495923  0.78278801
   0.81253291]
 ..., 
 [ 0.68301262  0.68593651  0.67145614 ...,  0.81065089  0.78050922
   0.76519266]
 [ 0.67341718  0.6622986   0.65331611 ...,  0.80436681  0.77483099  0.75      ]
 [ 0.63973799  0.62396514  0.66731813 ...,  0.7094648   0.70005244
   0.74574523]]

Max NDVI: 0.904601300891
Mean NDVI: 0.708813395381
Median NDVI: 0.731919521479
Min NDVI: 0.0947030497592

This looks correct.

Speaking of looking correct, the next chapter will demonstrate how you can visualize your results using actual plots!