1 INTRODUCCIÓN

Las transformaciones ráster implican la manipulación matemática de imágenes para generar otras nuevas. Estas imágenes transformadas o bien contienen información de interés o bien resaltan las características singulares de las imágenes originales, aunque también se realizan para preservar la información esencial en un número reducido de dimensiones o capas.

Las herramientas más utilizadas para transformar las imágenes de satélite son las siguientes:

#library(devtools)
#install_github("bleutner/RStoolbox")

En primer lugar, deben importarse todos los paquetes necesarios para realizar este ejercicio:

library(raster)             # Importación y transformación de imágenes ráster
## Loading required package: sp
library(rasterVis)          # Visualización de datos ráster
## Loading required package: lattice
library(RStoolbox)          # Transformación y visualización de datos ráster.
library(ggplot2)            # Visualización de datos ráster con RStoolbox
library(sf)
## Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE

En segundo lugar, debe establecerse la carpeta de trabajo

setwd("D:/G174_OCW/Lab_5_transformaciones_matematicas")

A continuación se creará una lista para posteriormente agrupar todas las bandas en un único Rasterbrick.

lista <- list.files("./datos/originales/", 
                    pattern = ".*B[1234567]\\.tif$",            
                    ignore.case=TRUE,                           
                    full.names = TRUE)
stack<- stack(lista)
imagen <- brick(stack)

A partir de este momento, se prescinde de aquellos objetos que no serán necesarios.

rm(lista, stack)

Los nombres originales de las bandas son:

names(imagen)                                                         
## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"

Para facilitar el trabajo con las 7 bandas, se cambiarán los nombres de las bandas individuales en la imagen multibanda importada por los correspondientes a sus respectivas longitudes de onda (Landsat 8 OLI ).

names(imagen) <- c("coastal aerosol", "blue", "green", "red", "NIR", "SWIR1", "SWIR2")
names(imagen)                                                        
## [1] "coastal.aerosol" "blue"            "green"           "red"            
## [5] "NIR"             "SWIR1"           "SWIR2"

Como se ha realizado en otras ocasiones, es conveniente conocer los atributos de ese objeto raster

imagen                                  
## class      : RasterBrick 
## dimensions : 1245, 1497, 1863765, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : coastal.aerosol,          blue,         green,           red,           NIR,         SWIR1,         SWIR2 
## min values :    0.0964179114,  0.0748399049,  0.0425921641,  0.0208406653,  0.0008457669, -0.0078721829, -0.0050529451 
## max values :       0.7346282,     0.7177562,     0.6924697,     0.7861769,     1.0124315,     1.0432045,     1.1179360

También es conveniente una representación gráfica en falso color.

ggRGB(imagen, r=5, g=4, b=3, stretch ="lin", geom_raster = TRUE) + 
  ggtitle("Fichero original") +                                                 # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +                                      # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,                                   # Título centrado
                                  size =20),                                    # Tamaño del título
        axis.title = element_text(size =10))                                    # Tamaño de las etiquetas de los ejes

2 CÁLCULO DE ÍNDICES ESPECTRALES

Un índice espectral es una medida cuantitativa de ciertas características de las imágenes, generalmente formada por la combinación de varias bandas espectrales. Los valores de estas bandas se someten a ciertas operaciones matemáticas (suma, resta, multiplicación y división) que resultan en un único parámetro, conocido como índice espectral. Puede considerarse un indicador mejorado de alguna de las propiedades de las imágenes, y suministra información valiosa sobre características singulares, en comparación con la información original de las bandas.

Los índices espectrales mejoran la respuesta radiométrica de una característica particular y al mismo tiempo reducen la dimensión de los datos. Facilitan la comparación tanto entre diferentes características o como entre características similares presentes en diferentes ubicaciones o momentos en el tiempo. Estos índices también reducen el ruido causado por las variaciones en la iluminación, la influencia de la atmósfera etc…

Existen varios índices espectrales disponibles en la literatura con sus respectivas ventajas y limitaciones. Para ilustrar las diversas maneras de cálculo de estos índices en R se utilizará como ejemplo el índice NVDI. El Índice de vegetación de diferencia normalizada, también conocido como NDVI por sus siglas en inglés, es uno de los índices espectrales más famosos y utilizados. Con este índice se puede estimar la cantidad, calidad y desarrollo de la vegetación. La fórmula matemática que calcula este índice es la siguiente

Fórmula de cálculo del índice NDVI

En R existen diferentes posibilidades de cálculo.

2.1 Método 1: cálculo directo

Este método usa directamente las capas de un objeto ráster mediante sus índices [[]] o sus nombres (b1…). Los operandos típicos (p. ej., +, -, /, *) también se pueden usar como funciones (p. ej., sqrt, log, cos). Sin embargo, dado que el procesamiento ocurre a la vez en la memoria, debe asegurarse de que sus datos quepan en la RAM. Un ejemplo de cálculo empleando las capas entre índices es el siguiente

ndvi <- (imagen[[5]] - imagen[[4]]) / (imagen[[5]] + imagen[[4]])
ndvi
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : -0.959918, 0.883263  (min, max)

Una vez calculado el índice, se puede proceder a una representación gráfia de calidad por medio de la función ggR..

ggR(ndvi, geom_raster = TRUE) + 
  ggtitle("Índice NDVI") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10),         # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),               # size of the legend
        legend.title = element_text(size = 20),      #size of Legend tittle
        legend.text = element_text (size = 15)) +    # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores NDVI", colours = c("white", "darkgreen"))

También podemos verificar la distribución de los valores NDVI en el ráster trazando su histograma.

hist(ndvi, main = "Distribution of NDVI",              # Title of histograma
     xlab ="NDVI", ylab = "Frecuencia",                # Etiquetas de los ejes
     col ="sky blue",                                  # Colores de los bins
     xlim = c(-1,1), breaks = 20,                      # Limites e intervalos del eje x
     xaxt = "n")                                       # Deshabilitar el eje predeterminado
# Eje personalizado del histograma
axis(side = 1, at = seq(-1, 1, 0.1), labels = seq(-1, 1, 0.1))

El histograma muestra que la mayoría de los píxeles de la imagen ráster tienen valores NDVI que oscilan entre 0,2 y 0,8, si bien la mayoría se sitúa a partir de 0,1.

Como se ha señalado anteriormente, existen otras alternativas para este tipo de transformación. Por ejemplo, susitutir el número correspondiente de la banda por su nombre.

ndvi_2 <- (imagen[["NIR"]] - imagen[["red"]]) / (imagen[["NIR"]] + imagen[["red"]])
ndvi_2
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : -0.959918, 0.883263  (min, max)
plot(ndvi_2)

2.2 Método 2: cálculo a través de funciones.

