Introducción

En este trabajo, vamos a testar diferentes metodologías de visualización de archivos raster satel·litales y a estudiar un indicador muy importantes para conocer el estado de la cobertura vegetal -muy utilizado en el seguimiento de las sequías y la prevención de incendios: el NDVI-, detallando los pasos para extraerlo de los archivos de imágenes satelitales que arcGis obtiene de Landsat-8.

Qué es el NDVI

El Índice de vegetación de diferencia normalizada, también conocido como NDVI por sus siglas en inglés, es el índice de vegetación más utilizado, empleado para estimar la cantidad, calidad y desarrollo de la vegetación con base a la medición, por medio de sensores remotos instalados comúnmente desde una plataforma espacial, de la intensidad de la radiación de ciertas bandas del espectro electromagnético que la vegetación emite o refleja

Este índice se basa en cómo refleja la energía y la luz. Usa las bandas infrarrojas cercana y roja del espectro electromagnético para estimar un indicador adimensional entre -1 y 1.

Una planta sana -con mucha clorofila y estructuras celulares- absorbe activamente la luz roja y refleja NIR (Near Infrared = Infrarojo Cercano) cuando ocurre la fotosíntesis. La planta se desarrolla y crece y contiene más estructuras celulares. Con una planta no saludable ocurre exactamente lo contrario.

El índice NDVI detecta y cuantifica la presencia de vegetación verde viva utilizando esta luz reflejada en las bandas visible e infrarroja cercana.

En cuanto a cómo se calcula el índice NDVI, se utiliza la fórmula estándar de NDVI:

\[NDVI=\frac{(NIR-red)}{NIR+red}\]

El rango del NDVI va desde -1 hasta 1 y nos muestra el vigor vegetativo.

Valores cercanos a 1: cuanto más intenso es el verde, más vigorosidad existe en la vegetación y la cubierta vegetal.

Valores cercanos a 0: corresponden a zonas con muy poca vegetación, primeras fases del cultivo, suelos desnudos o zonas no productivas.

Valores negativos: suele estar asociado a suelo desnudo, zonas de agua, nieve o nubes.

Así que podemos establecer la clasificación siguiente:

Principales librería de R que hay que instalar

En el ámbito de este trabajo vamos a utililizar: sf, para tratar datos vectoriales, arcgislayers, para obtención de datos de la plataforma arcGis i tratamiento de rasters, terra, para manipular i visualizar archivos raster, ggplot2 de tidyverse, para visualización i tmap, para obtener visualizaciones interactivas de varias capas.

# install.packages("wdpar", repos = "https://cran.rstudio.com/")

libs <- c(
    "wdpar", # Para obtener datos cartográficos de áreas protegidas
    "arcgislayers",
    "tidyverse",
    "sf",
    "terra",
    "classInt", # genera intervalos para las leyendas
    "tmap"
)

invisible(
    lapply(
        libs, library,
        character.only = TRUE
    )
)
terraOptions(progress=0) # Per a no veure la barra de progrés

Definición de la zona a mapear

En nuestro caso serà el área protegida del Cap de Creus, por su especial valor medioambiental i su vulnerabilidad a los incendios. Los datos del perimetro se pueden obtener en Protected Planet, una plataforme que diapone de datos cartográficos de zonas protegidas de todo el planeta.

# POLIGON DEL CAP DE CREUS
# https://www.protectedplanet.net/349127

temp <- tempfile()

unzip(zipfile = 
        '../repo/WDPA_WDOECM_Nov2024_Public_349127_shp/WDPA_WDOECM_Nov2024_Public_349127_shp_0.zip', # Poligon descarregat Cap de Creus
      exdir = temp)

capdc_sf <- sf::read_sf(temp)

capdc_sf <- st_transform(capdc_sf, crs = 4326)

terra::plot(
    sf::st_geometry(
        capdc_sf, add=TRUE)
    )
mtext("Protected Planet: The World Database on Protected Areas (WDPA), Cambridge, UK: UNEP-WCMC and IUCN. \nAvailable at: www.protectedplanet.net.",
      side=1, line=1, , adj=0.6, cex=0.7)

capdc_bbox <- sf::st_bbox(capdc_sf)

