Análisis de dos indices de vegetación (NDVI y ARVI) en cuatro provincias del departamento de Cundinamarca, a partir de imágenes Landsat 7 para el año 2012
Nadia Luque Sanabria, Ing. Agronoma. Estudiante maestria en geomática. Universidad Nacional de Colombia
28.05.2021
1. Introducción
La percepción remota permite obtener información sobre la dinámica de la vegetación, proporcionando información muy útil para realizar el monitoreo de las diferentes áreas identificadas (Xue, J. & Su, B.,2017) A través de índices de vegetación, calculados a partir de valores de reflectividad en distintas longitudes de onda, se puede obtener información sobre el estado y las características de la vegetación (Carvacho, L. & Sanchez, M., 2010).
En este informe se muestra el estado de la vegetación y su distribución en un momento especifico, en una zona del departamento de Cundinamarca que presenta una gran variedad en sus paisajes y vegetación. Para este fin, se aplicaron dos índices de vegetación: el Índice de Vegetación de Diferencia Normalizada - NDVI (por su nombre en inglés Normalized Difference Vegetation Index) y el Índice de vegetación resistente a la atmosfera - ARVI (por su nombre en inglés Atmospherically resistant vegetation index), en cuatro provincias de Cundinamarca. De esta forma se identifica la distribución de la vegetación y su concordancia con la clasificación de cobertura del suelo realizada por el Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM) en el año 2010-2012 para Colombia.
1.1 Problema
El departamento de Cundinamarca presenta diversos pisos térmicos, y por ende, con gran variedad de vegetación. Los índices de vegetación se han estudiado y aplicado ampliamente en diferentes zonas, pero aún son pocos los informes o publicaciones dedicados a la zona de Cundinamarca. El interés principal es realizar observación sobre la distribución de la vegetación en cuatro provincias del departamento para el año 2012, a partir de dos índices de vegetación y como referencia los datos de cobertura del suelo de Colombia 2010-2012, proporcionados por el IDEAM.
1.2 Objetivo
El objetivo de este estudio es analizar y comparar dos índices de vegetación en una zona del departamento de Cundinamarca, a partir de imágenes Landsat 7 adquiridas en el 2012 y teniendo como referencia la cobertura de la tierra durante el periodo 2010-2012, basado en la interpretación visual de imágenes de satélite con verificación en campo.
1.3 Conceptos teóricos
Índice de Vegetación de Diferencia Normalizada - NDVI
El índice de vegetación normalizada NDVI es el más utilizado para monitorear la cobertura y producción de vegetación, es un índice que en plantas mide que tan verde es el área y de esta manera se estima la cantidad de vegetación que existe en un área específica y su estado (Meneses, C., 2011). Es sensible a los efectos de la atmosfera y el suelo.
NDVI = (NIR Band-Red Band)/(NIR Band+Red Band)
Índice de vegetación resistente a la atmosfera - ARVI
El índice de vegetación resistente a la atmosfera corrige los efectos de la dispersión atmosférica, combinando las bandas roja y azul, las cuales cuentan con las propiedades de autocorrección del efecto atmosférico. También mejora las superficies con presencia de partículas de aerosol de tamaño pequeño a moderado (Kaufman and Tanré, 1992).
ARVI = Near_Infrared−Red−y(Red−Blue)Near_Infrared+Red−y(Red−Blue)
Cobertura del suelo
La cobertura del suelo o de la tierra hace referencia a lo que se observa sobre la superficie de la tierra, discriminando entre vegetación, cuerpos de agua, afloramientos rocosos, centros urbanos, entre otras características que se puedan diferenciar en la superficie. La cobertura también se puede definir como unidades delimitadas a partir de respuestas espectrales determinadas por características en común (IDEAM, 1997).
2. Metodos
2.1 Zona de estudio
Cundinamarca tiene un área de 24.210 Km2 y está compuesto por 116 municipios divididos en 15 provincias (Figura 1): Almeidas, Alto Magdalena, Bajo Magdalena, Gualivá, Guavio, Magdalena Centro, Medina, Oriente, Rionegro, Sabana Centro, Sabana Occidente, Soacha, Sumapaz, Tequendama y Ubaté. Se ubica en la cordillera oriental y presenta variedad de topografía incluyendo relieves bajos, planos y montañosos, incluyendo pisos térmicos desde los 300 hasta más de 3.500 msnm, permitiendo gran variedad de climas (Gobernación de Cundinamarca, 2018).
library(tidyverse)
library(dplyr)
library(ggplot2)
library(sf)
library(RColorBrewer)
library(leaflet)
library(readxl)Cund = st_read("D:/Datos referencia SHP/Provincias de cundinamarca/CundProv_Proj.shp")
ggplot(Cund)+geom_sf()Figura 1. Provincias de Cundinamarca.
Este análisis se focalizó en 4 provincias de Cundinamarca: Alto Magdalena, Sabana Occidente, Tequendama y Soacha (figura 2). La provincia del alto Magdalena está dividida en ocho municipios, cuenta con un área de 1.187 Km2 y la temperatura promedio se encuentra entre 23 a 35°C.La provincia de Sabana occidente, cuenta con ocho municipios, una extensión de 914 km2. La provincia de Tequendama cuenta con diez municipios y un área de 1168 km2. Por último la provincia de Soacha, la más pequeña, cuenta con dos municipios y un área de 307 km2.
Cund = st_read("D:/Datos referencia SHP/Provincias de cundinamarca/Provincias_5_Cund.shp")
names(Cund)
Cund.gcs <- st_transform(Cund, "+proj=longlat +datum=WGS84")
st_crs(Cund.gcs)
m<-leaflet() %>%
addTiles() %>%
addPolygons(data = Cund.gcs ,stroke = TRUE, fillOpacity = 0.8,
highlight = highlightOptions(weight = 3,
color = "red",
bringToFront = TRUE))mFigura 2. Mapa de las provincias de Alto Magdalena, Sabana Occdiente, Tequendama y Soacha (área de estudio).
2.2 Descripción de datos
El Instituto de Hidrología, Meteorología y Estudios Ambientales - IDEAM cuenta con un geo portal en el que se encuentran datos abiertos a nivel nacional. Cuenta con mapa temático que representa la cobertura de la tierra en Colombia para el periodo 2010-2012 a escala 1:100.000, generado a partir de la metodología Corine Land Cover. Este documento se redujo a las cuatro provincias de Cundinamarca y cuenta con la siguiente información: tejido urbano continuo, tejido urbano discontinuo, zonas industriales y comerciales, aeropuertos, zonas de extracción minera, zonas de disposición de residuos, zonas verdes urbanas, instalaciones recreativas, otros cultivos transitorios, cereales, hortalizas, tubérculos, cultivos permanentes arbustivos, cultivos agroforestales, cultivos confinados, pastos limpios, pastos enmalezados, mosaico de cultivos, mosaico de pastos y cultivos, mosaico de cultivos, pastos y espacios naturales, mosaicos de pastos con espacios naturales, mosaico de cultivos con espacios naturales, bosque denso, bosque fragmentado, bosque de galería y ripario, plantación forestal, herbazal, arbustal, vegetación secundaria o en transición, afloramientos rocosos, tierras desnudas y degradadas, zonas quemadas, zonas pantanosas, vegetación acuática sobre cuerpos de agua, ríos, lagunas, lagos y ciénagas naturales, cuerpos de agua artificiales y nubes.
Para facilidad de interpretación y visualización de los datos de referencia, se combinaron algunas clases para reducir la lista que originalmente estaba incluida, quedando las siguientes clases: Arbustos, vegetación herbácea, agricultura, urbano, desnudo, agua permanente, humedal herbáceo, bosque cerrado, bosque abierto y nubes.
Adicionalmente, se descargó una imagen de Landsat 7 nivel dos para el año 2012, del portal Earth Explorer de Estados Unidos (USGS). En esta imagen aparece gran parte del departamento de Cundinamarca y se tiene despejada el área de las cuatro provincias en estudio. La imagen obtenida por Landsat 7 corresponde a la zona path 008 row 057 del año 2012. Se escogieron las cuatro provincias que estaban más despejadas en la imagen y que tenían completa la información de cobertura del suelo para los años 2010-2012.
3. Resultados
3.1 Imagen landsat 7 de 2012:
La imagen landsat 7 de 2012 está compuesta por las bandas: Banda 1= blue, banda 2=green , banda 3= red, banda 4= Near infrared (NIR), banda 5= Short-wave Infrared (SWIR), banda 6= Thermal infrarred (TIR) y banda 7= Mid-Infrared (MIR).
library(raster)
library(terra)
LS02 = stack ("D:/Datos referencia SHP/Image Landsat8/LE07_L2SP_008057_20120221_02_T1/sebkcstack.TIF")
names(LS02)
names(LS02) <- c('blue', 'green', 'red', 'NIR', 'SWIR', 'TIR', 'MIR')El primer acercamiento para analizar la imagen es conocer sus características descriptivas, con la cuál se obtiene información como: valores mínimos, máximos, medidas de tendencia central, medidas de dispersión e histogramas. Esta información permite explorar las imágenes y establecer los procesos que se deben desarrollar para mejorar las características de la imagen y realizar otro tipo de cálculos. Los valores obtenidos muestran, para la mayoría de las bandas, un valor maximo de 65635, que corresponde a una imagen de 16 bits. Aún se encuentra en valores digitales, por lo que se debe convertir a valores de reflectancia.
#Recorto raster con un shapefile de referencia
library(rgdal)
Prov4cund=readOGR(dsn='D:/Datos referencia SHP/Provincias de cundinamarca/Provincias_5_Cund.shp',layer='Provincias_5_Cund')## OGR data source with driver: ESRI Shapefile
## Source: "D:\Datos referencia SHP\Provincias de cundinamarca\Provincias_5_Cund.shp", layer: "Provincias_5_Cund"
## with 4 features
## It has 8 fields
## Integer64 fields read as strings: OBJECTID(LS02crop <- crop(LS02, extent(Prov4cund)))## class : RasterBrick
## dimensions : 3227, 2996, 9668092, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 511755, 601635, 464385, 561195 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : r_tmp_2021-05-28_160529_22312_97267.grd
## names : blue, green, red, NIR, SWIR, TIR, MIR
## min values : 6667, 6521, 6210, 7538, 6664, 34940, 6732
## max values : 65535, 65535, 65535, 65535, 65535, 54325, 65535writeRaster(LS02crop, filename="LS02crop_Prov4Cund.TIF", overwrite=TRUE)
summary(LS02crop)## blue green red NIR SWIR TIR MIR
## Min. 6829 6969 7001 7677.00 7271.00 35259.00 7190
## 1st Qu. 8284 8816 8502 16021.89 12988.75 43120.00 9831
## Median 8613 9283 9121 17989.75 14367.00 44340.33 10870
## 3rd Qu. 9423 10163 10447 19843.00 16322.04 45515.00 12812
## Max. 65535 65535 65535 65535.00 65535.00 53316.05 65535
## NA's 0 0 0 0.00 0.00 123462.00 03.2 Cálculo de reflectancia aparente TOA:
Para realizar el análisis de la imagen, es necesario transformar sus niveles digitales a valores de radiancia, teniendo en cuenta los coeficientes de calibración del sensor y un factor de escala. Posteriormente, los niveles se normalizan, llevándolos de energía reflejada a valores de reflectividad (Tristán, Wainschenker y Doorn, 2008), que expresan en porcentaje, la relación entre el flujo incidente y el reflejado por una superficie (Camargo, C., Pacheco, C., & López, R., 2021). Al realizar el cálculo cambian los valores de niveles digitales a valores entre -1 y 1. Para esta imagen, que es de nivel 2, el proceso ya se ha realizado antes de descargar la imagen, solo necesita un factor de conversión para ver los valores de reflectancia.
TOALS02crop2<- -0.2 + LS02crop*0.0000275
TOALS02crop2## class : RasterBrick
## dimensions : 3227, 2996, 9668092, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 511755, 601635, 464385, 561195 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : blue, green, red, NIR, SWIR, TIR, MIR
## min values : -0.0166575, -0.0206725, -0.0292250, 0.0072950, -0.0167400, 0.7608500, -0.0148700
## max values : 1.602213, 1.602213, 1.602213, 1.602213, 1.602213, 1.293938, 1.602213names(TOALS02crop2) <- c('blue', 'green', 'red', 'NIR', 'MIR', 'TIR', 'SWIR')
writeRaster(TOALS02crop2, filename="TOALS02crop_Prov4Cund2.TIF", overwrite=TRUE)
#https://www.usgs.gov/faqs/how-do-i-use-a-scale-factor-landsat-level-2-science-products?qt-news_science_products=0#qt-news_science_products3.3 Visualización color natural y falso color:
En la figura 3 se muestran dos imágenes: a) composición natural (r = 4, g = 3, b = 2) y b) falso color (r = 5, g = 4, b = 3). La imagen tiene presencia de aproximadamente el 40% de nubes. Para esta zona de Cundinamarca es un buen porcentaje, ya que, por su posición geográfica, varias épocas del año permanecen con un porcentaje mayor al 50% de nubes. En la imagen de falso color (b) se aprecian mejor los detalles que en la imagen natural (a), resaltando las masas de agua y el suelo descubierto.
par(mfrow=c(1,2))
plotRGB(TOALS02crop2, r = 3, g = 2, b = 1, main = "a.Composición natural", stretch="hist", colNA="white")
plotRGB(TOALS02crop2, r = 4, g = 3, b = 2, main = "b.Composición falso color", stretch="hist", colNA="white")
Figura 3. Imagen de landsat 7, recortada en el área de estudio a. Composición natural y b. composición falso color
3.4 Matriz de covarianza:
La matriz de covarianza contiene la relación lineal entre las bandas, es decir la covarianza y los valores de variación con respecto a la media. Se pueden leer a lo largo de la diagonal de la matriz. Al tener todos los valores positivos, indica que las bandas varían hacia la misma dirección, es decir, que cuando una banda aumenta la otra también. Entre las bandas del espectro visible no hay mucha covarianza, pero cuando se observan las bandas TIR y SWIR si se ve una diferencia entre los valores.
memory.limit(size=40000)## [1] 40000covmat <- cov(as.matrix(TOALS02crop2), use = "pairwise")
round(covmat*1000,2)## blue green red NIR MIR TIR SWIR
## blue 90.88 77.02 87.91 47.92 34.60 -8.96 29.73
## green 77.02 71.70 76.09 43.97 31.12 -7.66 26.39
## red 87.91 76.09 87.51 47.23 34.83 -8.42 29.96
## NIR 47.92 43.97 47.23 36.35 22.67 -4.60 17.44
## MIR 34.60 31.12 34.83 22.67 20.01 -2.22 16.00
## TIR -8.96 -7.66 -8.42 -4.60 -2.22 3.42 -2.39
## SWIR 29.73 26.39 29.96 17.44 16.00 -2.39 14.103.5 Matriz de correlación:
La correlación mide la fuerza y la dirección de la relación lineal entre dos bandas. Los valores de correlación se encuentran entre -1 y 1. Para éste caso, la mayoría se encuentra entre 0 y 1. Al ser estos positivos, indica que la media varía en la misma dirección para las dos bandas. La banda TIR es la única con valores negativos, es decir que cuando las otras bandas aumenta esta disminuye y viceversa. Que los valores estén tan cerca de uno (1) indican que existe una alta relación entre ellas, p.ej., las bandas del espectro visible entre si.
library(PerformanceAnalytics)
tmp1 <- as.matrix(sample(TOALS02crop2,500000))
chart.Correlation(tmp1, histogram=TRUE, pch=19)
Figura 4. Matriz de correlación entre las 7 bandas que componen la imagen.
3.6.Análisis de componentes principales:
El análisis de componentes principales tiene la ventaja de sintetizar la información de variación de “n” bandas en una, principalmente (PC1) y máximo tres (PC2 y PC3), que pueden recoger hasta el 90% de la variación de las bandas. Las nuevas bandas que se generan corresponden a la misma cantidad de bandas con las que se trabaja, en este caso, siete (Figura 5). El primer componente, PC1, representa la mayor proporción de varianza del conjunto de bandas estudiado.
library(RStoolbox)
library(reshape2)
library(factoextra)
LS02crop= rast('D:/Dropbox/NYO/Nadia/Maestria en geomatica/Remote sensing/Laboratorio1 remote sensing/LS02crop_Prov4Cund.TIF')
names(LS02crop) <- c('blue', 'green', 'red', 'NIR', 'SWIR', 'TIR', 'MIR')
TOALS02crop2= rast('D:/Dropbox/NYO/Nadia/Maestria en geomatica/Remote sensing/Laboratorio1 remote sensing/TOALS02crop_Prov4Cund2.TIF')
names(TOALS02crop2) <- c('blue', 'green', 'red', 'NIR', 'SWIR', 'TIR', 'MIR')
TOALS02crop_2= stack('D:/Dropbox/NYO/Nadia/Maestria en geomatica/Remote sensing/Laboratorio1 remote sensing/TOALS02crop_Prov4Cund2.TIF')
names(TOALS02crop_2) <- c('blue', 'green', 'red', 'NIR', 'SWIR', 'TIR', 'MIR')Se puede observar en la figura 6, que el primer componente, PC1, es el que mayor grado de variación contiene y al combinarlo con PC2 y PC3, se obtiene más información sobre lo que contiene la imagen de landsat 7.
set.seed(1)
sr <- spatSample(TOALS02crop2, 10000)
ACP <- prcomp(na.omit(sr), scale = TRUE)
ACP## Standard deviations (1, .., p=7):
## [1] 2.3533675 0.9205920 0.5296482 0.4931157 0.2259588 0.1706498 0.1015032
##
## Rotation (n x k) = (7 x 7):
## PC1 PC2 PC3 PC4 PC5 PC6
## blue 0.4094720 -0.06808091 0.3935754 -2.064305e-01 0.418734772 -0.031678817
## green 0.4095091 -0.03847383 0.3450824 -5.435309e-02 -0.833713098 -0.109284629
## red 0.4119139 -0.02915215 0.3546419 -2.330458e-01 0.328345527 0.001697604
## NIR 0.3821847 0.06709540 0.0150785 8.706077e-01 0.095915349 0.286111920
## SWIR 0.3857883 0.34513472 -0.4591008 3.937663e-02 0.074978181 -0.717010085
## TIR -0.2242105 0.89735256 0.3728429 5.574362e-05 -0.005839695 0.071975649
## MIR 0.3870321 0.25332493 -0.5002328 -3.749847e-01 -0.083170256 0.621213694
## PC7
## blue -0.6736766050
## green -0.0420153118
## red 0.7359226950
## NIR 0.0137157558
## SWIR -0.0003516556
## TIR -0.0161732978
## MIR -0.0486056721screeplot(ACP, col="blue")
Figura 5. Análisis de componentes principales
# función para restringir la predicción a los primeros dos componentes principales
pca_predict2 <- function(model, data, ...) {
predict(model, data, ...)[,1:2]
}
pci <- predict(TOALS02crop2, ACP, fun=pca_predict2)
plot(pci)
Figura 6. Representación de los dos primeros componentes, PC1 y PC2.
#Realizamos el raster con los 3 primeros componentes
pca <- rasterPCA(TOALS02crop_2)
PCA_123 = stack(pca$map$PC1, pca$map$PC2, pca$map$PC3)
writeRaster(PCA_123,"ACP_Cund4prov2012.tif", drivername="Gtiff", overwrite=TRUE)En la figura 7, se compara la imagen en falso color (Bandas R=4, G=3, B=2), con la imagen obtenida con los tres primeros componentes (Bandas R=PC1, G=PC2, B=PC3). Se observa en la imagen de componentes principales un mayor detalle y discriminación de lo que se encuentra en la imagen de falso color. En tonos blancos y amarillos resalta las nubes, en tonos verdes claros sugiere zonas con menor vegetación y mayor contenido de agua o humedad y en tonos negros correspondería a zonas urbanas.
par(mfrow = c(1,2))
plotRGB(TOALS02crop2,4,3,2, stretch = "hist", main = "Land7_2012 RGB432")
plotRGB(pca$map,1,2,3, stretch = "lin", main = "ACP (PC -123)")
Figura 7. Comparación imagen falso color con imagen de componentes principales
3.7 Transformación Tasseled Cap:
En la transformación Tasseled Cap se forman capas o bandas nuevas que tienen un significado físico con respecto a la imagen. En la figura 8 se observa la banda brightness, relacionada con el suelo descubierto o zona urbana. Lo primero que se resalta es la zona urbana de la provincia de sabana occidente, que es la zona urbana más grande de las cuatro provincias. En la parte baja se encuentra suelo descubierto, todo en color verde. La banda greenness revela una gran cantidad de vegetación en toda la zona de estudio. Lo que se observa en color amarillo o verde más claro, corresponde a suelo descubierto y zona urbana. La banda wetness, muestra la humedad y en la imagen parece ser homogéneo el contenido de humedad en la zona de estudio; hay una leve diferenciación donde hay nubes.
library(RStoolbox)
lsat_tc <- tasseledCap(TOALS02crop_2[[c(1:5,7)]], sat = "Landsat7ETM")
lsat_tc## class : RasterBrick
## dimensions : 3227, 2996, 9668092, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 511755, 601635, 464385, 561195 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : r_tmp_2021-05-28_161038_22312_55345.grd
## names : brightness, greenness, wetness
## min values : 0.01717249, -1.52384836, -1.90471245
## max values : 3.5705307, 0.3886665, 0.7499328plot(lsat_tc)
Figura 8. Diferentes capas de la transformación Tasseled Cap
3.8. Perfil espectral
El perfil espectral permite visualizar a partir de datos de referencia, en este caso polígonos de clasificación de cobertura del suelo, la información espectral que se encuentra de cada clase en la imagen en cada una de las bandas. Cada tipo de cobertura se caracteriza por reflejar la radiación de manera específica. Por ejemplo, la vegetación presenta valores altos en la banda del infrarrojo.
# Cargo los datos de referencia
library(sp)
Cund4prov <- shapefile('D:/Datos referencia SHP/Provincias de cundinamarca/CoverLand10-12_4Cund.shp')Para poder construir los perfiles, se debe tomar una muestra de valores de reflectancia para cada clase de paisaje en cada una de las bandas y se promedia la muestra para obtener un valor por banda y de esta manera generar el perfil espectral. Como es una muestra, puede que no tome valores de todas las clases, esto depende de la cantidad de datos que se tenga por clase y del tipo de muestra que se realice; que para este caso, es regular.
#Se toma la muestra de los poligonos de referencia para extraer los valores de reflactancia de la imagen
set.seed(1)
ptCund4prov <- spsample(Cund4prov, 200, type='regular')
ptCund4prov$class <- over(ptCund4prov,Cund4prov)$Clase
ptCund4prov <- vect(ptCund4prov)
xy <- as.matrix(geom(ptCund4prov)[,c('x','y')])
dfm <- extract(TOALS02crop2, xy)
head(dfm)set.seed(1)
pe <- aggregate(dfm, list(ptCund4prov$class), mean)
# en lugar de la primera columna, usamos nombres de fila
rownames(pe) <- pe[,1]
pe <- pe[,-1]
pe## blue green red NIR SWIR
## Agricultura 0.08268605 0.08736537 0.09142126 0.3086822 0.22167405
## Arbustos 0.02926292 0.03437028 0.02988778 0.1587971 0.08566185
## Bosque abierto 0.44985251 0.44938272 0.44658918 0.5528881 0.26646645
## Bosque cerrado 0.08472889 0.09325024 0.09083092 0.2816411 0.19198398
## Nubes 0.32879750 0.31678001 0.31742625 0.3853787 0.26717000
## Urbano 0.09521250 0.11612900 0.13733700 0.1980020 0.22949500
## Vegetación 0.03165336 0.04652780 0.04048763 0.2863603 0.18138541
## Vegetación herbácea 0.03657563 0.05499375 0.06019812 0.2094337 0.18979875
## TIR MIR
## Agricultura NaN 0.12514511
## Arbustos 0.9582323 0.04227144
## Bosque abierto 0.9414352 0.18509396
## Bosque cerrado 1.0019723 0.11540259
## Nubes 0.8934447 0.20654625
## Urbano 1.0823415 0.20602100
## Vegetación 1.0425871 0.08480023
## Vegetación herbácea 1.0237294 0.11381625Para la muestra que se realizó de la clase que contiene polígonos de cuerpo de agricultura, arbustos, bosque abierto, bosque cerrado, nubes, urbano, vegetación y vegetación herbácea, tomo solo ocho clases (figura 9). La banda que mejor las discrimina es la del infrarrojo cercano (4) y la del infrarrojo termal (6). Al ver el comportamiento de la reflectancia en la banda NIR (4) se aprecia que en esta banda se obtuvo una diferenciación entre las clases muy marcada, que muestra que el promedio de valores para cada clase está bien definido.
set.seed(1)
# Cree un vector de color para las clases de cobertura terrestre para usarlo en el trazado.
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
#transformar ps de un data.frame a una matriz
pe <- as.matrix(pe)
# Primero crea una grafica vacía
plot(0, ylim=c(0,1), xlim = c(1,7), type='n', xlab="Bandas", ylab = "Reflectancia")
#plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bandas", ylab = "Reflectancia")
# agregar las diferentes clases
for (i in 1:nrow(pe)){
lines(pe[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Title
title(main="Perfil espectral para cuatro provincias de Cundinamarca (2010-2012)", font.main = 2)
# Legend
legend("topleft", rownames(pe),
cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")
Figura 9. Perfil espectral para el área de estudio, con datos de referencia de cobertura del suelo proporcionados por el IDEAM para los años 2010-2012.
Se pueden usar algoritmos de árboles de decisión, como el árbol de regresión y clasificación CART. Los algoritmos en general son de carácter descriptivo, que muestran una segmentación jerárquica (Parra, 2019).
set.seed(1)
ptCund4prov1 <- spsample(Cund4prov, 100, type='random')
ptCund4prov1$class <- over(ptCund4prov1,Cund4prov)$Clase
ptCund4prov1 <- vect(ptCund4prov1)
xy <- as.matrix(geom(ptCund4prov1)[,c('x','y')])
dfm <- extract(TOALS02crop2, xy)
head(dfm)
sampdata <- data.frame(class = ptCund4prov1$class, dfm)library(rpart)
library(rpart.plot)
cartmodel <- rpart(as.factor(class)~., data = sampdata, method = 'class', minsplit = 5)
prp(cartmodel, faclen = 0, cex = 0.8, extra = 1)
Figura 10. Árbol de regresión y clasificación CART, para la zona de estudio, a partir de los datos de referencia de cobertura del suelo 2010-2012
En el árbol para la muestra de datos obtenida, la banda de mayor relevancia es la TIR. Cuando la reflectancia de esta banda es menor a 0.097 y la banda roja es menor a 0.029 se encuentra agricultura. Si TIR es mayor a 0.097, la banda green es mayor a 0.061, la banda MIR menor a 0.089 y la banda blue menor a 0.057 se encuentra bosque cerrado. Sin embargo, la cantidad de datos y clasificaciones puede ser poca, o se pueden encontrar más datos en una sola clasificación que en otra, tal vez por eso la mayoría de terminaciones muestra a la clase agricultura.
3.9. Indice de vegetación NDVI y ARVI
Los resultados obtenidos en el índice NDVI se relacionan con el dosel, la radiación fotosintéticamente activa, índice de área foliar, evapotranspiración y biomasa vegetal (Borowik, T., et al 2013). Se diferencia en el rango de -0.4 a 0 suelo descubierto, masas de agua, nubes y zonas urbanas. La vegetación se encuentra en los rangos mayores a 0.6, siendo la vegetación más densa y sana cercana a 1.
El índice de vegetación resistente a la atmosfera ARVI, hace un proceso de autocorrección del efecto atmosférico en la banda roja, obteniendo una nueva combinación de banda rojo azul que minimiza los efectos de la dispersión atmosférica causada por aerosoles en la banda roja. ARVI tiene similitudes con NDVI, pero es menos sensible a los efectos atmosféricos (Abdou, B. et al 1996). La imagen obtenida es muy similar a la de NDVI, sin embargo, se diferencia un poco mejor el suelo descubierto.
par(mfrow=c(1,2))
vi <- function(img, k, i) {
bk <- img[[k]]
bi <- img[[i]]
vi <- (bk - bi) / (bk + bi)
return(vi)
}
ndvi <- vi(TOALS02crop2,4,3)
ndvi[ndvi>1] <- 1; ndvi[ndvi< (-1)] <- -1
plot(ndvi, col=rev(terrain.colors(10)), main = "NDVI")
rvi <- function(img, k, i, e) {
bk <- img[[k]] #NIR
bi <- img[[i]] #Red
be <- img[[e]] #Blue
rvi <- (bk - (2*bi)+be) / (bk + (2*bi)+be)
return(rvi)
}
arvi <- rvi(TOALS02crop2, 4, 3, 1)
arvi[arvi>1] <- 1; arvi[arvi< (-1)] <- -1
plot(arvi, col=rev(terrain.colors(10)), main = "ARVI")
Figura 11. Representación de índices de vegetación NDVI y ARVI para la zona de estudio.
La matriz de correlación (Figura 12) confirma que existe una fuerte relación entre estos dos índices y que ambos se mueven en la misma dirección, la tendencia es lineal, aunque se ven unos datos por fuera de la línea de tendencia.
indices = c(ndvi,arvi)
names(indices) <- c('NDVI','ARVI')
pairs(indices)
Figura 12. Matriz de correlación para los índices NDVI y ARVI en el área de estudio.
3.11 Histograma NDVI y ARVI
El histograma permite visualizar cual es la distribución de los datos para cada uno de los índices. Los histogramas de los dos índices (Figura 13) son bimodales, ya que presentan sutilmente unos picos dentro de la gráfica. En el histograma de NDVI se aprecian cuatro picos, uno entre -0.5 y -0.3, el otro entre -0.2 y 0.05, el tercero entre 0.1 y 0.2, el cuarto entre 0.3 y 1 aproximadamente. El histograma del índice ARVI presenta tres picos uno entre -0.3 y -0.1, el segundo entre 0 y 0.05, el tercero entre 0.2 y 1 aproximadamente. Estos picos pueden estar relacionados con clases bien diferenciadas, que podrían estar asociadas a vegetación y zona urbana con suelo descubierto.
par(mfrow=c(1,2))
hist(ndvi,
main = "Distribucion de valores NDVI",
xlab = "NDVI",
ylab= "Frecuencia",
col = "wheat",
xlim = c(-0.6, 1),
breaks = 30,
xaxt = 'n', maxcell=1000000)
axis(side=1, at = seq(-0.6,1, 0.05), labels = seq(-0.6,1, 0.05))
hist(arvi,
main = "Distribucion de valores ARVI",
xlab = "ARVI",
ylab= "Frecuencia",
col = "wheat",
xlim = c(-0.6, 1),
breaks = 30,
xaxt = 'n', maxcell=1000000)
axis(side=1, at = seq(-0.6,1, 0.05), labels = seq(-0.6,1, 0.05))
Figura 13. Histogramas para índices NDVI y ARVI.
4. Discusión
Las provincias de Alto Magdalena, Sabana Occidente, Tequendama y Soacha, capturadas en la imagen de landsat 7 para el año 2012, muestran un porcentaje de nubes menor al 40%, permitiendo realizar la observación de algunas características de la zona. En la imagen natural y falso color, muestra hacia el centro de la imagen una zona montañosa que en su gran parte presenta vegetación. Se discriminan zonas urbanas y suelo desnudo de la vegetación, aunque es muy poco detalle el que se puede observar con solo estas dos imágenes. Al realizar el análisis de componentes principales, se ve con mejor detalle lo que se encuentra en la imagen, siendo mas claras las zonas urbanas, algunos cuerpos de agua y suelo desnudo. Las nubes y la vegetación en las tres imágenes son evidentes.
Para la tranformación Tasseled CAP, las bandas greenness y wetness no diferencian bien lo que se aprecia en la imagen de componentes principales. Ya que en greenness no se diferencia bien el suelo desnudo de la vegetación, en wetness no son tan claros los cuerpos de agua. Al consultar los datos de referencia proporcionados por el Ideam, las clases que tienen mayor representación son las relacionadas con la vegetación, seguidas por la clase urbana. En cuanto a la clase agua permanente son muy pocos los polígonos identificados en los datos de referencia, esto podría explicar porque no son tan evidentes en la capa wetness.
Al observar el perfil espectral y el árbol de decisiones CART, se identifica que predominan las clasificaciones relacionadas con vegetación, seguido de las clases urbano y nubes. El de mayor predominio es Agricultura, que se diferencia rápidamente en el árbol generado por CART, con la banda TIR y la banda roja.
Los índices de vegetación NDVI y ARVI (figura 14), muestran diferencias en la vegetación sana, es decir, la que se encuentra en el umbral cercano a 1, siendo mayor para el índice NDVI. Probablemente el índice ARVI, al realizar las correcciones identifique mejor el suelo desnudo y las áreas que no tienen vegetación densa. Teniendo en cuenta que la clasificación de referencia es en su mayoría agricultura y esta es bastante diversa dependiendo del cultivo, el estado fenológico y sus asociaciones entre diferentes especies, puede ser que el índice ARVI lo discrimine mejor.
par(mfrow=c(1,2))
vn <- c(-Inf, 0.5,NA, 0.6, 1, 1, 1, Inf,NA)
rcl1 <- matrix(vn, ncol = 3, byrow = TRUE)
vegn <- classify(ndvi, rcl1)
plot(vegn, col=colorRampPalette(c("white", "darkgreen"))(255), main = "Vegetación NDVI")
va <- c(-Inf, 0.5, NA, 0.6, 1, 1, 1, Inf, NA)
rcl2 <- matrix(va, ncol = 3, byrow = TRUE)
vega <- classify(arvi, rcl2)
plot(vega,col=colorRampPalette(c("white", "darkgreen"))(255), main = "Vegetación ARVI")
Figura 14. Comparación umbral de vegetación (0.6 - 1) para índices NDVI y ARVI.
5. Conclusiones
Al comparar los dos índices en estudio, se aprecia que ARVI discrimina mejor la vegetación y el suelo desnudo que NDVI. Para las cuatro provincias se encuentra que la clase con mayor frecuencia es vegetación y según la clasificación de cobertura, está relacionada en su mayoría con agricultura. Sin embargo, sería bueno estudiar un área más grande que permitiera abarcar más puntos de muestreo para cada clase o discriminar mejor las clases, para mejorar su representación y análisis.
El índice de vegetación NDVI y el índice ARVI se comportaron de manera similar. Para este caso, en el que queremos explorar como se encontraba la vegetación de acuerdo con los índices y a los datos de referencia, se recomienda hacerlo con el análisis ARVI que discriminó mejor la vegetación del suelo desnudo.
Se recomienda completar dicha comparación de índices con una clasificación supervisada y no supervisada, mejorando las clases y sus respectivas muestras, calculando su exactitud, de esta manera se podrá hacer un análisis completo.
Referencias
Abdou, B., Morin, D., Bonn, F. & Huete, A. (1996). A review of vegetation indices. Remote Sensing Reviews. 13. 95-120. 10.1080/02757259509532298.
Borowik, T., Pettorelli, N., Sönnichsen, L., & Jedrzejewska, B. (2013). Normalized difference vegetation index (NDVI) as a predictor of forage availability for ungulates in forest and field habitats. European Journal of Wildlife Research. 59. 10.1007/s10344-013-0720-0.
Camargo, C.E. Pacheco, C.E. & López, R. (2021). Evaluación de métodos de corrección atmosférica y sombreado topográfico en imagen Landsat 8 OLI sobre un área montañosa semiárida. UD y la Geomática (16), 23-39. DOI: https://doi.org/10.14483/23448407.17040
IDEAM (1997). Instituto de hidrologia, meterología y estudios ambientales. Geosistemas de alta Montaña Colombiana. Bogotá: Universidad Nacional de Colombia.
Kaufman, YJ, Tanre, D (1992).Atmospherically resistant vegetation index (ARVI) for EOS-MODIS.IEEE Transactions on Geoscience and Remote Sensing, 30(2), 261-270, doi: 10.1109/36.134076.
Lizarazo, I, 2021. A tutorial on land cover classification from satellite imagery using R. Available at https://rpubs.com/ials2un/cart_landsat8.
Meneses, C., 2011. NDVI as indicator of degradation.Food and Agriculture Organization (FAO), Unasylva 238, Vol. 62.
Parra, F. 2019. Estadística y Machine Learning con R. Libro electronico. https://econometria.wordpress.com/2019/10/29/estadistica-y-machine-learning-con-r-libro-electronico/.
Tristan, P., Wainschenker, R. y Doorn, J. (2008). Normalización de imágenes satelitales en el análisis multitemporal. En X Workshop de Investigadores en Ciencias de la Computación (pp. 315–319). General Pico, Argentina: UNLPAM.
6. Reproducibilidad
sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_Europe.1252 LC_CTYPE=English_Europe.1252
## [3] LC_MONETARY=English_Europe.1252 LC_NUMERIC=C
## [5] LC_TIME=English_Europe.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rpart.plot_3.0.9 rpart_4.1-15
## [3] factoextra_1.0.7 reshape2_1.4.4
## [5] RStoolbox_0.2.6 PerformanceAnalytics_2.0.4
## [7] xts_0.12.1 zoo_1.8-9
## [9] rgdal_1.5-23 terra_1.2-10
## [11] raster_3.4-10 sp_1.4-5
## [13] readxl_1.3.1 leaflet_2.0.4.1
## [15] RColorBrewer_1.1-2 sf_0.9-8
## [17] forcats_0.5.1 stringr_1.4.0
## [19] dplyr_1.0.6 purrr_0.3.4
## [21] readr_1.4.0 tidyr_1.1.3
## [23] tibble_3.1.2 ggplot2_3.3.3
## [25] tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-1 ellipsis_0.3.2 class_7.3-19
## [4] fs_1.5.0 rstudioapi_0.13 proxy_0.4-25
## [7] ggrepel_0.9.1 prodlim_2019.11.13 fansi_0.4.2
## [10] lubridate_1.7.10 xml2_1.3.2 codetools_0.2-18
## [13] splines_4.1.0 doParallel_1.0.16 knitr_1.33
## [16] jsonlite_1.7.2 pROC_1.17.0.1 caret_6.0-88
## [19] broom_0.7.6 dbplyr_2.1.1 rgeos_0.5-5
## [22] compiler_4.1.0 httr_1.4.2 backports_1.2.1
## [25] assertthat_0.2.1 Matrix_1.3-3 cli_2.5.0
## [28] htmltools_0.5.1.1 tools_4.1.0 gtable_0.3.0
## [31] glue_1.4.2 Rcpp_1.0.6 cellranger_1.1.0
## [34] vctrs_0.3.8 nlme_3.1-152 iterators_1.0.13
## [37] crosstalk_1.1.1 timeDate_3043.102 gower_0.2.2
## [40] xfun_0.23 rvest_1.0.0 lifecycle_1.0.0
## [43] XML_3.99-0.6 MASS_7.3-54 scales_1.1.1
## [46] ipred_0.9-11 hms_1.1.0 parallel_4.1.0
## [49] yaml_2.2.1 geosphere_1.5-10 stringi_1.6.1
## [52] highr_0.9 foreach_1.5.1 e1071_1.7-6
## [55] lava_1.6.9 rlang_0.4.11 pkgconfig_2.0.3
## [58] evaluate_0.14 lattice_0.20-44 recipes_0.1.16
## [61] htmlwidgets_1.5.3 tidyselect_1.1.1 plyr_1.8.6
## [64] magrittr_2.0.1 R6_2.5.0 generics_0.1.0
## [67] DBI_1.1.1 pillar_1.6.1 haven_2.4.1
## [70] withr_2.4.2 units_0.7-1 survival_3.2-11
## [73] nnet_7.3-16 modelr_0.1.8 crayon_1.4.1
## [76] KernSmooth_2.23-20 utf8_1.2.1 rmarkdown_2.8
## [79] grid_4.1.0 data.table_1.14.0 ModelMetrics_1.2.2.2
## [82] reprex_2.0.0 digest_0.6.27 classInt_0.4-3
## [85] stats4_4.1.0 munsell_0.5.0 quadprog_1.5-8