2.2.1 Función calc()`

Este es un procedimiento apropiado para objetos raster de gran tamaño (por ejemplo una escena completa), ya que realiza los cálculos por tramos, ahorrando memoria. Esto significa que para que el resultado de la función definida sea correcto, no debe depender del acceso simultáneo a todos los valores de una vez.

calcNDVI_1 <- function(x) return((x[[5]] - x[[4]]) / (x[[5]] + x[[4]]))

ndvi_3 <- calc(imagen,
                  fun = calcNDVI_1)
plot(ndvi_3)

2.2.2 Función Overlay()

La función Overlay() combinar dos (o más) objetos Raster*, siendo eficiente cuando se aplica sobre grandes conjuntos de datos ráster que no se pueden cargar en la memoria de una vez (de manera similar a calc).

calcNDVI_2 <- function(x, y) return((x - y) / (x + y))

ndvi_4 <- overlay(x = imagen[[5]], 
                 y = imagen[[4]], 
                 fun = calcNDVI_2)
plot(ndvi_4)

Otra manera de calcular la función es la siguiente

vi <- function(x, y) {
  (x - y) / (x + y)
  }

ndvi_5 <- overlay(imagen[[5]], imagen[[4]], 
                  fun = vi)
plot(ndvi_5)

En general, no debe utilizarse el primer método con ficheros de gran tamaño, recomendándose métodos que requieren menos memoria, como el 2º o el 3º. Además, como regla general, si un cálculo necesita usar varias capas individuales por separado (a veces en diferentes objetos), será más fácil configurarlo con la función Overlay() que con la función calc.

rm(ndvi_2, ndvi_3, ndvi_4, ndvi_5, calcNDVI_1, calcNDVI_2, vi)

2.2.3 Opción alternativa

También se puede definir una función directamente para la imagen ráster en R sin usar el comando overlay. En el siguiente ejemplo de código R, se define una otra función para el cálculo del Índice de Vegetación Mejorado (EVI) que tiene la siguiente fórmula:

Fórmula de cálculo del índice NDVI

Esta función incorpora el nombre del objeto ráster para que se pueda aplicar directamente el procedimiento matemático correspondiente.

evi <- function(imagen, i, j, k) {
  bi =imagen[[i]]
  bj =imagen[[j]]
  bk =imagen[[k]]
  evi =2.5 *((bi-bj)/(bi+(6*bj)-(7.5*bk)+1))
  return(evi)
}

Una vez creada la función, se aplica proporcionando a R el nombre y el número de la banda

evi <- evi(imagen = imagen, 5,4,2)
evi                                                            # Propiedades de la capa
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : -3.791283, 1.24265  (min, max)

Como en el caso anterior, se puede representar gráficamente este índice espectral con ggR.

ggR(evi, geom_raster = TRUE) + 
  ggtitle("Índice EVI") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10),         # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),               # size of the legend
        legend.title = element_text(size = 20),      #size of Legend tittle
        legend.text = element_text (size = 15)) +    # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores EVI", colours = c("white", "darkgreen"))

También se puede dibujar el histograma para representar la distribución de valores de EVI en la imagen raster, mostrando que la mayoría de los píxeles tienen valores que van desde -0.5 hasta 1.

histogram(evi, 
          col = "skyblue", 
          main = "Histograma de EVI", 
          xlab = "values de EVI", ylab = "Frecuencia")

2.3 MÉTODO 3: cálculo a través del paquete RSToolbox.

En el paquete RStoolbox contiene una función conocida como spectrallndices() que calcula directamente 23 índices espectrales. El resultado de esta función es el valor de ese índice en forma de una capa ráster (si sólo se calcula un único índice) o un ráster brick (si se calculan varios índices).

La lista de índices calculados por la función spectralIndice con sus respectivas fórmulas, fuente y bandas espectrales aparecen en la siguiente tabla.

ÍNDICE DESCRIPCIÓN FUENTE BANDAS Formula
CLG Green-band Chlorophyll Index Gitelson, 2003 redEdge3, green redEdge3/green - 1
CLRE Red-edge-band Chlorophyll Index Gitelson, 2003 redEdge3, redEdge1 redEdge3/redEdge1 - 1
CTVI Corrected Transformed Vegetation Index Perry, 1984 red, nir (NDVI + 0.5)/sqrt(abs(NDVI + 0.5))
DVI Difference Vegetation Index Richardson, 1977 red, nir s * nir - red
EVI Enhanced Vegetation Index Huete, 1999 red, nir, blue G * ((nir - red)/(nir + C1 * red - C2 * blue + L_evi))
EVI2 Two-band Enhanced Vegetation Index Jiang , 2008 red, nir G * (nir - red)/(nir + 2.4 * red + 1)
GEMI Global Environmental Monitoring Index Pinty, 1992 red, nir (((nir^2 - red^2) * 2 + (nir * 1.5) + (red * 0.5))/(nir + red + 0.5)) * (1 - ((((nir^2 - red^2) * 2 + (nir * 1.5) + (red * 0.5))/(nir + red + 0.5)) * 0.25)) - ((red - 0.125)/(1 - red))
GNDVI Green Normalised Difference Vegetation Index Gitelson, 1998 green, nir (nir - green)/(nir + green)
KNDVI Kernel Normalised Difference Vegetation Index Camps-Valls, 2021 red, nir tanh(((nir - red)/(nir + red)))^2
MCARI Modified Chlorophyll Absorption Ratio Index Daughtery, 2000 green, red, redEdge1 ((redEdge1 - red) - (redEdge1 - green)) * (redEdge1/red)
MNDWI Modified Normalised Difference Water Index Xu, 2006 green, swir2 (green - swir2)/(green + swir2)
MSAVI Modified Soil Adjusted Vegetation Index Qi, 1994 red, nir nir + 0.5 - (0.5 * sqrt((2 * nir + 1)^2 - 8 * (nir - (2 * red))))
MSAVI2 Modified Soil Adjusted Vegetation Index 2 Qi, 1994 red, nir (2 * (nir + 1) - sqrt((2 * nir + 1)^2 - 8 * (nir - red)))/2
MTCI MERIS Terrestrial Chlorophyll Index DashAndCurran, 2004 red, redEdge1, redEdge2 (redEdge2 - redEdge1)/(redEdge1 - red)
1NBRI Normalised Burn Ratio Index Garcia, 1991 nir, swir3 (nir - swir3)/(nir + swir3)
NDREI1 Normalised Difference Red Edge Index 1 GitelsonAndMerzlyak, 1994 redEdge2, redEdge1 (redEdge2 - redEdge1)/(redEdge2 + redEdge1)
NDREI2 Normalised Difference Red Edge Index 2 Barnes, 2000 redEdge3, redEdge1 (redEdge3 - redEdge1)/(redEdge3 + redEdge1)
NDVI Normalised Difference Vegetation Index Rouse, 1974 red, nir (nir - red)/(nir + red)
NDVIC Corrected Normalised Difference Vegetation Index Nemani, 1993 red, nir, swir2 (nir - red)/(nir + red) * (1 - ((swir2 - swir2ccc)/(swir2coc - swir2ccc)))
NDWI Normalised Difference Water Index McFeeters, 1996 green, nir (green - nir)/(green + nir)
NDWI2 Normalised Difference Water Index Gao, 1996 nir, swir2 (nir - swir2)/(nir + swir2)
NRVI Normalised Ratio Vegetation Index Baret, 1991 red, nir (red/nir - 1)/(red/nir + 1)
REIP Red Edge Inflection Point GuyotAndBarnet, 1988 red, redEdge1, redEdge2, redEdge3 0.705 + 0.35 * ((red + redEdge3)/(2 - redEdge1))/(redEdge2 - redEdge1)
RVI Ratio Vegetation Index red, nir red/nir
SATVI Soil Adjusted Total Vegetation Index Marsett, 2006 red, swir2, swir3 (swir2 - red)/(swir2 + red + L) * (1 + L) - (swir3/2)
SAVI Soil Adjusted Vegetation Index Huete, 1988 red, nir (nir - red) * (1 + L)/(nir + red + L)
SLAVI Specific Leaf Area Vegetation Index Lymburger, 2000 red, nir, swir2 nir/(red + swir2)
SR Simple Ratio Vegetation Index Birth, 1968 red, nir nir/red
TTVI Thiam’s Transformed Vegetation Index Thiam, 1997 red, nir sqrt(abs((nir - red)/(nir + red) + 0.5))
TVI Transformed Vegetation Index Deering, 1975 red, nir sqrt((nir - red)/(nir + red) + 0.5)
WDVI Weighted Difference Vegetation Index Richardson, 1977 red, nir nir - s * red

Algunos índices enumerados anteriormente requieren coeficientes adicionales que no serán descritos aquí.

Si desea calcular cualquiera de los índices enumerados anteriormente, debe incluirse el nombre del índice en el argumento index = de la función Rstoolbox::spectralindices(). Los nombres de las bandas en la imagen ráster de entrada también deben incluirse en los argumentos de la función. Por ejemplo, para el cálculo del índice “Specific Leaf Area Vegetation Index (SLAVI)”.

  • Argumento índice = "SLIWI".
  • Argumentos red = "red", nir = "nir"y swir2 = "swir1" ya que se requieren las bandas del rojo, nir y swir2. La sustitución del argumento swir1 por swir2 se debe a que el rango de longitud de onda de SWIR1 en el sensor Landsat 8 OLI & TIRS es de 1570 a 1650 nm, lo que corresponde al rango de longitud de onda swir2.
slavi <- spectralIndices(imagen, 
                         red = "red", nir = "NIR", swir2 = "SWIR1", 
                         index = "SLAVI")
slavi               
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : SLAVI 
## values     : 0.02063482, 3.717375  (min, max)

La distribución de valores del índice es la siguiente.

histogram(slavi, 
          col = "wheat", 
          main = "HistogramA del índice SLAVI" , 
          xlab ="valores", 
          ylab = "Frecuencia")

En el histograma podemos ver que la mayoría de los píxeles en SLAVI están en el rango the de 0.5 - 0.8. Además, podermos ver la distribución espacial del índice representándolo gráficamente.

ggR(slavi, stretch ="lin", geom_raster = TRUE) + 
  ggtitle("Índice SLAVI") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10),         # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),               # size of the legend
        legend.title = element_text(size = 20),      #size of Legend tittle
        legend.text = element_text (size = 15)) +    # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores SLAVI",colours = c("white","darkgreen"))

Es posible calcular todos los índices espectrales disponibles en la función spectralIndices() al mismo tiempo, proporcionando todos los argumetnos y los coeficientes requeridos. Como no tenemos los valores de los coeficientes para ‘swir2ccc’ y ‘swir2coc’, por lo tanto, todos los índices espectrales excepto el Índice de Vegetación de Diferencia Normalizada Corregido (NDVIC) se calcularán utilizando el código R ~ dado.

SI <- spectralIndices(imagen,
                      blue = "blue", green = "green", red = "red", nir = "NIR", swir2 = "SWIR1", swir3 ="SWIR2", 
                      Coefs = list(L = 0.5, G = 2.5, L_evi = 1, C1 = 6, C2 = 7.5, s = 1))
SI    
## class      : RasterBrick 
## dimensions : 1245, 1497, 1863765, 23  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2023-05-07_230110.46043_3800_74316.grd 
## names      :          CTVI,           DVI,           EVI,          EVI2,          GEMI,         GNDVI,         KNDVI,         MNDWI,         MSAVI,        MSAVI2,          NBRI,          NDVI,          NDWI,         NDWI2,          NRVI, ... 
## min values : -0.6781725276, -0.1207067296, -0.6921978226, -0.9936967801, -2.2579093183, -0.9733607897,  0.0000000000, -0.6122067856, -0.4613329627, -0.1833715458, -0.6914897308, -0.9599179772, -0.7963592630, -0.7273798924, -0.8832629723, ... 
## max values :     1.1761220,     0.6577497,     0.9999864,     1.0000000,     1.0544429,     0.7963593,     0.7265791,     0.9960576,     0.7377701,     0.8439642,     0.9918455,     0.8832630,     0.9733608,     0.9866543,     0.9599180, ...

La salida de esta función es un RasterBrick en el que cada capa contiene valores de un índice específico. Se podrían mostrar los nombres de todas las bandas (índices espectrales) calculados mediante este procedimiento.

names(SI)                                 # Lista de índices calculados
##  [1] "CTVI"   "DVI"    "EVI"    "EVI2"   "GEMI"   "GNDVI"  "KNDVI"  "MNDWI" 
##  [9] "MSAVI"  "MSAVI2" "NBRI"   "NDVI"   "NDWI"   "NDWI2"  "NRVI"   "RVI"   
## [17] "SATVI"  "SAVI"   "SLAVI"  "SR"     "TTVI"   "TVI"    "WDVI"

En el procedimiento anterior se calcularon todos los índices espectrales disponibles en RStoolbox:spectralIndices excepto el Índice de Vegetación de Diferencia Normalizada Corregida (NDVIC), debido a que no se han proporcionado los coeficientes requeridos para su cálculo. Podría consultarse los atributos de un índice específico (capa ráster) seleccionándolo con el signo $.

SI$NDWI             # Propiedades de la capa ráster correspondiente
## class      : RasterLayer 
## band       : 13  (of  23  bands)
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2023-05-07_230110.46043_3800_74316.grd 
## names      : NDWI 
## values     : -0.7963593, 0.9733608  (min, max)

También podrían representarse gráficamente todas las bandas mediante el uso de la función rasterVis:levelplot del paquete rasterVis.

levelplot(SI, 
          main ="Mapa de los índices espectrales", 
          xlab = "Longitud (m)", ylab = "Latitud (m)")

Al mismo tiempo, se puede dibujar un único índice espectral de ese RasterBrick seleccionando su nombre acompañado del signo $.

ggR(SI$TVI, geom_raster = TRUE) + 
  ggtitle("Índice TVI") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10),         # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),               # Tamaño de la leyenda
        legend.title = element_text(size = 20),      # Tamaño del título de la leyenda
        legend.text = element_text (size = 15)) +    # Tamaño del texto de la leyenda
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores Índice TVI", colours = c("skyblue", "wheat", "darkgreen"))

3 ÁLGEBRA DE MAPAS SOBRE IMÁGENES RASTER

Bajo el término álgebra de mapas se agrupan un conjunto de técnicas y procedimientos que, operando sobre una o varias capas de una imagen en formato ráster, permiten obtener información derivada, generalmente en forma de nuevas capas. Al igual que sobre un fichero ráster clásico, se pueden realizar diferentes operaciones matemáticas sobre una imagen de satélite usando el paquete raster.

Un ejemplo sencillo es el siguiente: multiplicación de los valores del objeto NDVI por 1000, convirtiendo posteriormente esos valores numericos float a integer.

ndvi_integer <- as.integer(ndvi*10000)
ndvi_integer
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : -9599, 8832  (min, max)
plot(ndvi_integer)

Las operaciones matemáticas se pueden realizar simultáneamente sobre bandas diferentes. Por ejemplo, se pueden realizar varias operaciones sobre la bandas correspondientes a los índices NDVI y EVI para obtener un índice hipotético combinado.

algebra1 <- ((ndvi + evi)*1000) / 2
algebra1
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : -2133.383, 1062.632  (min, max)

La salida, también en formato ráster, se representará gráficamente para comprobar los resultados..

ggR(algebra1, geom_raster = TRUE) + 
  ggtitle("Índice adición NVDI + EVI") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10),         # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),               # size of the legend
        legend.title = element_text(size = 20),      #size of Legend tittle
        legend.text = element_text (size = 15)) +    # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores NDVI + EVI", colours = c("white", "brown", "green"))

De manera similar, es posible una una operación matemática que involucre varias bandas de una única imagen.

algebra2 <- (((imagen[[7]]-imagen[[5]]) / (imagen[[7]]+imagen[[5]])) *100)
algebra2                                                                                    
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : -341.3234, 69.14897  (min, max)

La representación gráfica de ese nuevo índice será la siguiente.

ggR(algebra2, geom_raster = TRUE) + 
  ggtitle("Índice adición varias bandas") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10),         # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),               # size of the legend
        legend.title = element_text(size = 20),      #size of Legend tittle
        legend.text = element_text (size = 15)) +    # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores complejos", colours = c("white", "darkgreen"))

Si desea realizar alguna operación matemática que involucre todos los píxeles de una capa ráster, es conveniente definir una función. En el siguiente ejemplo, se aplica álgebra de mapas para cambiar la escala de los valores NDVI de 0 a 100.

rescale <- function(x) {
  y = (x-min(x))/(max(x)-min(x))*100
       return(y)
}

A continuación, se aplicará esta función a la imagen raster

algebra3 <- calc(ndvi, fun = rescale)
algebra3
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 100  (min, max)

Puede ver en el gráfico ráster que los valores de NDVI se han reescalado a ese nuevo rango. A aquellas superficies con un valor de NDVI más bajo (como cuerpos de agua) se les ha asignado un valor de 0 y a aquellas con un valor más elevado (vegetación sana y densa) se les asigna el valor 100 valor.

ggR(algebra3, geom_raster = TRUE) + 
  ggtitle("Índice adición varias bandas") +                   # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +                    # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,                 # Título centrado
                                  size =20),                  # Tamaño del título
        axis.title = element_text(size =10),                  # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),                        # Tamaño de la leyenda
        legend.title = element_text(size = 20),               # Tamaño del título de la leyenda
        legend.text = element_text (size = 15)) +             # Tamaño del texto de la leyenda

  # Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores complejos", colours = c("white", "darkgreen"))

4 RECLASIFICACIÓN MEDIANTE UMBRALIZACIÓN Y DENSITIY SLICING

El establecimiento de umbrales sobre los valores de un índice espectral es un proceso por el cual se retienen los píxeles cuyos valores se encuentran dentro de unos límites especificados por el usuario en la imagen entrada.

Umbralizar un objeto ráster es un tipo especial de cuantificación que segmenta los valores de los píxeles en dos clases, dependiendo de un valor. La finalidad es identificar qué píxeles conforman los fenómenos espaciales buscados, asignándolos a cierto grupo, llamado comúnmente “segmento”, y qué píxeles son sólo el entorno de estos objetos.

El objeto a segmentar está compuesta por valores numéricos, por lo que la pertenencia de un píxel a un cierto segmento se decide mediante la comparación de su valor con el valor umbral.

La clave de este método es la elección de dicho valor umbral. Este valor se puede elegir de manera manual, o apoyándose en ciertas herramientas estadísticas.

4.1 Métodos manuales

Un enfoque tradicional que se emplea en diferentes estudios es el valor cero, esto significa que los pixeles que tengan un valor de un índice mayor que cero se clasificarán como perteneciente a una determinada clase, por ejemplo, vegetación o agua. Sin embargo, este método puede introducir errores debido a cambios espacio temporales en el brillo y contraste de las imágenes de teledetección. Si se quisiera identificar las zonas con vegetación sana en un ráster de salida, se podría usar ese umbral con la función raster::reclassify`, que identificaría como tales todos los píxels con valores superiores a 0, mientras que se asignaría NA a los píxeles con valores inferiores a 0.