Obtención de los datos de Landsat 8 a través de arcGis

# 3. IMAGEN LANDSAT
#-----------------
# Landsat8 és el que disposa d'imatges més al dia (cada 16 dies)

url2 <- "https://landsat2.arcgis.com/arcgis/rest/services/Landsat8_Views/ImageServer"

landsat_data <- arcgislayers::arc_open(
    url2
)

# Per a extrerure la data d'adquisició
fecha <- landsat_data[[9]][3]$timeExtent[2]
fecha <- format(as.POSIXct(fecha / 1000, 
                  origin = "1970-01-01", tz = "UTC"), 
       "%Y-%m-%d %H:%M")

capdc_raster <- arcgislayers::arc_raster(
    x = landsat_data,
    xmin = capdc_bbox[["xmin"]],
    xmax = capdc_bbox[["xmax"]],
    ymin = capdc_bbox[["ymin"]],
    ymax = capdc_bbox[["ymax"]],
    crs = sf::st_crs(capdc_sf),
    width = 800,  # El valor per defecte és 400 pixels que és una mica petit
    height = 800  # si dona error, tornar al valor original
)

plot(capdc_raster$NearInfrared)
title(
  main = paste0("Imagen aproximada del raster \nobtenida el ",fecha),
  outer = TRUE, line = -1.5, cex.main=1)
mtext('Information System. QGIS Association. http://www.qgis.org',
      side=1, line=1, , adj=0.6, cex=0.8)

QGIS.org, 2024. QGIS Geographic Information System. QGIS Association. http://www.qgis.org

Cálculo del NDVI y su primera visualización

Se utiliza la fórmula que hemos visto arriba: \[NDVI=\frac{(NIR-red)}{NIR+red}\]

Para ello se van a utilizar los valores de la banda Infrarojo cercano y de la banda Infrarojo, que como se puede ver, son la cuarta i la quinta variable del raster obtenido:

names(capdc_raster)
##  [1] "CoastalAerosol"     "Blue"               "Green"             
##  [4] "Red"                "NearInfrared"       "ShortWaveInfrared1"
##  [7] "ShortWaveInfrared2" "Cirrus"             "QA"                
## [10] "ThermalInfrared1"   "ThermalInfrared2"
# Descripció
#  {"name": "NDVI Colorized", "description": 
# "Normalized difference vegetation index (NDVI) with color map. Dark green is thick vigorous vegetation and brown represents sparse vegetation.", "help": "" }, 
# { "name": "Normalized Difference Moisture Index Colorized", "description": "Normalized Difference Moisture Index with color map computed as (b5 - b6) / (b5 + b6). 
# Wetlands and moist areas are blues, and dry areas in deep yellow and brown.", "help": "" }

# Són per tant importants:
# Band 4: Red
# Band 5: NearInfrared
# 4. EL CÀLCUL DEL NDVI
#--------

# NDVI = (b5 - b4) / (b5 + b4)

# És un índex usat per estimar la quantitat, qualitat i desenvolupament 
# de la vegetació amb base al mesurament, per mitjà de sensors remots 
# instal·lats comunament des d'una plataforma espacial, de la intensitat 
# de la radiació de certes bandes de l'espectre electromagnètic que 
# la vegetació emet o reflecteix.

#El NDVI varia com a conseqüència entre -1,0 i +1,0. 


b5 <- capdc_raster[[5]]
b4 <- capdc_raster[[4]]

ndvi <- (b5 - b4) / (b5 + b4)

plot(ndvi)
title(
  main = "Imagen aproximada del NDVI",
  outer = TRUE, line = -1.5, cex.main=1)
mtext('Information System. QGIS Association. http://www.qgis.org',
      side=1, line=1, , adj=0.6, cex=0.8)

# PUNT DUBTÓS: LES FONTS CONSULTADES POSEN EL NDVI ENTRE -1 I +1.
# Teòricament no caldria utilitzar terra::clamp() perque els rangs són OK.

#ndvi_clamped <- terra::clamp(
#    x = ndvi,
#    lower = 0,
#    upper = 1,
#    values = TRUE
#)

#plot(ndvi_clamped)
#plot(sf::st_geometry(city_sf), add= TRUE)

