Cálculo del NDVI de una imagen Landsat8 con Python3 y Rasterio - Tutorial

El índice de vegetación NDVI es muy utilizado para la evaluación del medio ambiente, cultivos agrícolas y el cambio de uso de suelo. Este procedimiento es muy sencillo y popular en softwares de Sistemas de Información Geográfica (SIG) como QGIS, sin embargo el procedimiento es efectivo y eficiente cuando se trata de una sola imagen. Pero que pasaría si tenemos que analizar una serie de imágenes?, o si tenemos una imagen muy grande?, o si tenemos que aplicar algún filtro especial a la imagen, es allí donde nos planteamos la necesidad de analizar la imagen de una manera más abstracta y más eficiente con los recursos computacionales disponibles.

Las imágenes satelitales son georasters, es decir, son un arreglo regular de columnas y filas (matriz) con una georeferenciacion. Python es un lenguaje de programación y análisis muy versátil para el análisis de matrices con la librería Numpy, sin embargo hasta ahora no era muy eficiente para el procesamiento de georasters hasta la llegada del paquete Rasterio.

Rasterio es una librería de Python que permite la lectura, inspección, visualización y la escritura de raster geospaciales. La librería usa rasters en formato GeoTIFF y otros formatos y es capaz de trabajar con imágenes satelitales, modelos de elevación digital, productos de imágenes, e imágenes procesadas de drone. Rasterio permite importar un raster geoespacial de una banda y multibanda y permite trabajar en un entorno interactivo como un Jupyter notebook. La librería conserva la “dualidad” del georaster, eso significa que puede manejar los parámetros de localización y resolución del raster tanto como los valores matriciales de los elementos grillados.

Este tutorial muestra el procedimiento para el cálculo del NDVI de un imagen satelital Landsat 8 utilizando Python y Rasterio. El tutorial se realiza en un entorno interactivo llamado Jupyter Notebook. Para la instalación del software requerido para este ejercicio se recomienda seguir este tutorial.

Ventajas de usar Python para el análisis de raster

Por que cambiar de entorno de trabajo si el cálculo de NDVI se puede hacer fácilmente en QGIS? Esta es una buena pregunta para los analistas espaciales. El trabajo en Python con Rasterio es más abstracto, requiere un mayor dominio de los comandos y no es tan intuitivo, pero tiene las siguientes ventajas:

• Análisis masivo de rasters

• Automatización de procesos

• Menor tiempo de cálculo

• Ambiente con menos errores

Tutorial

Código

Este es el código en Python utilizado para este análisis.

#import required libraries
import rasterio
from rasterio import plot
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

import os
os.listdir('../Landsat8/')

#import bands as separate 1 band raster
imagePath = '../Landsat8/'
band4 = rasterio.open(imagePath+'LC08_L1TP_042035_20180603_20180615_01_T1_B4.TIF') #red
band5 = rasterio.open(imagePath+'LC08_L1TP_042035_20180603_20180615_01_T1_B5.TIF') #nir

#number of raster bands
band4.count

#number of raster columns
band4.width

#number of raster rows
band4.height

#plot band 
plot.show(band4)

#type of raster byte
band4.dtypes[0]

#raster sytem of reference
band4.crs

#raster transform parameters
band4.transform

#raster values as matrix array
band4.read(1)

#multiple band representation
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
plot.show(band4, ax=ax1, cmap='Blues')
plot.show(band5, ax=ax2, cmap='Blues')
fig.tight_layout()

#generate nir and red objects as arrays in float64 format
nir = band5.read(1).astype('float64')
red = band4.read(1).astype('float64')

nir[4000:4005,4000:4005]

#ndvi calculation, empty cells or nodata cells are reporte as 0
ndvi = np.where(
                (nir+red)==0.,
                0,
                (nir-red)/(nir+red)
               )

#export ndvi image
ndviImage = rasterio.open('../Output/ndviImage.tiff', 'w', driver='Gtiff',
                          width=band4.width, height=band4.height,
                          count=1,
                          crs=band4.crs,
                          transform=band4.transform,
                          dtype='float64'                  
                         )
ndviImage.write(ndvi,1) #ndvi
ndviImage.close()

#plot ndvi
ndvi = rasterio.open('../Output/ndviImage.tiff')
plot.show(ndvi, cmap='Greens')

Datos de ingreso

Una parte recortada de la imagen Landsat 8 y los códigos de procesamiento se pueden descargar de este enlace.