m1 <- c(-Inf, 0, NA,  0, Inf, 2)
matriz_clas_veg1 <- matrix(m1, ncol=3, byrow=TRUE)
matriz_clas_veg1
##      [,1] [,2] [,3]
## [1,] -Inf    0   NA
## [2,]    0  Inf    2
veg1 <- reclassify(ndvi, matriz_clas_veg1)
veg1
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 2, 2  (min, max)

A continuación, se puede representar gráficamente la imagen resultado de esta reclasificación.

ggR(veg1, geom_raster = TRUE) + 
  ggtitle("Vegetación") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10),         # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),               # size of the legend
        legend.title = element_text(size = 20),      #size of Legend tittle
        legend.text = element_text (size = 15)) +    # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores NDVI", colours = c("lightgreen", "darkgreen"))

Para determinar el valor de ese umbral se puede explorar también la distribución de valores incluidos en el objeto ráster analizado mediante la función hist(), que produce un histograma a partir del cual se pueden identificar tanto valores erróneos, valores extremos pero también cambios en la distribución de frecuencias que pudieran sugerir discontinuidades (umbrales).

hist(ndvi,
     main = "Distribución de valores del NDVI según píxeles",
     xlab = "NDVI",
     ylab= "Frecuencia",
     col = "wheat",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

De acuerdo con el anterior ejemplo, si consideramos que los valores teóricos posibles de NDVI oscilan entre -1 y + 1, los píxeles con valores superiores a 0,25 probablemente correspondan a vegetación. Si se quisiera identificar las zonas con vegetación sana en un ráster de salida, se podría usar ese umbral con la función raster::reclassify`, que identificaría como tales todos los píxels con valores superiores a 0.25, mientras que se asignaría NA a los píxeles con valores inferiores a 0.25.

m2 <- c(-Inf, 0.25, NA,  0.25, Inf, 2)
matriz_clas_veg2 <- matrix(m2, ncol=3, byrow=TRUE)
matriz_clas_veg2
##      [,1] [,2] [,3]
## [1,] -Inf 0.25   NA
## [2,] 0.25  Inf    2
veg2 <- reclassify(ndvi, matriz_clas_veg2)
veg2
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 2, 2  (min, max)

A continuación, se puede representar gráficamente la imagen resultado de esta reclasificación.

ggR(veg2, geom_raster = TRUE) + 
  ggtitle("Vegetación") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10),         # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),               # size of the legend
        legend.title = element_text(size = 20),      #size of Legend tittle
        legend.text = element_text (size = 15)) +    # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores NDVI", colours = c("lightgreen", "darkgreen"))

Si desea retener píxeles con valores NDVI que van desde 0,1 a 0,3, probablemente representando suelos sin vegetación, se puede asignar el valor ‘1’ a todos los píxeles dentro del rango de valores que va de 0.1 a 0.3, y NA a los píxeles por encima y por debajo de esos límites

m3 <- c(-Inf, 0.1, NA,0.1, 0.3, 1, 0.3, Inf, NA)
matriz_clas_suelo <- matrix(m3, ncol=3, byrow=TRUE)
matriz_clas_suelo
##      [,1] [,2] [,3]
## [1,] -Inf  0.1   NA
## [2,]  0.1  0.3    1
## [3,]  0.3  Inf   NA
suelo <- reclassify(ndvi, matriz_clas_suelo)

De nuevo, se representa gráficamente

ggR(suelo, geom_raster = TRUE) + 
  ggtitle("Suelo desnudo") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10),         # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),               # size of the legend
        legend.title = element_text(size = 20),      #size of Legend tittle
        legend.text = element_text (size = 15)) +    # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores", colours = c("white", "brown", "white"))

Un proceso algo más complejo es la definición de sucesivos umbrales, que conducen a la aparición de más clases, técnica que en el ámbito del análisis de imágenes se conoce habitualmente con el nombre de Density slicing. Este procedimiento divide los valores de los píxeles en diferentes rangos y asigna un valor o color específico a cada rango en el ráster de salida.

Este procedimiento se realiza también mediante la función raster::reclassify. En esta función, los argumentos principales son la imagen ráster de entrada y la matriz que sirve de guía para la reclasificación. Esta matriz tiene tres columnas.

  • La primera columna y la segunda columna contienen los límites más bajo y más alto de un rango particular.
  • La tercera columna contiene el valor que se asignará a ese rango.

Un ejemplo de esta matriz de reclasificación es la siguiente:

ds_matrix <- matrix(c(-Inf, 0.1, 1, 
                      0.1, 0.3, 2, 
                      0.3, 0.5, 3, 
                      0.5, Inf, 4),
                      ncol = 3,                                                 # Matriz con 3 columnas
                      byrow = TRUE)                                             # Lectura en horizontal
ds_matrix
##      [,1] [,2] [,3]
## [1,] -Inf  0.1    1
## [2,]  0.1  0.3    2
## [3,]  0.3  0.5    3
## [4,]  0.5  Inf    4
ndvi_2 <- round(ndvi, 
                digits = 2)                                        # Redondea a 2 dígitos el raster con valores ndvi

ds_ndvi <- reclassify(ndvi_2, 
                      ds_matrix)
ds_ndvi
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 4  (min, max)
ggR(ds_ndvi, geom_raster = TRUE) + 
  ggtitle("Reclasificación en 4 categorías") +        # Título de la imagen.
  labs(x="Longitud(m)", y="Latitud (m)") +            # Etiquetas del eje.
  theme(plot.title = element_text(hjust =0.5,         # Título centrado.
                                  size =20),          # Tamaño del título.
        axis.title = element_text(size =10),          # Tamaño de las etiquetas de los ejes.
        legend.key.size =unit(1,"cm"),                # Tamaño de la leyenda.
        legend.title = element_text(size = 20),       # Tamaño del título de la leyenda.
        legend.text = element_text (size = 15)) +     # Tamaño del texto de la leyenda.
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores", 
                     colours = c("turquoise2", "burlywood", "green3", "darkgreen"))

###Análisis estadístico de la distribución de categorías

Una vez definidas umbrales correspondientes a las categorías, podría ser interesante extraer algún tipo de información estadística sobre la frecuencia con la que aparecen en la escena analizada. Para realizarlo, es necesario transformar las categorías en un factores, para posterioremente representarlas gráficamente y calcular algún valor adicional.

El primer paso consiste en transformar los valores correspondientes a cada categoría en un factor, con la función raster::ratify.

ds_ndvi_factor <- ratify(ds_ndvi)                            # Define los valores del ráster como factores (variable categórica)
ds_ndvi_factor
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 4  (min, max)
## attributes :
##  ID
##   1
##   2
##   3
##   4

A continuación, se crea la rat o tabla de atributos del ráster utilizando para ello la función `base::levels

rat <- levels(ds_ndvi_factor)[[1]]
rat
##   ID
## 1  1
## 2  2
## 3  3
## 4  4

A continuación, se añade a cada categoría el siguiente texto, y dicho texto se extiende al raster con las categorías numéricas.

rat$legend  <- c("Clase 1", "Clase 2", "Clase 3", "Clase 4") 
levels(ds_ndvi_factor) <- rat

Con la siguiente sintaxis se crean los colores que van a servir para identificar cada categoría en el diagrama de barras

my_col=c("darkolivegreen","yellow2", "orange2", "red")

Finalmente, se puede construir con toda esta información un diagrama de barras

options(scipen=999)                                                             #Sustituye la notación científica por 0.
barplot(ds_ndvi_factor,
        main = "Distribución de categorías de superficie",
        col = my_col,
        names.arg = c("C1", "C2", "C3", "C4"), 
        las=0)
## Warning in .local(height, ...): a sample of 53.6% of the raster cells were used
## to estimate frequencies

El gráfico anterior se puede mejorar, añadiendo una etiqueta con el número de casos por categoría.

vals <- getValues(ds_ndvi_factor)

pix_1 <- length(subset(vals, vals == 1))
pix_2 <- length(subset(vals, vals == 2))
pix_3 <- length(subset(vals, vals == 3))
pix_4 <- length(subset(vals, vals == 4))

pix_values <- c(pix_1,pix_2,pix_3,pix_4)

bp <- barplot(ds_ndvi_factor,
        main = "Número de píxeles según categoría",
        col = my_col,
        names.arg = c("C1", "C2", "C3", "C4"), 
        las=0)
## Warning in .local(height, ...): a sample of 53.6% of the raster cells were used
## to estimate frequencies
text(bp, y=70000, labels=pix_values) #y chosen manually

Otra variante es incluir el porcentaje respecto al número de píxeles de la imagen.

pix_total <- length(vals)
area_1_perc <- round(pix_1/pix_total*100, digits = 5)
area_2_perc <- round(pix_2/pix_total*100, digits = 5)
area_3_perc <- round(pix_3/pix_total*100, digits = 5)
area_4_perc <- round(pix_4/pix_total*100, digits = 5)

area_perc_covered <- c(area_1_perc,area_2_perc,area_3_perc,area_4_perc)

bp <- barplot(ds_ndvi_factor,
        main = "Porcentaje de píxeles ocupados por cada categoría",
        col = my_col,
        names.arg = c("C1", "C2", "C3", "C4"), 
        las=0)
## Warning in .local(height, ...): a sample of 53.6% of the raster cells were used
## to estimate frequencies
text(bp, y=70000, labels=area_perc_covered)

Finalmente, es posible calcular ese porcentaje, pero respecto a la extensión del ráster en kilómetros cuadrados

ex <- extent(ds_ndvi_factor)
x <- (ex[2]-ex[1])/1000
y <- (ex[4]-ex[3])/1000
area_total <- x*y

area_1 <- round((area_1_perc/100)*area_total, digits = 5)
area_2 <- round((area_2_perc/100)*area_total, digits = 5)
area_3 <- round((area_3_perc/100)*area_total, digits = 5)
area_4 <- round((area_4_perc/100)*area_total, digits = 5)

area_covered <- c(area_1,area_2,area_3,area_4)

bp <- barplot(ds_ndvi_factor,
        main = "Área de cada categoría",
        col = my_col,
        names.arg = c("C1", "C2", "C3", "C4"), 
        las=0,
        ylim = c(0, 800000))
## Warning in .local(height, ...): a sample of 53.6% of the raster cells were used
## to estimate frequencies
text(bp, y=75000, labels=area_covered) #y chosen manually

Finalmente,

tabla <- as.data.frame(rbind(pix_values, area_perc_covered, area_covered))
names(tabla) <- c("Cl1", "Cl2", "Cl3", "Cl4")

4.2 Procedimiento guiado

El método de Otsu, llamado así en honor a Nobuyuki Otsu en 1979 (Otsu, 1979), utiliza técnicas estadísticas, para resolver el problema. En concreto, utiliza la varianza, una medida de la dispersión de valores. Su método calcula el valor umbral de forma que la dispersión dentro de cada segmento sea lo más pequeña posible, pero al mismo tiempo la dispersión sea lo más alta posible entre segmentos diferentes. Para ello se calcula el cociente entre ambas varianzas y se busca un valor umbral para el que este cociente sea máximo.

A continuación se mostrará una aplicación práctica de esta metodología, aplicándose a un raster en el que se ha calulado el índice MNDWI, cuyos valores oscilan entre -1 a 1: los valores negativos corresponden a terreno continental y los positivos a agua.

\[NDWI = (G - NIR)/(G + SWIR)\]

El primer procedimiento consiste en calcular el citado índice

imagen
## class      : RasterBrick 
## dimensions : 1245, 1497, 1863765, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : coastal.aerosol,          blue,         green,           red,           NIR,         SWIR1,         SWIR2 
## min values :    0.0964179114,  0.0748399049,  0.0425921641,  0.0208406653,  0.0008457669, -0.0078721829, -0.0050529451 
## max values :       0.7346282,     0.7177562,     0.6924697,     0.7861769,     1.0124315,     1.0432045,     1.1179360
ndwi <- (imagen$green - imagen$NIR) / (imagen$green + imagen$NIR)

A continuación se deberá cargar y activar el paquete colorspace para obtener diferentes paletas de colores.

library(colorspace)
## 
## Attaching package: 'colorspace'
## The following object is masked from 'package:raster':
## 
##     RGB

La paleta de colores elegida es divergente.

paleta <- diverging_hcl(n= 11,
                        palette= "Broc", 
                        rev=T)

A continuación se utiliza esa paleta para colorear el mapa:

plot(ndwi, 
    col=colorRampPalette(paleta)(16), 
    axes=F,
    box=F,
    legend=F,
    zlim=c(-1,1))

Para umbralizar, destacando y extrayendo los cuerpos de agua de mayor tamaño utilizando el método de Otsu se recurre a una función del paquete EBImage. EBImage proporciona una funcionalidad de propósito general para el procesamiento y análisis de imágenes (Pau et al. [2010]) y es un paquete creado y mantenido por la iniciativa Bioconductor. Para la mantención e instalación de sus propios paquetes instalamos en primer lugar el paquete para gestión, BiocManager (Morgan [2021]).

#install.packages("BiocManager")

Y para activar en la sesión actual:

library(BiocManager)

A continuación, se puede instalar el paquete deseado (por ejemplo EBImage) utilizando la función install() del paquete BiocManager.

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("EBImage", force = TRUE)
## Bioconductor version 3.17 (BiocManager 1.30.20), R 4.3.0 (2023-04-21 ucrt)
## Installing package(s) 'EBImage'
## package 'EBImage' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\Domingo\AppData\Local\Temp\RtmpY3B5JG\downloaded_packages
## Installation paths not writeable, unable to update packages
##   path: C:/Program Files/R/R-4.3.0/library
##   packages:
##     class, KernSmooth, MASS, nnet
## Old packages: 'cachem', 'checkmate', 'evaluate', 'fs', 'gdata', 'later',
##   'profvis', 'rlang', 'sass', 'servr', 'testthat', 'tseries', 'units',
##   'viridis', 'viridisLite', 'vroom', 'waldo', 'wk', 'xml2'

Y luego cargar en la sesión actual:

library(EBImage)
## 
## Attaching package: 'EBImage'
## The following objects are masked from 'package:raster':
## 
##     flip, rotate

Las funciones del paquete requieren la información de nuestro modelo ráster en formato de array o matriz, para ello utilizamos la función correspondiente as.array(x) que toma un raster x y devuelve como array

arr <- as.array(ndwi)
str(arr)
##  num [1:1245, 1:1497, 1] -0.465 -0.469 -0.497 -0.545 -0.564 ...

La función EBIimage::otsu() devuelve un valor umbral o de corte para binarizar la imagen

corte <- otsu(x = arr, 
              range = c(-1,1), 
              levels=256)                              # Number of grayscale levels
corte
## [1] -0.08984375

Aplicamos el valor de corte, pero antes duplicaremos la imagen original para trabajar sobre ella.

agua <- ndwi                                           
agua[agua >= corte] <- 1
agua[agua < corte] <- NA
plot(agua, 
     col= "blue", 
     axes=F,
     box=T, 
     legend=F)

4.3 Procesamiento de Ráster Binario

El resultado del proceso anterior es una imagen binaria, valor 1 para los píxeles asignados como cuerpo de agua y el resto NA.

agua
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 1  (min, max)

A continuación vamos a revisar algunos post procesos, operaciones o filtros sobre la imagen de este tipo.

4.3.1 Parches

Cada grupo de píxeles asignados con valor uno adyacentes o en contacto forman un “parche”.

El proceso de identificar cada uno de ellos lo realizaremos con la función raster::clump del paquete raster. Esta función detecta grupos (parches) de píxeles conectados. Cada grupo obtiene una identificación única. NA y cero se utilizan como valores de fondo (es decir, estos valores se utilizan para separar grupos). Para la definición de las celdas adjacentes se puede usar el kernel de vecinos tipo reina o torre , usando el argumento directions. Cuando todos los vecinos correspondientes al tipo pertenecen al parche, el píxel central se asigna a ese id.

Si el argumento gaps=TRUE los ID no son consecutivos (por ejemplo, puede haber parches con los IDs 1, 2, 3 and 5, pero no el 4). Si se establece como FALSE, los números serán consecutivos.

También es necesario instalar el paquete igraph para usar esta función.

# install.packages("igraph")
library(igraph)
## 
## Attaching package: 'igraph'
## The following object is masked from 'package:EBImage':
## 
##     normalize
## The following object is masked from 'package:raster':
## 
##     union
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
agua_id <- clump(x = agua, 
                directions = 8,
                gaps = FALSE)

El argumento range() de la función cellStats identifica el máximo ID devuelto por la función anterior (que coincide con el número de ‘parches’ encontrados)

cellStats(x = agua_id, 
          stat = range)
## [1]    1 1818

Puede observarse que se han detectado 1818 “parches” con píxeles identificados como agua. Su numeración corresponde a un orden aleatorio de definición y no tiene relación con la ubicación geográfica, ni con su tamaño.

paleta <- qualitative_hcl(n= 11, 
                          palette= "Cold", 
                          rev=F)
plot(agua, 
    col=colorRampPalette(paleta)(64), 
    axes = F,
    box = T, 
    legend= F)

4.3.2 Cálculo de Áreas

El paquete bfastSpatial (Dutrieux and DeVries, 2014) pre procesa datos de series de tiempo en formato ráster para que puedan ser analizados con algoritmos de detección de cambios como bfast1. Utiliza clases del paquete ráster e incluye utilidades para 1 BFAST, (Breaks for Additive Season and Trend), es un algoritmo que integra la descomposición de series de tiempo en componentes de tendencia, temporada y resto con métodos para detectar y caracterizar cambios dentro de series de tiempo, ejecutar dichos algoritmos y post procesar los resultados. El paquete se aloja en el sitio github.com y existe un instalador de paquetes especialmente diseñado para esta modalidad. Wickham et al. [2021] desarrollan el paquete devtools con una colección de herramientas para el desarrollo de paquetes.

# install.packages("devtools")
library(devtools)
## Loading required package: usethis
## 
## Attaching package: 'devtools'
## The following object is masked from 'package:BiocManager':
## 
##     install

A continuación, instalamos el paquete bfastSpatial (y su dependencia gdalUtils) usando la función install_github():

#devtools:::install_github("gearslaboratory/gdalUtils")
#devtools:::install_github('dutri001/bfastSpatial')

Y agregamos a nuestra sesión actual:

library(bfastSpatial)
## Loading required package: parallel
## Loading required package: bfast
## Loading required package: strucchangeRcpp
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
## Attaching package: 'bfast'
## The following object is masked _by_ '.GlobalEnv':
## 
##     ndvi
## Loading required package: gdalUtils
## 
## Attaching package: 'gdalUtils'
## The following object is masked from 'package:sf':
## 
##     gdal_rasterize
## Loading required package: stringr
## 
## Attaching package: 'stringr'
## The following object is masked from 'package:strucchangeRcpp':
## 
##     boundary
## Loading required package: rgdal
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-6, (SVN revision 1201)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.6.2, released 2023/01/02
## Path to GDAL shared files: D:/Documentos/R/win-library/4.3/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.2.0, March 1st, 2023, [PJ_VERSION: 920]
## Path to PROJ shared files: D:/Documentos/R/win-library/4.3/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Loading required package: bitops

Para calcular el tamaño de los grupos de píxeles en un RasterLayer. El área del grupo se asigna a cada grupo y se genera un nuevo ráster con las dimensiones idénticas al original, con los valores que representan el tamaño por grupo. La definición de un grupo de píxeles sigue la lógica utilizada en la función clump(), donde las diagonales se incluyen o se ignoran al definir los grupos (diferencia en vecindad reina o torre). Se puede proporcionar un factor f de conversión opcional basado en el tamaño de píxel.

Para nuestro caso, las imágenes Landsat que hemos utilizado, poseen un pixel de 30 m de lado, o también lo podemos definir como un píxel de 900 m2 (30 m × 30 m) y si consideramos que una hectárea son 10.000 \(m^2\) el factor lo podemos calcular como:

\[f = 900/10000\]

areas_agua <- clumpSize(x = agua, 
                        f = 900/10000)
paleta <- sequential_hcl(n= 11, 
                         palette= "Blues 2", 
                         rev=T) 
plot(areas_agua, 
    col=colorRampPalette(paleta)(64), 
    axes = F, 
    box = T, 
    legend=F)

Otra forma de explorar el resultado es mediante el uso de la función hist() para visualizar la distribución de áreas de los polígonos.

areas_agua
## class      : RasterLayer 
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : clumps 
## values     : 0.09, 16321.95  (min, max)
hist(areas_agua,
     breaks = 10,
     xlim = c(16000, 18000))

De la inspección podemos buscar un valor de área para separar los polígonos que deseamos mantener, en nuestro ejemplo seleccionamos el valor de 20.000 ha o de \(2 × 10^8\) m. La función areaSieve(x, thresh) del paquete bfastSpatial aplica un umbral de área thresh a una capa ráster x. Se proporciona un umbral en metros cuadrados (se generalizará en el futuro) y se eliminan todos los ‘grupos’ de píxeles más pequeños que este umbral (Figura 9.40).

agua_mayores <- areaSieve(x = areas_agua, 
                             thresh = 15000)
id_agua_mayores <- clump(agua_mayores)

plot(id_agua_mayores, 
     col= "blue", 
     axes= F, 
     box = T, 
     legend=F)

4.3.3 Conversión a una capa Vectorial

La función raster::boundaries() detecta bordes (bordes son pixeles que tienen más de una clase en la vecindad denominadas ‘reina’ o ‘torre’), asignando un valor cero a aquellos pixeles dentro del área y un valor uno a los pixeles límite o borde.

agua_borde <- boundaries(x = id_agua_mayores, 
                           type="inner", 
                           classes=F, 
                           directions=8)

En el siguiente paso vamos a signar los píxeles de borde (1) a NA para dejar solo una clase (cero):

agua_borde[agua_borde == 1] <- NA

Y finalmente la función rasterToPolygons() convierte de modelo ráster a modelo vectorial.

vector_agua <- rasterToPolygons(agua_borde, 
                                dissolve = TRUE)

spplot(vector_agua, 
       colorkey=FALSE, 
       col.regions="gray95")

La conversión directa genera polígonos con un exceso de vértices (Figura 9.41) y una posible solución es la contenida en el paquete rmapshaper (Teucher and Russell, 2021) cuya principal ventaja es la disponibilidad del algoritmo de simplificación con reconocimiento topológico en la función ms_simplify().

#install.packages("rmapshaper")
library(rmapshaper)

Esto significa que los límites compartidos entre polígonos adyacentes siempre se mantienen intactos, sin espacios ni superposiciones, incluso con altos niveles de simplificación (Figura 9.42).

vector_cuerpo_agua <- ms_simplify(vector_agua,
                                  explode = T,
                                  keep = 0.1)

class(vector_cuerpo_agua)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
spplot(vector_cuerpo_agua,
       colorkey=FALSE,
       col.regions="gray95")

Finalmente, es posible exportar ese vector en formato shapefile para usar posteriormente como máscara, previa transformación en objeto sf:

vector_cuerpo_agua_sf <- st_as_sf(vector_cuerpo_agua)

st_write(vector_cuerpo_agua_sf, 
         "D:/G174_OCW/Lab_5_transformaciones_matematicas/datos/vector_agua.shp", 
         append = TRUE)
## Updating layer `vector_agua' to data source `D:/G174_OCW/Lab_5_transformaciones_matematicas/datos/vector_agua.shp' using driver `ESRI Shapefile'
## Updating existing layer vector_agua
## Writing 166 features with 1 fields and geometry type Polygon.

O como en el punto 3.12 se define el uso del paquete mapview para crear sesiones interactivas de cartografía (Figura 9.43) a la cual agregamos los polígonos recién creados.

library(mapview)
mapview(vector_cuerpo_agua)
mapview(vector_cuerpo_agua, 
        map.types = c("OpenTopoMap"), 
        legend=F, 
        alpha.regions=0.35, 
        col.regions="deeppink", 
        color="deeppink2",
        lwd=3)

ACTIVIDAD INDIVIDUAL Cada alumno deberá elegir dos índices espectrales que sean adecuados para el tema elegido. A continuación, deberá realizar una búsqueda en internet para obtener información sobre la fórmula de cálculo de esos índices, las bandas implicadas, los valores recomendados de esos índices así como su significado. Una vez realizado, deberá proceder al análisis de su escena, calculando los índices señalados y obteniendo posibles umbrales empíricos aplicando las técnicas de obtención de umbrales y Density Slicing

ACTIVIDAD INDIVIDUAL

5 TRANSFORMACIÓN “TASSELED CAP”

La transformación Tasseled cap es una transformación (en la que se obtienen otros datos no correlacionados entre sí) de datos en formato ráster en forma de un nuevo sistema de coordenadas para representar la línea del suelo y la vegetación.Es también un método de compresión que reduce el volumen de los datos espectrales de 6 bandas, en unas pocas, lo que que permite comprender fenómenos importantes del desarrollo de cultivos. La región triangular del gráfico muestra la vegetación en varias etapas de crecimiento. Por lo tanto, esta transformación es muy útil para el monitoreo y cartografía de la vegetación.

Relación entre la banda del rojo y el infrarrojo próximo

La transformación Tasseled Cap es una combinación lineal de información la información espectral proporcionada por las imágenes originales en un conjunto menor de bandas, normalmente organizadas en torno a tres ejes o componentes principales conocidos como brillo, verdor y humedad. Esta simplificación de los componentes espectrales originales permite tratarlas como variables físicas.

Por lo tanto, este procedimiento es muy beneficioso para detectar y comparar variaciones en la vegetación, el suelo y las características hechas por el hombre en diferentes períodos de tiempo.

Una de las mayores limitaciones de la transformación Tasseled-cap es que solo se puede aplicar a sensores para los que ya se han desarrollado coeficientes de transformación. Actualmente, esto incluye Landsat, MODIS, ASTER, Ikonos, SPOT y, además, Sentinel 2.

En R esta transformación se puede realizar mediante función RStoolbox::tasseledCap(). Debe incluir el nombre del sensor mediante el argumento sat, mientras que las unidades del ráster de entrada son de reflectancia.

Satélite Bandas
Landsat4TM 1,2,3,4,5,7
Landsat5TM 1,2,3,4,5,7
Landsat7ETM 1,2,3,4,5,7
Landsat8OLI 2,3,4,5,6,7
MODIS 1,2,3,4,5,6,7
QuickBird 2,3,4,5
Spot5 2,3,4,5
RapidEye 1,2,3,4,5

En el siguiente ejemplo de código, efectuaremos una transformación Tasseled cap a una imagen raster.

imagen_tc <- tasseledCap(imagen[[2:7]], 
                          sat ="Landsat8OLI" )
imagen_tc
## class      : RasterBrick 
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      :  brightness,   greenness,     wetness 
## min values :  0.04875321, -0.23735039, -0.80978684 
## max values :   1.9836461,   0.4419796,   0.2490443

Entonces, podemos ver a partir de los atributos de la salida que esta función ha calculado los tres componentes en forma de capas raster separadas de un RasterBrick. A su vez, estas imágenes pueden ser dibujadas para observar las variaciones de los valores de brillo, verdor y humedad correspondientes al suelo, la vegetación y la humedad.

ggR(imagen_tc$brightness, geom_raster = TRUE) + 
  ggtitle("Brightness Tasseled cap") +                  # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +              # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,           # Título centrado
                                  size =20),            # Tamaño del título
        axis.title = element_text(size =10),            # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),                  # size of the legend
        legend.title = element_text(size = 20),         # size of Legend tittle
        legend.text = element_text (size = 15)) +       # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores", colours = c("white", "green", "brown"))

