Lesson 2. Calculate NDVI Using NAIP Remote Sensing Data in the Python Programming Language


Learning Objectives

After completing this tutorial, you will be able to:

  • Calculate NDVI using NAIP multispectral imagery in Python.
  • Export or write a raster to a .tif file from Python.

Calculate NDVI in Python

Sometimes you can download already calculated NDVI data products from a data provider.

However, in this case, you don’t have a pre calculated NDVI product from NAIP data. You need to calculate NDVI using the NAIP imagery / reflectance data that you have downloaded from Earth Explorer.

How to Derive the NDVI Vegetation Index From Multispectral Imagery

The normalized difference vegetation index (NDVI) uses a ratio between near infrared and red light within the electromagnetic spectrum. To calculate NDVI you use the following formula where NIR is near infrared light and red represents red light. For your raster data, you will take the reflectance value in the red and near infrared bands to calculate the index.

(NIR - Red) / (NIR + Red)

You can perform this calculation using matrix math with the numpy library.

To get started, load all of the required Python libraries.

import os
import matplotlib.pyplot as plt
import numpy as np
import rasterio as rio
import geopandas as gpd
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep

# Download data and set working directory
data = et.data.get_data('cold-springs-fire')
os.chdir(os.path.join(et.io.HOME, 'earth-analytics'))
Downloading from https://ndownloader.figshare.com/files/10960109
Extracted output to /root/earth-analytics/data/cold-springs-fire/.

Next, open up the NAIP data that you wish to calculate NDVI with. You will use the data from 2015 for this example that you downloaded for week 7 of this course:

data/cold-springs-fire/naip/m_3910505_nw_13_1_20150919/crop/m_3910505_nw_13_1_20150919_crop.tif

naip_data_path = os.path.join("data", "cold-springs-fire", 
                              "naip", "m_3910505_nw_13_1_20150919", 
                              "crop", "m_3910505_nw_13_1_20150919_crop.tif")

with rio.open(naip_data_path) as src:
    naip_data = src.read()

# View shape of the data
naip_data.shape
(4, 2312, 4377)

Calculate NDVI using regular numpy array math. In this case, the bands you are subtracting come from the same data file. The tif file format requires that all layers are in the same CRS and of the same size so you assume the data line up. Thus you do not need to test the data for equal shape, crs and extent.

naip_ndvi = es.normalized_diff(naip_data[3], naip_data[0])

Finally plot the data. Note below that the vmin= and vmax= arguments are used to stretch the colorbar across the full possible range of NDVI values (-1 to 1).

ep.plot_bands(naip_ndvi, 
              cmap='PiYG',
              scale=False,
              vmin=-1, vmax=1,
              title="NAIP Derived NDVI\n 19 September 2015 - Cold Springs Fire, Colorado")
plt.show()
Plotting the NDVI calculation of the 2015 NAIP data with a colorbar that reflects the data.
Plotting the NDVI calculation of the 2015 NAIP data with a colorbar that reflects the data.

View Distribution of NDVI Values

Using a histogram, you can view the distribution of pixel values in your NDVI output.

ep.hist(naip_ndvi,
        figsize=(12, 6),
        title=["NDVI: Distribution of pixels\n NAIP 2015 Cold Springs fire site"])

plt.show()
Histogram of NDVI values derived from 2015 NAIP data.
Histogram of NDVI values derived from 2015 NAIP data.

Optional - Export a Numpy Array to a Raster Geotiff in Python

When you are done, you can export your NDVI raster data so you could use them in QGIS or ArcGIS or share them with your colleagues. To do this, you use the rio.write() function.