Visualización con ggplot2

# 5. BREAKS AND COLORS
#---------------------

colors <- c('#A67B5C','#D0C59A','#6D9775','#104730')

#breaks <- c(-1,0,0.33,0.66,1)

breaks <- classInt::classIntervals( # Crea els breaks a partir dels valors
    terra::values(ndvi),             # basant-se en intervals naturals
    n = 4,
    fixedBreaks= c(-1,0,0.33,0.66,1),
    style = "fixed"
)$brks



#colors <- hcl.colors(
#    n = length(breaks),
#    palette = "Green-Brown",
#    rev = TRUE
#)

# 6. 2D MAP OF NDVI
#------------------

ndvi_df <- ndvi |>
    as.data.frame(xy = T) # xy indica que el df té long i lat i que hia ha 
                          # un tercer valor: l'NDVI

# head(ndvi_df)
names(ndvi_df)[3] <- "value"

annotation <- 'Information System. QGIS Association. http://www.qgis.org'

p <- ggplot() +
    geom_raster(
        data = ndvi_df,
        aes(
            x = x, y = y,
            fill = value)
        ) +
    scale_fill_gradientn(
        name = "NDVI",
        colors = colors,
        breaks = breaks,
        labels = round(breaks, 2)
    ) +
#    geom_text(
#      x = (max(ndvi_df$x)-min(ndvi_df$x))*.5,
#      y = min(ndvi_df$y),
#      label = 'Information System. QGIS Association. http://www.qgis.org',
#      size = 4.5 # font size
#    )+
#    guides(
#        fill = guide_colorbar(
#            direction = "horizontal",
#            barheight = unit(1.5 , "mm"),
#            barwidth = unit(100 , "mm"),
#            title.position = "top",
#            title.hjust = .5,
#            label.hjust = .5,
#            label.position = "bottom",
#            nrow = 1,
#            byrow = TRUE)
#        ) +
    #labs(caption = 'Information System. QGIS Association. http://www.qgis.org') +
    theme_minimal() +
    theme(
        axis.line = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        legend.position = "top",
        legend.title = element_text(
            size = 11, color = "grey10"
        ),
        legend.text = element_text(
            size = 8, color = "grey10"
        ),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        plot.background = element_blank(),
        legend.background = element_rect(
            fill = "white", color = NA
        ),
        plot.margin = margin(0, 0, 0, 0, "cm")
    )

library(grid)

grob <- grobTree(textGrob('Information System. QGIS Association. http://www.qgis.org', 
                          x=.3,  y=0.08, hjust=0,
  gp=gpar(col="black", fontsize=12, fontface="italic")))

p + annotation_custom(grob)

Visualización con tmap en modo ‘view’

Clave de interpretación de los valores de NDVI

  • -1 - 0: Vegetación muerta, suelo desnudo, o agua;
  • 0 - 0.33: Vegetación en sufrimiento;
  • 0.33-0.66: Vegetación medianamente sana;
  • 0.66- 1: Vegetación muy sana;
tmap_mode("view")

tm_shape(ndvi) +
  tm_raster(col = "NearInfrared",
            alpha = 0.75,
            style = "fixed",
            n= 4,
            breaks = breaks,
            palette = colors,
            legend.format = list(text.separator= '-'),
            title = "NDVI") +
tm_shape(capdc_sf)+
  tm_polygons(id= "NOMMUNI",
              alpha = 0,
popup.vars= NULL, # c()
#              popup.format = list(AREAM5000= list(big.mark = ".",
#                                                  decimal.mark= ",",
#                                                  digits= 1)),
              border.col = "steelblue",
              lwd= 2) +
tmap_options(basemaps = "Esri.WorldImagery",
               basemaps.alpha = 1)

QGIS.org, 2024. QGIS Geographic Information System. QGIS Association. http://www.qgis.org

Conclusión

Es factible utilizar imágenes satelitares del NDVI, aunque vemos que existe un lapso de tiempo y la publicación en arcGis.

Por otro lado, las visualizaciones con tmap nos parecen más útiles para visualizar este indicador, dando la posibilidad al ususario de acercar/alejar la imagen y de añadir y quitar capas, facilitándole así la localización topográfica de los valores.

Referencias