ggR(imagen_tc$greenness, geom_raster = TRUE) + 
  ggtitle("Greenness Tasseled cap") +                  # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +              # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,           # Título centrado
                                  size =20),            # Tamaño del título
        axis.title = element_text(size =10),            # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),                  # size of the legend
        legend.title = element_text(size = 20),         # size of Legend tittle
        legend.text = element_text (size = 15)) +       # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores", colours = c("white", "darkgreen"))

ggR(imagen_tc$wetness, geom_raster = TRUE) + 
  ggtitle("Wetness Tasseled cap") +                  # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +              # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,           # Título centrado
                                  size =20),            # Tamaño del título
        axis.title = element_text(size =10),            # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),                  # size of the legend
        legend.title = element_text(size = 20),         # size of Legend tittle
        legend.text = element_text (size = 15)) +       # size of legend text
# Dibuja un gradiente de colores 
scale_fill_gradientn(name = "valores", colours = c("skyblue", "darkblue"))

6 ANÁLISIS EN COMPONENTES PRINCIPALES

El Análisis de Componentes Principales (PCA) otra transformación matemática que convierte cierto número de variables originales (en este caso, bandas) en combinaciones lineales de componentes ortogonales, combinaciones conocidas como Componentes Principales, cuya importancia decrece progresivamente. Este método transforma los valores de brillo de las bandas de entrada en un espacio de atributo multivariante cuyos ejes se rotan para representar el varianza en cada variable de los datos. El primer eje explica la varianza máxima. El segundo eje es ortogonal al primer eje y explica la varianza restante que no fue explicada por el primer eje, y así sucesivamente. Estos ejes se conocen como componentes principales.