Exporting a raster in Python is a bit different from what you may have learned using another language like R. In Python, you need to:

  1. Create a new raster object with all of the metadata needed to define it. This metadata includes:
    • the shape (rows and columns) of the object
    • the coordinate reference system (crs)
    • the type of file (you will export a geotiff (.tif) in this lesson
    • and the type of data being stored (integer, float, etc).

Lucky for you, all of this information can be accessed from the original NAIP data that you imported into Python using attribute calls like:

.transform and .crs

To implement this, below you will create a rasterio object to grab the needed spatial attributes.

naip_data_path = os.path.join("data", "cold-springs-fire", 
                              "naip", "m_3910505_nw_13_1_20150919", 
                              "crop", "m_3910505_nw_13_1_20150919_crop.tif")

with rio.open(naip_data_path) as src:
    naip_data = src.read()
    naip_meta = src.profile

naip_meta
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32768.0, 'width': 4377, 'height': 2312, 'count': 4, 'crs': CRS.from_wkt('PROJCS["UTM Zone 13, Northern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'), 'transform': Affine(1.0, 0.0, 457163.0,
       0.0, -1.0, 4426952.0), 'tiled': False, 'compress': 'lzw', 'interleave': 'band'}
naip_transform = naip_meta["transform"]
naip_crs = naip_meta["crs"]

# View spatial attributes
naip_transform, naip_crs
(Affine(1.0, 0.0, 457163.0,
        0.0, -1.0, 4426952.0),
 CRS.from_wkt('PROJCS["UTM Zone 13, Northern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'))

You can view the type of data stored within the ndvi array using .dtype. Remember that the naip_ndvi object is a numpy array.

type(naip_ndvi), naip_ndvi.dtype
(numpy.ndarray, dtype('float64'))

Use rio.open() to create a new blank raster ‘template’. Then write the NDVI numpy array to to that template using dst.write().

Note that when we write the data we need the following elements:

  1. the driver or type of file that we want to write. ‘Gtiff’ is a geotiff format
  2. dtype: the structure of the data that you are writing. We are writing floating point values (values with decimal places)
  3. the heigth and width of the ndvi object (accessed using the .shape attribute)
  4. the crs of the spatial object (accessed using the rasterio NAIP data)
  5. the transform information (accessed using the rasterio NAIP data)

Finally you need to specify the name of the output file and the path to where it will be saved on your computer.

Export a Numpy Array to a Raster Geotiff Using the Spatial Profile or Metadata of Another Raster

You can use the naip_meta variable that you created above. This variable contains all of the spatial metadata for naip data.

In this case, the

  1. number of bands (we have only one band vs 4 in the color image) and
  2. the data format (we have floating point numbers - numers with decimals - vs integers)

have changed. Update those values then write out the image.

naip_meta
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32768.0, 'width': 4377, 'height': 2312, 'count': 4, 'crs': CRS.from_wkt('PROJCS["UTM Zone 13, Northern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'), 'transform': Affine(1.0, 0.0, 457163.0,
       0.0, -1.0, 4426952.0), 'tiled': False, 'compress': 'lzw', 'interleave': 'band'}
# Change the count or number of bands from 4 to 1
naip_meta['count'] = 1

# Change the data type to float rather than integer
naip_meta['dtype'] = "float64"
naip_meta
{'driver': 'GTiff', 'dtype': 'float64', 'nodata': -32768.0, 'width': 4377, 'height': 2312, 'count': 1, 'crs': CRS.from_wkt('PROJCS["UTM Zone 13, Northern Hemisphere",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]'), 'transform': Affine(1.0, 0.0, 457163.0,
       0.0, -1.0, 4426952.0), 'tiled': False, 'compress': 'lzw', 'interleave': 'band'}

Note below that when we write the raster, we use **naip_meta.

The two ** tells Python to unpack all of the values in the naip_meta object to use as arguments when writing the geotiff file. We already updated the elements that we needed to above (count and dtype). So this naip_meta object is ready to be used for the NDVI raster.

naip_ndvi_outpath = os.path.join("data", "cold-springs-fire", 
                                 "outputs", "naip_ndvi.tif")

# Write your the ndvi raster object
with rio.open(naip_ndvi_outpath, 'w', **naip_meta) as dst:
    dst.write(naip_ndvi, 1)

Leave a Comment