Cómo calcular el índice NDVI con R

En este artículo vamos a tratar de mostrar cómo calcular el NDVI con R paso a paso. El cálculo del índice NDVI de vegetación mediante teledetección resulta muy útil para múltiples campos por sus múltiples aplicaciones.

Como veremos, el hecho de realizar el cálculo del índice de vegetación mediante el lenguaje de programación R supone una gran ventaja por la rapidez de la ejecución, así como por la sencillez del lenguaje.

El proceso lo realizaremos en base a una serie de comandos de R que comentaremos en detalle. Asimismo, para calcular el NDVI deberemos utilizar múltiples paquetes (packages) para este lenguaje de programación. Como base de cálculo usaremos imágenes de satélite provenientes del Sentinel-2.

Antes de ello, daremos unas pincelas teóricas para ponernos en contexto y comprender el procedimiento utilizado para el cálculo del índice NDVI mediante scripts en R.

¿Qué es el índice NDVI?

El NDVI se considera un índice de vegetación. Un índice de vegetación es un valor único que cuantifica la salud o estructura de la vegetación.

Las matemáticas asociadas al cálculo de un índice de vegetación se derivan de la física de la reflexión y la absorción de luz a través de las bandas espectrales.

Por ejemplo, se sabe que la vegetación saludable refleja la luz con fuerza en la banda del NIR (infrarrojo cercano, en inglés near Infrared) y mucho menos en la parte visible del espectro.

De este modo, si se crea una relación entre la luz reflejada en el NIR y la luz reflejada en la banda del rojo, se identificarán áreas que potencialmente tengan una vegetación saludable.

NDVI, vegetación y reflectancia

El NDVI (Normalized Difference Vegetation Index) es un índice cuantitativo de verdor que va desde -1 a 1, donde -1 representa un verdor mínimo o inexistente y 1 representa el verdor máximo.

Este índice se utiliza a menudo para una medida aproximada de salud, cobertura y fenología de la vegetación (etapa del ciclo de vida) en grandes extensiones territoriales.

También es habitual su uso en agricultura de precisión para estimar la salud de las plantaciones en determinadas parcelas.

El cálculo del NDVI se realiza a partir del canal visible y el NIR reflejado por la vegetación.

• La vegetación saludable (izquierda de la imagen) absorbe la mayor parte de la luz visible, y refleja una gran cantidad de energía de infrarrojo cercano.

• La vegetación poco saludable o escasa (derecha de la imagen) refleja una luz más visible y una luz menos infrarroja.

Puedes ampliar información en el artículo de la NASA sobre el NDVI, el cálculo de este índice de vegetación y sus aplicaciones. También puede encontrar información relevante en la web de EOS – Earth Observing System, desde donde visualizar y descargar imágenes satélite y NDVI.

Fórmula para el cálculo del NDVI e interpretación

La fórmula para el cálculo del índice NDVI es la siguiente:

Donde,

• NIR es la reflectancia del Infrarrojo Cercano (b8 Sentinel, b5 en Landsat 8, b2 MODIS).
• RED es la reflectancia del rojo (b4 Sentinel, b4 en Landsat 8, b1 MODIS).
• El rango de valores de salida estarán entre -1 y 1.

Utilizar R para calcular el NDVI

El lenguaje R es uno de los lenguajes de programación más usados en GIS para geoprocesos (junto con Python) que se encuentra en auge desde hace unos años por sus cualidades.

R permite realizar cálculos estadísticos y espaciales con capas tanto vectoriales como ráster cómodamente mediante scripts, haciendo uso de distintos paquetes específicos para tal fin.

El tratamiento de imágenes de satélite con R se convierte en una tarea relativamente fácil de realizar para usuarios con ciertos conocimientos. Bajo determinadas circunstancias, el uso de R puede ofrecer más ventajas que inconvenientes.

Entre dichas ventajas cabe mencionar, por ejemplo:

• la velocidad de cálculo en los geoprocesos y análisis espaciales
• la capacidad de ejecutar tareas de geoproceso sin una interfaz gráfica explícita
• la posibilidad de replicar la tarea de cálculo para múltiples capas o imágenes utilizando el mismo script.