Fórmula de cálculo del índice NDVI

Los componentes principales son ortogonales y no están correlacionados entre sí. La dirección del componente principal es el vector propio y su magnitud es el valor propio. El ángulo de rotación del eje es el ángulo utilizado en la transformación. Estos vectores propios, valores propios y matriz de covarianza calculada (matriz de correlación en el caso de bandas normalizadas) del ráster de entrada se utilizan para desarrollar una fórmula lineal que define el cambio y rotación. Esta fórmula se aplica para transformar los valores de todos los píxeles en el ráster de entrada en relación con el nuevo eje para crear un ráster de salida transformado. El análisis de componentes principales se realiza principalmente para eliminar la redundancia en los datos, ya que los primeros componentes principales explican la mayor parte de la variación en los datos ráster de entrada multibanda.

En R el análisis de componentes principales (PCA) en imágenes ráster multibanda se puede realizar utilizando la función Rstoolbox::rasterPCA().

imagen_pca <- rasterPCA(imagen, 
                        nComp = nlayers(imagen))
imagen_pca
## $call
## rasterPCA(img = imagen, nComp = nlayers(imagen))
## 
## $model
## Call:
## princomp(cor = spca, covmat = covMat[[1]])
## 
## Standard deviations:
##       Comp.1       Comp.2       Comp.3       Comp.4       Comp.5       Comp.6 
## 0.1473147197 0.0845584282 0.0384519375 0.0182971395 0.0073537457 0.0031112840 
##       Comp.7 
## 0.0009613602 
## 
##  7  variables and  1863765 observations.
## 
## $map
## class      : RasterBrick 
## dimensions : 1245, 1497, 1863765, 7  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      :         PC1,         PC2,         PC3,         PC4,         PC5,         PC6,         PC7 
## min values : -0.39137375, -0.40236374, -0.33008778, -0.05978816, -0.10655715, -0.06080169, -0.05798591 
## max values :  1.49606648,  0.68809466,  0.82989370,  0.53319468,  0.22249235,  0.13013959,  0.04241485 
## 
## 
## attr(,"class")
## [1] "rasterPCA" "RStoolbox"

