El Blog colaborativo sobre Geotecnologías, GIS y Mapas para la comunidad GEO
  • Política de cookies
  • Más información
  • Contacto
Geotecnologías, GIS y Mapas
  • Proyecto
  • Contacto
Sin resultados
Ver todos los resultados
geomapik

Cómo calcular el índice NDVI con R

Marc Lemus por Marc Lemus
31 marzo, 2019
en Análisis GIS, Desarrollo GIS, Tutoriales
Home Análisis GIS
2.9k
Vistas
CompartirTwitter

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.

calcular NDVI con R
Esquema interpretativo de la reflectancia de la vegetación y su estado de verdor. Fuente: EOS.

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:

NDVI indice vegetación

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.

cómo calcular NDVI con R
Imagen resultado obtenida del cálculo del índice NDVI con R

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.

Etiquetas: AnálisisDesarrollo GISNDVIProgramaciónRTeledetecciónTutorial
Compartir92Tweet57
Post anterior

El plugin de GBIF en QGIS para descargar datos de distribución de especies

Post siguiente

Cómo obtener valores de una capa ráster en QGIS

Marc Lemus

Marc Lemus

Post siguiente
qgis extraer valores raster

Cómo obtener valores de una capa ráster en QGIS

interpolacion espacial qgis

Interpolación espacial en QGIS: métodos, procesos y evaluación

Comentarios 6

  1. Mariel Vanin says:
    1 año atrás

    Excelente!

    Responder
    • Raúl Estévez says:
      1 año atrás

      Gracias Mariel. Saludos!

      Responder
  2. Andres says:
    1 año atrás

    hola.. queria consultarle siguiente.. me da error:

    # 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)

    El erro que sale: Error in which(names(raster) == band) : object 'img_sentinel2' not found

    Me estaria faltando la img_sentinel2??

    Responder
    • Raúl Estévez says:
      1 año atrás

      Hola Andres,

      por el error que devuelve al ejectuar el script, es muy probable que sea ese el origen del problema. Revisa las variables especificadas y las rutas de acceso a la imagen de la que surge ese error.

      Un saludo.

      Responder
  3. Luis says:
    8 meses atrás

    Hola! Una pregunta, para qué es el #Reescalado para evitar outliers?
    También se aplica el reescalado si utilizó imágenes de Landsat 8?

    Responder
  4. gustavo says:
    6 meses atrás

    Gracias,
    al replicar tu ejemplo con otros raster, al hacer:
    sentinel2 <- stack(dir)
    me tira el error: # Error in compareRaster(rasters) : different extent

    Todos los raster tienen la misma ptoyeccion , ¿a que se debera ese error?

    gracias nuevamente

    Responder

Deja una respuesta Cancelar la respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Posts Recomendados

kepler gl mapbox mapas analisis datos

Kepler.gl de Mapbox: mapas interactivos y visualización de datos

13 abril, 2020
Cómo extraer y representar curvas de nivel con QGIS

Cómo extraer y representar curvas de nivel con QGIS

11 diciembre, 2018
Geotecnologías GIS Mapas Cartografía Digital

Geotecnologías, GIS y Mapas

Información, recursos de aprendizaje y novedades del sector GEO para los usuarios especializados en el manejo de la información geográfica, su análisis y su representación.

Posts Nuevos

que es la geoestadistica

¿Qué es la Geoestadística? Análisis geoestadístico con R, ArcGIS y QGIS

8 mayo, 2020
libros R Data Science GIS

Libros gratis de R para GIS y Análisis de datos espaciales

2 mayo, 2020
productos ESRI GIS

Principales productos de ESRI para GIS

26 abril, 2020
como cargar wms wmts qgis

Cómo cargar servicios WMS y WMTS en QGIS: proceso y diferencias

15 abril, 2020
bases de datos espaciales

¿Por qué trabajar con bases de datos espaciales en GIS?

12 abril, 2020
Geotecnologías GIS Mapas Cartografía Digital

El blog colaborativo sobre Geotecnologías, GIS y Mapas para la comunidad GEO

Información, recursos de aprendizaje y novedades del sector GEO para los usuarios especializados en el manejo de la información geográfica, su análisis y su representación.

Todo reunido en un blog colaborativo y abierto. ¿Te animas a colaborar con geomapik?

Categorías

  • Análisis GIS
  • Desarrollo GIS
  • Novedades GIS
  • Recursos GIS
  • Spatial Data Science
  • Tutoriales
  • Uncategorized
  • Webmapping
Sin resultados
Ver todos los resultados
Spatial Data Science

¿Qué es la Geoestadística? Análisis geoestadístico con R, ArcGIS y QGIS

8 mayo, 2020
Spatial Data Science

Libros gratis de R para GIS y Análisis de datos espaciales

2 mayo, 2020
Novedades GIS

Principales productos de ESRI para GIS

26 abril, 2020
Tutoriales

Cómo cargar servicios WMS y WMTS en QGIS: proceso y diferencias

15 abril, 2020

Etiquetas

Análisis (15) Análisis de datos (4) aprender GIS (5) ArcGIS (5) Bases de Datos Espaciales (5) Big Data (2) CARTO (2) Cartografía (27) catastro (2) Data Analysis (4) Data Science (4) Datos (9) Desarrollo GIS (15) Diseño cartográfico (6) ESRI (3) Geoinformática (7) Geomarketing (2) Geotecnología (4) GIS (37) Javascript (4) Jupyter Notebook (3) leaflet (2) Mapas (19) Mapbox (2) mapeo colaborativo (2) MDT (3) OpenLayers (2) OpenStreetMap (2) OSGEO (2) plugins (7) PostGIS (6) PostgreSQL (6) Programación (12) Python (7) QGIS (24) R (5) ráster (2) Sistemas de Información Geográfica (29) SQL (6) Teledetección (2) terreno (2) Tutorial (12) visor de mapas (3) visualización (3) Webmapping (11)

© 2018 geomapik - Geotecnologías, GIS y Mapas

Sin resultados
Ver todos los resultados
  • Proyecto
  • Contacto

© 2018 geomapik - Geotecnologías, GIS y Mapas

Este sitio web utiliza cookies para que usted tenga la mejor experiencia de usuario. Si continúa navegando está dando su consentimiento para la aceptación de las mencionadas cookies y la aceptación de nuestra política de cookies, pinche el enlace para mayor información.plugin cookies

ACEPTAR
Aviso de cookies