Home - News ( United States | United Kingdom | Italy | Germany ) - Football scores

Showing content from http://ceholden.github.io/open-geo-tutorial/python/chapter_2_indices.html below:

chapter_2_indices

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 [6]:

# 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 range(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.

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 [4]:

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.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]]
0.904601300891

In Python 2 an integer divided by an integer produces an integer, even if the division would have produced a float point number. Python 3 changed this behavior, but if we run the NDVI calculation with Python 2 we would end up with all of our NDVI values equal to 0 because our input image is an integer datatype (int16). See documentation for division in NumPy for more information.

While we don't necessarily need to change anything in Python 3, it is generally useful to be explicit with the datatypes involved in our calculations for the sake of clarity. Additionally, we generally also want code written using Python 3 to work with Python 2.

In order to ensure we perform our calculation using floating point numbers, we can cast the demoninator or numerator of the calculation to a float:

In [13]:

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.9046013008913515
Mean NDVI: 0.7088133953809207
Median NDVI: 0.7319195214790647
Min NDVI: 0.09470304975922954

This looks correct.

Speaking of looking correct, the next chapter (link to webpage or Notebook) will demonstrate how you can visualize your results using actual plots!

RetroSearch is an open source project built by @garambo | Open a GitHub Issue

Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo

HTML: 3.2 | Encoding: UTF-8 | Version: 0.7.4