Este procedimiento ha creado un RasterBrick de 10 componentes principales. La desviación estándar y la proporción de varianza explicada por cada componente principal se muestra con la función summary y model.

summary(imagen_pca$model)
## Importance of components:
##                           Comp.1     Comp.2     Comp.3     Comp.4      Comp.5
## Standard deviation     0.1473147 0.08455843 0.03845194 0.01829714 0.007353746
## Proportion of Variance 0.7062085 0.23267753 0.04811462 0.01089449 0.001759778
## Cumulative Proportion  0.7062085 0.93888602 0.98700065 0.99789514 0.999654918
##                              Comp.6        Comp.7
## Standard deviation     0.0031112840 0.00096136024
## Proportion of Variance 0.0003150068 0.00003007551
## Cumulative Proportion  0.9999699245 1.00000000000

De acuerdo con los resultados, el primer componente principal explica el 99,6 % de la varianza en los datos ráster. La proporción de varianza disminuye en los componentes sucesivos. Se puede mostrar la proporción de varianza explicada por cada componente principal en forma gráfica usando el siguiente código.

screeplot(imagen_pca$model, 
          main = "Scree plot de los componentes principales" )

El screeplot muestra claramente que el primer componente principal acumula la máxima varianza. Este primer componente se puede representar gráficamente.