Paquetes de R para calcular el NDVI

En el siguiente ejemplo vamos a calcular el NDVI en un área agrícola y forestal típicamente mediterránea. Para realizar tal ejercicio, necesitaremos una serie de librerías o paquetes de R que nos permitan leer, calcular y visualizar nuestros resultados.

En este ejemplo usaremos las librerías (packages):

• RStoolbox, específica para el tratamiento de imágenes satélite.
• raster, necesaria para el tratamiento de datos grid.
• ggplot2, muy conocida por su gran versatilidad a la hora de visualizar datos.
• gridExtra, que finalmente nos permitirá mostrar dos imágenes en una sola visualización.

Recuerda que si nunca has usado estas librerías, primero tendrás que instalarlas con el comando install.packages(“nombre_librería”). Si ya las tienes instaladas, simplemente usando el comando library() cargarás la librería para su uso. En el script que mostramos a continuación aparecen comentadas.

Análisis del código de R para calcular el NDVI

Fíjate que el código es sencillo. Primeramente, usamos list.files() para listar todas las bandas ráster contenidas en nuestra carpeta de trabajo, que hemos configurado previamente con la función setwd.

#install.packages("RStoolbox")
#install.package("raster")
#install.packages("gridExtra")
#install.packages("ggplot2")
library(RStoolbox)
library(raster)
library(ggplot2)
library(gridExtra)


# Configuramos nuestra ruta de trabajo (Carpeta con las bandas Sentinel) --

setwd("Ruta/de/trabajo")

# Lectura y apilado de las bandas del sensor Sentinel 2 -------------------

dir <- list.files(pattern = "S2",
                  full.names = T)

sentinel2 <- stack(dir)
names(sentinel2)<- c("Blue", "Green", "Red", "NIR")

Una vez están listadas, apilamos estas bandas en un solo objeto RasterStack, mediante la función stack de la librería raster. Para el cálculo del NDVI, simplemente, replicamos la ecuación descrita al principio del post, pero con las bandas Sentinel.

# Cálculo del NDVI --------------------------------------------------------


nir <- sentinel2$NIR # Infrarojo cercano
red <- sentinel2$Red # Rojo

ndvi <- (nir-red)/(nir+red) # NDVI

ndvi[ndvi>1] <- 1; ndvi[ndvi< (-1)] <- -1 #Reescalado para evitar outliers


# Visualización del NDVI e imagen falso color -----------------------------

# Visualización NDVI guardada en el objeto ndvi_plot

ndvi_plot <- ggR(ndvi, geom_raster = TRUE,alpha = 1)+
  scale_fill_gradientn(colours = rev(terrain.colors(100)), 
                       name = "NDVI") + 
  theme(legend.positio = "bottom")

Por último, lo que hacemos es crear los “plots” representando el NDVI y la combinación de falso color infrarrojo cercano – verde – azul.

# Visualización falso color guardada en el objecto falso_color

falso_color <- ggRGB(img_sentinel2, r= "NIR" , g="Green" , b="Blue",
                    stretch = "lin")

# Representación final

grid.arrange(ndvi_plot, falso_color, ncol=2)

Dicha representación se lleva a cabo con la función ggR y ggRGB, ambas contenidas en la librería RStoolbox y que se alimentan directamente de ggplot2 para realizar la visualización de los resultados. La función grid.arrange nos permite visualizar ambos mapas en una sola vista.

Interpretación del NDVI calculado con R

En la imagen resultado obtenida, se percibe claramente, que los valores más altos del índice hacen referencia a los campos de cultivo y los bosques, mientras que alrededor de 0 encontramos el tejido urbano y, por debajo, las masas de agua.

Sin embargo, también podemos detectar que los campos de cultivo están más sanos (mayor verdor) que los bosques.

Los datos usados en este ejercicio para calcular el NDVI se pueden descargar haciendo clic aquí.

Ante cualquier duda, deja un comentario a continuación en este post y lo responderemos lo antes posible.