LABORATORIO 6. TRANSFORMACIONES MATEMÁTICAS: ÍNDICES ESPECTRALES
G174: FOTOINTERPRETACIÓN Y TELEDETECCIÓN
Domingo F. Rasilla
2025-04-14
💡 OBJETIVOS:
Comprender el concepto y la importancia de los índices espectrales.
Crear índices a partir de bandas de imágenes de satélite
Evaluar los cambios en el territorio utilizando las diferencias entre los índices en diferentes momentos del tiempo.
💡 MATERIALES PARA LA ACTIVIDAD:
Los materiales para el desarrollo de esta actividad son los siguientes.
► Imágenes satélite. Es la misma imagen multibanda utilizada en la actividad previa.
INTRODUCCION
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…
DESARROLLO DE LA ACTIVIDAD
En primer lugar, deben importarse todos los paquetes necesarios para realizar este ejercicio:
library(terra) # Importación y transformación de imágenes ráster## Warning: package 'terra' was built under R version 4.3.3## terra 1.8.21library(RStoolbox) # Transformación y visualización de datos ráster.## Warning: package 'RStoolbox' was built under R version 4.3.3## This is version 1.0.0 of RStoolboxlibrary(ggplot2) # Visualización de datos ráster con RStoolbox## Warning: package 'ggplot2' was built under R version 4.3.3library(rasterVis) # Visualización de datos ráster## Warning: package 'rasterVis' was built under R version 4.3.2## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.3.2En primer lugar se deberá establecer el directorio de trabajo.
En el presente apartado trabajaremos con el objeto
imagen que ha sido preparado en la introducción a este
apartado. Si no se dispusiera de él deberá importarse otra vez.
imagen <- rast("D:/G174_2025/LABORATORIO_6_Transformaciones_matematicas/datos/imagen.tif")
names(imagen)## [1] "coastal aerosol" "blue" "green" "red"
## [5] "NIR" "SWIR1" "SWIR2"Existen múltiples í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
En R existen varias posibilidades de cálculo para obtener los diferentes índices espectrales.
MÉTODO 1: CÁLCULO DIRECTO
Este método usa 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
en la memoria de la máquina, debe asegurarse que el volumen de datos sea
manejable por la RAM de su máquina.
La fórmula para el cálculo del índice NDVI es la siguiente:
\[NDVI = \frac{NIR-RED}{NIR+RED}\]
En el caso de Landsat 5
\[NDVI (Landsat 5) = \frac{(B4 – B3)} {(B4 + B3)}\] En el caso de Landsat 8/9
\[NDVI (Landsat 8/9) = \frac{(B5 – B4)} {(B5 + B4)}\]
En el caso de Sentinel 2
\[NDVI (Sentinel 2) = \frac{(B8 – B4)} {(B8 + B4)} \]
Un ejemplo de cálculo empleando los nombres de las capas entre corchetes es el siguiente:
ndvi <- (imagen[["NIR"]] - imagen[["red"]]) / (imagen[["NIR"]] + imagen[["red"]])
ndvi## class : SpatRaster
## dimensions : 1245, 1497, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source(s) : memory
## varname : imagen
## name : NIR
## min value : -0.959918
## max value : 0.883263Una vez calculado el índice NDVI, se puede proceder a una
representación gráfica mediante la función
terra::plot().
plot(ndvi,
main = "Índice NDVI")
Como alternativa gráfica se sugiere el uso de la función
RStoolbox::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"), # 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
scale_fill_gradientn(name = "valores NDVI", colours = c("white", "darkgreen"))+ # Dibuja un gradiente de colores
theme_bw()
✅ SINTAXIS ALTERNATIVA:
Como se ha señalado anteriormente, existe una sintaxis alternativa,
por ejemplo, sustituyendo el nombre de la banda por su posición dentro
del objeto imagen. En caso de una escena de Landsat 8/9 que
contenga las bandas 1 a 7, la fórmula sería la siguiente:
ndvi_2 <- (imagen[[5]] - imagen[[4]]) / (imagen[[5]] + imagen[[4]])También se podría utilizar una sintaxis basada en el nombre (reducido) de las bandas, por ejemplo, si así las hubiéramos llamado:
ndvi_3 <- (imagen[[“B5”]] - imagen[[“B4”]]) / (imagen[[“B5”]] + imagen[[“B4”]])
Para comprobar que todos estos procedimientos resultan en objetos ráster iguales, se procede a su representación gráfica.
par(mfrow = c(2,2))
plot(ndvi, main = "NDVI nombres de capas")
plot(ndvi_2, main = "NDVI índices de capas")
par(mfrow = c(1,1))
A continuación podemos borrar este último objeto.
rm(ndvi_2)📝 ACTIVIDAD 1 ÍNDICES ESPECTRALES:
Calcula y representa gráficamente los siguientes índices espectrales
orientados hacia el análisis de la vegetación al objeto ráster
multibanda imagen. .
- Índice SAVI (Soil Adjusted Vegetation Index). Este índice es similar al NDVI, pero suprime los efectos de los píxeles que contienen suelo desnudo, incorporando un factor de ajuste (canopy background adjusting factor) denominado L. L es una función de la densidad de la vegetación, que puede basarse en el conocimiento previo su magnitud. En ausencia de tal conocimiento, Huete (1988) sugiere utilizar un valor de L = 0,5 para tener en cuenta las variaciones del fondo del suelo de primer orden. Este índice se utiliza mejor en áreas con vegetación relativamente escasa donde el suelo es visible. Su formulación es la siguiente:
\[SAVI = \frac{(L+1) * (NIR – RED)} {(NIR + RED + L)}\]
- Índice MSAVI (modified soil-adjusted vegetation index). Este índice de vegetación también aborda algunas de las limitaciones del NDVI cuando se aplica a áreas con suelos expuestos. El problema con el índice de vegetación ajustado al suelo (SAVI) original es que requería especificar el factor de corrección del brillo del suelo (L) mediante prueba y error, basado en la cantidad de vegetación en el área de estudio o usando un factor de corrección fijo de 0,5. Para calcular de manera más confiable y sencilla un factor de corrección de brillo del suelo, Qi et al. (1994a) desarrollaron MSAVI (Qi et al. 1994b).
\[MSAVI = NIR + 0.5 - (0.5 * (\sqrt{2NIR+1)^{2}}-8 * (NIR - (2*RED)))\]
- Índice EVI (Enhanced Vegetation Index). Este índice también se desarrolló para abordar algunas de las limitaciones del NDVI, en este caso para áreas con abundante biomasa, lo que causa un problema para el NDVI. Reduce la influencia de las condiciones atmosféricas en los valores del índice de vegetación y corrige las señales de fondo del dosel. En consecuencia, el EVI tiende a ser más sensible que el NDVI a las diferencias del dosel de las plantas, como el índice de área foliar (LAI), la estructura del dosel y la fenología de las plantas.
Los valores propuestos para los subíndices son:
L_evi= 1
C1 = 6
C2 = 7.5
G = 2.5
Índice EVI2 (two-band Enhanced Vegetation Index). Desarrollado recientemente sin la reflectancia de la banda azul como indicador de la fenología, cantidad y actividad de la vegetación. EVI2 tiene varias ventajas sobre NDVI, incluida la capacidad de resolver diferencias de LAI para vegetación con diferentes valores de reflectancia del suelo de fondo.
Para seguir con la actividad se borrarán los siguientes componentes
rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI)## Warning in rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI): objeto 'C1' no
## encontrado## Warning in rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI): objeto 'C2' no
## encontrado## Warning in rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI): objeto 'EVI' no
## encontrado## Warning in rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI): objeto 'EVI2' no
## encontrado## Warning in rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI): objeto 'G' no
## encontrado## Warning in rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI): objeto 'L' no
## encontrado## Warning in rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI): objeto 'L_evi' no
## encontrado## Warning in rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI): objeto 'MSAVI' no
## encontrado## Warning in rm(C1, C2, EVI, EVI2, G, L, L_evi, MSAVI, SAVI): objeto 'SAVI' no
## encontradoMÉTODO 2: CÁLCULO A TRAVÉS DE FUNCIONES.
Si la tarea consiste en obtener un mismo índice de múltiples imágenes, es conveniente crear una función y replicarla posteriormente tantas veces como imágenes. La función puede ser la siguiente:
funcion <- function(x, y) {
(x - y) / (x + y)
}La función lapp()[https://rspatial.github.io/terra/reference/lapp.html]
aplica cualquier función a todas las capas de un objeto
SpatRaster, o a subconjunto de ese objeto.
ndvi2 <- lapp(c(imagen[["NIR"]], imagen[["red"]]),
fun = funcion)
ggR(ndvi2, 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"), # 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
scale_fill_gradientn(name = "valores NDVI", colours = c("red", "yellow", "green", "darkgreen"))+ # Dibuja un gradiente de colores
theme_bw()
Una variante puede ser la siguiente:
ndvi3 <- lapp(imagen[[5:4]], fun = funcion)Como no se van a utilizar posteriormente, estos objetos se eliminan.
rm(nir, red, ndvi2, ndvi3, funcion)## Warning in rm(nir, red, ndvi2, ndvi3, funcion): objeto 'nir' no encontrado## Warning in rm(nir, red, ndvi2, ndvi3, funcion): objeto 'red' no encontradoMÉ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 multibanda (si se calculan simultáneamente varios índices).
La lista de índices calculados por la función
spectralIndices() 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 | Dash And Curran, 2004 | red, redEdge1, redEdge2 | (redEdge2 - redEdge1)/(redEdge1 - red) |
| 1NBRI | Normalised Burn Ratio Index | Garcia, 1991 1 nir, swir3 | (nir - swir3)/(nir + swir3) | |
| NDREI1 | Normalised Difference Red Edge Index 1 | Gitelson And Merzlyak, 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 | Guyot And Barnet, 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 |
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)”.
Argumentos
red = "red",nir = "nir",blue = "blue"yswir2 = "swir1"ya que se requieren las bandas del rojo, nir y swir2 (la sustitución del argumentoswir1porswir2se 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).Argumento
indices = "SLAVI".
► Por ejemplo, el índice NDVI en Landsat 8 se calcularía así, dado que en B4 = Red y B5 = NIR.
ndvi2 <- spectralIndices(imagen,
red = 4, nir = 5,
indices = "NDVI")
ndvi2## class : SpatRaster
## dimensions : 1245, 1497, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source(s) : memory
## name : NDVI
## min value : -0.959918
## max value : 0.883263► El cálculo del índice MSAVI sería de la siguiente manera
msavi <- spectralIndices(imagen,
red = 4, nir = 5,
indices = "MSAVI")
msavi## class : SpatRaster
## dimensions : 1245, 1497, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source(s) : memory
## name : MSAVI
## min value : -0.4613330
## max value : 0.7377701Algunos índices enumerados anteriormente requieren coeficientes adicionales.
Además, es posible calcular todos los índices espectrales disponibles
en la función spectralIndices() simultáneamente,
proporcionando todos los argumentos y los coeficientes requeridos. Como
no tenemos los valores de los coeficientes para ‘swir2ccc’ y ‘swir2coc’,
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 = 2, green = 3, red = 4, nir = 5, swir2 = 6, swir3 = 7)La salida de esta función es un objeto multibanda en el que cada banda 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"Podría consultarse los atributos de un índice específico (capa
ráster) escribiendo el nombre de la capa precedido del símbolo
$.
SI$NDVI # Propiedades de la capa ráster correspondiente## class : SpatRaster
## dimensions : 1245, 1497, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source(s) : memory
## name : NDVI
## min value : -0.959918
## max value : 0.883263Para comprobar que ambos procedimientos resultan en sendos objetos ráster iguales, se procede a su representación gráfica.
par(mfrow = c(1,2))
plot(ndvi, main = "NDVI")
plot(SI$NDVI, main = "NDVI según RStoolbox")
par(mfrow = c(1,1))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)",
par.settings=BuRdTheme())
Todos los índices señalados anteriormente están normalizados, es decir, sus límites máximo y mínimo son (1, -1). A partir de este gráfico se puede detectar que alguno(s) de los índices posee valores superiores a esos límites. Para detectar cuales, se puede recurrir a calcular sus valores extremos.
minmax(SI)## CTVI DVI EVI EVI2 GEMI GNDVI KNDVI
## min -0.6781725 -0.1207067 -0.6921978 -0.9936968 -2.257909 -0.9733608 0.0000000
## max 1.1761220 0.6577497 0.9999864 1.0000000 1.054443 0.7963593 0.7265791
## MNDWI MSAVI MSAVI2 NBRI NDVI NDWI NDWI2
## min -0.6122068 -0.4613330 -0.1833715 -0.6914897 -0.959918 -0.7963593 -0.7273799
## max 0.9960576 0.7377701 0.8439642 0.9918455 0.883263 0.9733608 0.9866543
## NRVI RVI SATVI SAVI SLAVI SR
## min -0.883263 0.06198658 -0.3632643 -0.2401769 0.02063482 0.02045087
## max 0.959918 48.89768126 0.3275060 0.7925146 3.71737512 16.13252462
## TTVI TVI WDVI
## min 0.0001877213 0.01141001 -0.1207067
## max 1.1761220057 1.17612201 0.6577497Por ejemplo, este es el caso del índice SR. Por ello, es recomendable verificar la distribución de sus valores dibujando su histograma (aunque se emplea sólo una muestra del total de celdas disponibles). Para llamarlo se menciona su nombre acompañado del signo $.
hist(SI$SLAVI,
main = "Distribución del índice SLAVI", # Title of histograma
xlab ="SLAVI", ylab = "Frecuencia", # Etiquetas de los ejes
col ="wheat") # Colores de los bins## Warning: [hist] a sample of 54% of the cells was used
La función terra::clamp() ajusta los valores de un
raster a sendos umbrales, máximo y mínimo, definidos por el investigador
(en este caso, -1 y 1). Además, el argumento values = TRUE
convierte los valores fuera de rango a valores equivalentes a los
umbrales.
SI$SLAVI <- clamp(SI$SLAVI, -1, 1, values = TRUE)
SI$SLAVI## class : SpatRaster
## dimensions : 1245, 1497, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 594090, 639000, 4190190, 4227540 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source(s) : memory
## name : SLAVI
## min value : 0.02063482
## max value : 1.00000000Al mismo tiempo, se puede dibujar un único índice espectral de ese raster multibanda
ggR(SI$SLAVI, 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"), # 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 = "Índice SLAVI", colours = c("skyblue", "wheat", "darkgreen"))
Para su posterior análisis, los índices NDVI, MNDWI(Normalized Difference Water Index) y NBRI (Normalized Difference Burn Index) fueron convertidos en objetos ráster
writeRaster(SI$NDVI, "D:/G174_2025/LABORATORIO_6_Transformaciones_matematicas/NDVI.tif", filetype = "GTiff", overwrite = TRUE)
writeRaster(SI$MNDWI, "D:/G174_2025/LABORATORIO_6_Transformaciones_matematicas/MDWI.tif", filetype = "GTiff", overwrite = TRUE)
writeRaster(SI$NBRI, "D:/G174_2025/LABORATORIO_6_Transformaciones_matematicas/NBRI.tif", filetype = "GTiff", overwrite = TRUE)📝 ACTIVIDAD 2 ÍNDICES ESPECTRALES:
Cada alumno deberá elegir dos índices espectrales que sean adecuados para el tragajo final que haya elegido.
A continuación, deberá buscar en internet 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.
📝 ACTIVIDAD 3 ÍNDICES ESPECTRALES:
El siguiente fichero comprimido contiene varios ficheros ráster.
El fichero L7 corresponde a una imagen multibanda del
satélite Landsat 7. Esta imagen multibanda consta de 6 bandas como las
del siguiente cuadro:
| Band | Nombre |
|---|---|
| B1 | Blue |
| B2 | Green |
| B3 | Red |
| B4 | Near-Infrared (NIR) |
| B5 | Shortwave Infrared 1 |
| B7 | Shortwave Infrared 2 |
► Carga la imagen L7.
► Cambia los nombres de las bandas para facilitar los posibles cálculos
► Visualiza las bandas en forma de una imagen visible y otra en falso color.
► Busca en internet el índice que responde al acrónimo NDBI (Built-up Index), procede a su cálculo sobre la escena anterior y represéntalo gráficamente en un mapa. ¿Qué utilidad tiene dicho índice? ¿Existe algún valor numérico sensible?
► Busca en internet el índice que responde al acrónimo NDWI, procede a su cálculo sobre la escena anterior y represéntalo gráficamente en un mapa. ¿Qué utilidad tiene dicho índice? ¿Existe algún valor numérico sensible?
► Busca en internet el índice que responde al acrónimo BSI (barren soil index). Su fórmula es:
\[BSI = \frac{(SWIR2 + RED)−(NIR + BLUE)} {(SWIR2 + RED)+(NIR + BLUE)}\]
► Calcúlalo con la función RStoolbox::spectralIndices()
y represéntalo gráficamente en un mapa. ¿Qué utilidad tiene dicho
índice? ¿Existe algún valor numérico sensible?
📝 ACTIVIDAD 4 ÍNDICES ESPECTRALES:
El fichero L8.tif contiene una escena de Landsat 8, con
las primeras 7 bandas. Calcula los siguientes índices espectrales y
represéntalos gráficamente
► Normalized Difference Moisture Index (NDMI):
► Moisture Stress Index (MSI):
📝 ACTIVIDAD 5 ÍNDICES ESPECTRALES:
El fichero vancouver_S2_29_07_2020.tif se encuentra
dentro de la continene una imagen multiespectral del área que rodea la
ciudad de Vancouver, en Canadá, tomada por el sensor Sentinel-2 MSI, el
día 29 de julio de 2020. Este fichero contiene 10 bandas con una
resolución de 20 m.
| Banda | Denominación |
|---|---|
| Banda 1 | Blue |
| Banda 2 | Green |
| Banda 3 | Red |
| Banda 4 | Red Edge 1 |
| Banda 5 | Red Edge 2 |
| Banda 6 | Red Edge 3 |
| Banda 7 | NIR |
| Banda 8 | NIR (Narrow) |
| Banda 9 | SWIR1 |
| Banda 10 | SWIR2 |
A partir de esas bandas:
► Convierte el fichero *.tif en un objeto
SpatRaster.
► Renombra las bandas de la manera que consideres más oportuna.
► Busca en internet la fórmula y calcula a continuación los
siguientes índices espectrales (todos relacionados con la vegetación)
con RStoolbox::spectralIndices():
Normalized Difference Vegetation Index (NDVI):
Green Normalized Difference Vegetation Index (GNDVI):
Enhanced Vegetation Index (EVI):
Advanced Vegetation Index (AVI):
Soil Adjusted Vegetation Index (SAVI):
Green Coverage Index (GCI):
📝 ACTIVIDAD 6 ÍNDICES ESPECTRALES:
► Crea una función para el cálculo del Índice de Vegetación Mejorado
(EVI) para el objeto ráster multibanda imagen. Esta función
debe incluir un término para identificar el objeto ráster sobre el que
se aplicará directamente el procedimiento matemático
correspondiente.
► Una vez creada la función, se aplica proporcionando a R el nombre y el número de la banda
► Represéntala gráficamente con la función
Rstoolbox::ggR().