ggRGB(imagen_pca$map,
      1,2,3, 
      stretch="lin", q=0)

También se podrían representar los siguientes simultáneamente

levelplot(imagen_pca$map, 
          main ="Componentes principales", 
          xlab = "Longitud (m)", ylab ="Latitud (m)")

Dados los resultados, la imagen original se puede explicar de manera efectiva mediante el primer componente principal.

7 PAN-SHARPENING

La mayoría de las imágenes multiespectrales distribuidas de manera gratuita a menudo solo tienen una resolución espacial media. Algunos sensores (incluidos Landsat 7 y 8) poseen una única banda de mayor resolución espacial (15 m en el caso de Landsat) que se puede utilizar para mejorar las imágenes de resolución media. Estas bandas se llaman bandas pancromáticas (o pans), porque tienen un ancho de banda más amplio, cubriendo la mayor parte del visible y, hasta cierto punto, también el rango de longitud de onda NIR. Se podría pensar en ellos como la suma de los canales rojo, verde y azul con el doble de la resolución espacial.

La idea del enfoque panorámico es fusionar una imagen multiespectral de resolución media con una banda ráster pancromática de resolución más fina, creando un nuevo conjunto de datos ráster multibanda con las propiedades espectrales y la resolución espacial de la banda pancromática. Para ello se han propuesto diferentes métodos:

Otras alternativas populares al enfoque IHS son:

Mientras que el primer enfoque conserva un buen equilibrio de color, el segundo conduce a distorsiones de color, pero mejora la nitidez de los bordes y puede manejar más de tres bandas. Además, en el caso de la transformación de Brovey es fundamental que los valores pancromáticos y multiespectrales tengan la misma escala numérica (por ejemplo, valores de 0 a 1), pues de lo contrario los colores estarán apagados.

Como se puede deducir, estas transformaciones alteran considerablemente los valores de reflectancia reales y, por lo tanto, se utilizan principalmente para fines de visualización. No obstante, si se quisieran usar para una clasificación, normalmente es mejor usar la banda pancromática tal como está (si es que lo hace, ya que en su núcleo es solo la integral de la banda azul, verde y roja) o como la imagen de origen para el cálculo de las métricas de textura.

El enfoque panorámico requiere un registro conjunto muy preciso, de lo contrario sus resultados aparecerán borrosos. Incluso con bandas pancromáticas de la misma plataforma, a menudo es una buena idea ver si el registro conjunto automático puede mejorar el ajuste.

El “pan-sharpening” no se limita, por supuesto, a las bandas pancromáticas del mismo satélite, sino que también puede ser usada, por ejemplo, para mejorar la imagen Landsat 5 a partir de una imagen de un satélite con mayor resolución, como IKONOS.

Este procedimiento se puede realizar mediante la función RStoolbox::pansharpen. La imagen de entrada aparece en el argumen img y la imagen pancromática en pan. Aunque hay varios métodos disponibles, la función sólo incluye 3: el análisis de componentes principales, el método de Brovey y el de Intensidad - Tono - Método de saturación.

En esta sección, realizaremos un enfoque panorámico con cada método disponible en R. En primer lugar, importemos tanto el ráster multiespectral de 30 m de resolución espacial como el ráster pancromático de 15 m de resolución. Ambos conjuntos de datos ráster deben tener valores de píxeles en la misma escala.

Primero, importaremos como un fichero raster brick una imagen con la resolución de 30 m.

ps_input <- brick("./datos/originales/pan_sharpening_input.tif")
ps_input               
## class      : RasterBrick 
## dimensions : 333, 333, 110889, 4  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 400005, 409995, 2800005, 2809995  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=44 +datum=WGS84 +units=m +no_defs 
## source     : pan_sharpening_input.tif 
## names      : pan_sharpening_input_1, pan_sharpening_input_2, pan_sharpening_input_3, pan_sharpening_input_4 
## min values :                  10463,                   9852,                   9074,                   9300 
## max values :                  17455,                  17995,                  19205,                  23548

A continuación se importará el fichero con la banda pancromática de resolución espacial 15 m.

pan_band <- raster("./datos/originales/pan_band_15m.tif")
pan_band            
## class      : RasterLayer 
## dimensions : 666, 666, 443556  (nrow, ncol, ncell)
## resolution : 15, 15  (x, y)
## extent     : 399997.5, 409987.5, 2799998, 2809988  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=44 +datum=WGS84 +units=m +no_defs 
## source     : Pan_band_15m.tif 
## names      : Pan_band_15m 
## values     : 9376, 19793  (min, max)

La banda pancromática puede ser representada gráficamente.

ggR(pan_band, stretch = "lin") + 
  ggtitle("Imagen Pancromática") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +            # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,         # Título centrado
                                  size =20),          # Tamaño del título
        axis.title = element_text(size =10),          # Tamaño de las etiquetas de los ejes
        legend.key.size =unit(1,"cm"),                # size of the legend
        legend.title = element_text(size = 20),       #size of Legend tittle
        legend.text = element_text (size = 15))       # size of legend text

Los cambios en las resoluciones espaciales antes y después de la modificación se pueden ver con la función raster::res():

res(ps_input)
## [1] 30 30
res(pan_band)
## [1] 15 15

7.1 Método de análisis de componentes principales (PCA)

En este procedimiento, el primer componente principal se cambia por la banda pancromática. La sintaxis de R requiere el argumento method = pca en la función panSharpen.

ps_pca <- panSharpen(img = ps_input, 
                     pan = pan_band, 
                     method = "pca")
ps_pca
## class      : RasterBrick 
## dimensions : 666, 666, 443556, 4  (nrow, ncol, ncell, nlayers)
## resolution : 15, 15  (x, y)
## extent     : 399997.5, 409987.5, 2799998, 2809988  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=44 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : pan_sharpening_input_1_pan, pan_sharpening_input_2_pan, pan_sharpening_input_3_pan, pan_sharpening_input_4_pan 
## min values :                   8721.067,                   7228.862,                   4734.816,                   9755.690 
## max values :                   17155.96,                   17871.61,                   18922.92,                   22283.65

El aumento de la resolución espacial del ráster original de 30 a 15 m puede observarse comparando la imagen en color verdadero (TCC) tanto del ráster de entrada como del ráster de salida.

Imagen en color verdadero (TCC) del raster original con una resolución de 30 m

ggRGB(ps_pca, r = 3, g = 2, b = 1, stretch ="lin") + 
  ggtitle("Imagen en color verdadero original") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10))         # Tamaño de las etiquetas de los ejes

Imagen en color verdadero del raster original tras ser sometido a un “pan sharpening” según el método PCA.

ggRGB(ps_pca, r = 3, g = 2, b = 1, stretch ="lin") + 
  ggtitle("Pan sharpening según PCA") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10))         # Tamaño de las etiquetas de los ejes

7.2 Método Brovey

El método de Brovey realiza el pan-sharpening al multiplicar cada píxel de los datos multiespectrales de entrada con la proporción de la intensidad del píxel pancromático correspondiente a la suma de las intensidades multiespectrales de los píxeles en las bandas roja, verde y azul. También se puede modificar para la banda de infrarrojo cercano, pero solo está disponible para las bandas azul, verde y roja en el software R. Esto se puede aplicar asignando el argumento method = "brovery" en la función panSharpen(). Los números de banda de las longitudes de onda rojas, verdes y azules también se asignan en el argumento de esta función. Esta función proporciona la salida ráster de nitidez panorámica solo para estas tres bandas.

ps_brovey <- panSharpen(img = ps_input, 
                        pan = pan_band, 
                        r = 3, g = 2, b = 1, 
                        method = "brovey")
ps_brovey
## class      : RasterBrick 
## dimensions : 666, 666, 443556, 3  (nrow, ncol, ncell, nlayers)
## resolution : 15, 15  (x, y)
## extent     : 399997.5, 409987.5, 2799998, 2809988  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=44 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : pan_sharpening_input_1_pan, pan_sharpening_input_2_pan, pan_sharpening_input_3_pan 
## min values :                   3270.699,                   3122.848,                   2884.092 
## max values :                   6348.279,                   6556.499,                   6888.222

Imagen en color verdadero del raster original tras ser sometido a un “pan sharpening” según el método Brovery.

ggRGB(ps_brovey, r = 3, g = 2, b = 1, stretch ="lin") + 
  ggtitle("Pan sharpening según método Brovery") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                  size =20),         # Tamaño del título
        axis.title = element_text(size =10))         # Tamaño de las etiquetas de los ejes

7.3 Método de saturación de tono de intensidad (IHS)

El método Intensity Hue Saturation (HIS) de enfoque pancromático convierte la imagen multiespectral de RGB a intensidad, tono y saturación. La mayor parte de la información espacial del compuesto RGB está presente en el componente de intensidad, mientras que los componentes de tono y saturación tienen información espectral. El histograma de la banda pancromática se empareja con el histograma del componente de intensidad. Luego, el componente de intensidad se reemplaza con una imagen pancromática y se aplica una transformación inversa para obtener una imagen RGB con detalles espaciales de la banda pancromática de mayor resolución espacial. Este método solo es aplicable a las bandas roja, verde y azul de la imagen ráster multiespectral en el software R. Esto se puede realizar asignando method = ihs, argumento y números de banda correspondiente a las longitudes de onda rojas, verdes y azules en la función panSharpen.

ps_ihs <- panSharpen(img = ps_input, pan =pan_band, r =3, g =2, b = 1, method = "ihs")

ggRGB(ps_ihs, r = 3, g = 2, b = 1, stretch ="lin") + 
  ggtitle("Pan sharpening según método Brovery") +                          # Título de la imagen 
  labs(x="Longitud(m)", y="Latitud (m)") +           # Etiquetas del eje
  theme(plot.title = element_text(hjust =0.5,        # Título centrado
                                size =20),         # Tamaño del título
        axis.title = element_text(size =10))         # Tamaño de las etiquetas de los ejes

ACTIVIDAD INDIVIDUAL