1. Introducción

El departamento de Cundinamarca es muy diverso, ya que presenta una variedad en sus paisajes y vegetación. Conocer el estado de la vegetación y su distribución, permite realizar seguimiento y evaluación de las características y el comportamiento de dicha vegetación en un área específica. La percepción remota permite obtener información sobre la dinámica de la vegetación, proporcionando información muy útil para realizar 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). Para este fin se aplicarán 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 atmósfera - ARVI (por su nombre en inglés Atmospherically Resistant Vegetation Index), en una zona de Cundinamarca. De esta forma podremos identificar la distribución de la vegetación y si está acorde a la clasificación del suelo realizada en el año 2001 para el departamento de Cundinamarca.

El departamento de Cundinamarca presenta diversos pisos térmicos, lo que conlleva a contar 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 realizar observación sobre la distribución de la vegetación en el departamento, a partir de dos índices de vegetación y clasificación de suelos del departamento de Cundinamarca para 2001.

Objetivo

El objetivo de este estudio es analizar y comparar dos índices de vegetación en el departamento de Cundinamarca, a partir de imagenes landsat 7 adquiridas en el 2002 y teniendo como referencia una clasificación de uso del suelo realizada en el año 2001, que suministra información importante acerca del recurso suelo, que sirve como base para la determinación de las potencialidades y limitaciones de uso del suelo.

2. Metodos

Zona de estudio

El área de estudio de este informe es el departamento de Cundinamarca que tiene un área de 24.210 Km2 y está compuesto por 116 municipios divididos en 15 provincias: 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")
names(Cund)
ggplot(Cund) +
  geom_sf()

Para este análisis, se focalizo en 4 provincias Alto MAgadalena, Sabana Occdiente, Tequendama y Soacha. La provincia del alto Magdalena esta dividida en ocho municipios, cuenta con un área de 1.187 Km2 y la temepratura 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 mas 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))
m

Descripción de datos

El Instituto geográfico Agustín Codazzi cuenta con un geo portal en el que se encuentran datos abiertos a nivel nacional. Para Cundinamarca cuenta con mapa temático que representa la distribución de las características del suelo a escala 1:100.000 realizado en el año 2001. Este documento cuenta con información de paisaje, clima, tipo de relieve e información del suelo como litología, características, componentes, perfil y porcentaje de componentes. Cuando se recortó la imagen para las 4 provincias de interés, con respecto a las clases de paisaje que se trabajaran en el análisis, quedaron las siguientes: cuerpo de agua, montaña, pantano, planicie, valle, zona urbana.

Adicionalmente, se descargó una imagen de Landsat 7 para el año 2002, donde aparece gran parte del departamento de Cundinamarca y se tiene despejada el área de las cuatro provincias en estudio.

3. Resultados y discusión

3.1 Imagen landsat 2002:

La imagen obtenida en Landsat 7 corresponde a la zona 008057 del año 2002, que toma gran parte del departamento de Cundinamarca. Se escogieron las 4 provincias que estaban más despejadas en la imagen y que tenían completa la información de características del suelo para el año 2001. Las bandas de la imagen compuesta con las que se trabajaron son las siguientes: Banda 1= blue, banda 2= green, banda3= red, banda4= infrarrojo cercano (NIR), banda 5= infrarrojo medio rango 1.55 a 1.75 (MIR) y banda 7= infrarrojo medio rango 2.08 a 2.35 (SWIR).

library(raster)
library(terra)
LS02 = stack ("D:/Datos referencia SHP/Image Landsat8/LE07_L2SP_008057_20030101/LE07_L2SP_08057_200301_Rast6Band1.TIF")
names(LS02)
names(LS02) <- c('blue', 'green', 'red', 'NIR', 'MIR', 'SWIR')
summary (LS02)

Las imágenes obtenidas a través de satélites deben ser procesadas para obtener información apropiada de ellas. El primer acercamiento a dicha información es la estadística descriptiva, en la que obtenemos información como valores mínimos, máximos, medidas de tendencia central, medidas de dispersión e histogramas. Esta información nos permite explorar las imágenes y conocer los datos con los que contamos, que procesos debemos desarrollar para mejorar las características de la imagen y realizar otro tipo de procesos.

Para esta imagen tenemos un valor máximo de 65535 pixeles para las seis bandas, lo que corresponde a una imagen de 16 bits.

#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, Prov4cund))
## class      : RasterBrick 
## dimensions : 3227, 2996, 9668092, 6  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 511752.7, 601632.7, 464391.4, 561201.4  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +ellps=intl +units=m +no_defs 
## source     : C:/Users/nyluq/AppData/Local/Temp/RtmpkFEkuE/raster/r_tmp_2021-04-13_165136_35060_90721.grd 
## names      :  blue, green,   red,   NIR,   MIR,  SWIR 
## min values :  6796,  6924,  6823,  7426,  6994,  6755 
## max values : 65535, 65535, 65535, 65535, 65535, 65535
writeRaster(LS02crop, filename="LS02crop_Prov4Cund.TIF", overwrite=TRUE)
summary(LS02crop)
##          blue    green   red   NIR   MIR  SWIR
## Min.     6817  6924.00  6910  7571  7346  7166
## 1st Qu.  8080  8618.00  8325 15514 12619  9607
## Median   8375  9032.00  8902 17262 14008 10575
## 3rd Qu.  8841  9635.25  9820 18812 15677 12124
## Max.    65535 65535.00 65535 65535 65535 65535
## NA's        0     0.00     0     0     0     0

3.2 Calculo 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. En segundo lugar, 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 grandes 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.

TOALS02crop<-   -0.2 + LS02crop*0.0000275
TOALS02crop
## class      : RasterBrick 
## dimensions : 3227, 2996, 9668092, 6  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 511752.7, 601632.7, 464391.4, 561201.4  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +ellps=intl +units=m +no_defs 
## source     : memory
## names      :       blue,      green,        red,        NIR,        MIR,       SWIR 
## min values : -0.0131100, -0.0095900, -0.0123675,  0.0042150, -0.0076650, -0.0142375 
## max values :   1.602213,   1.602213,   1.602213,   1.602213,   1.602213,   1.602213
writeRaster(TOALS02crop, 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_products

3.3 Visualización color natural y falso color:

Tenemos un primer acercamiento a la imagen obtenida, observando sus características en composición natural (r = 3, g = 2, b = 1) y falso color (r = 4, g = 5, b = 3). La imagen tiene presencia de un 20% de nubes aproximadamente, pero para esta zona de Cundinamarca es un buen porcentaje, ya que, por su posición geográfica, varias épocas del año permanecen un porcentaje mayor de nubes. en la imagen de falso color se aprecian un poco mejor los detalles que en la imagen natural, resaltando un poco más las masas de agua y el suelo descubierto.

par(mfrow=c(1,2))
plotRGB(TOALS02crop, r = 3, g = 2, b = 1,  main='Composición natural', stretch='lin',  colNA='white')
plotRGB(TOALS02crop, r = 4, g = 5, b = 3,  main='Composición falso color', stretch='hist', colNA='white')

### 3.4 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, nos 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 varianza, pero cuando se observan las bandas NIR y SWIR si se ve una diferencia entre los valores.

memory.limit(size=40000)
## [1] 40000
covmat <- cov(as.matrix(TOALS02crop), use = "pairwise")
round(covmat*1000,2)
##        blue green   red   NIR   MIR  SWIR
## blue  43.53 37.75 41.70 15.44 25.61 17.42
## green 37.75 36.26 37.39 14.62 24.55 16.42
## red   41.70 37.39 41.69 15.18 26.01 17.77
## NIR   15.44 14.62 15.18 12.30 12.74  8.05
## MIR   25.61 24.55 26.01 12.74 23.80 15.20
## SWIR  17.42 16.42 17.77  8.05 15.20 12.81

3.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 nuestro caso, todos se encuentra entre 0 y 1, al ser todos positivos indica que la media varía en la misma dirección para las dos bandas. Al estar los valores tan cerca de uno, como por ejemplo las bandas del espectro visible entre si, quiere decir que existe una alta relación entre ellas. Para el caso de la banda NIR la relación con las demás bandas no es tan alta.

library(PerformanceAnalytics)
tmp1 <- as.matrix(sample(TOALS02crop,500000))
chart.Correlation(tmp1, histogram=TRUE, pch=19)

3.6.Análisis de componentes principales:

El análisis de componentes principales tiene la ventaja de sintetizar 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 seis. El primer componente denominado PC1 es el que representa la mayor proporción de varianza del conjunto de bandas estudiado.

library(RStoolbox)
library(reshape2)
library(factoextra)
LS02crop1= rast('D:/Maestria en geomatica/Remote sensing/Laboratorio1 remote sensing/LS02crop_Prov4Cund.TIF')
TOALS02crop1= rast('D:/Maestria en geomatica/Remote sensing/Laboratorio1 remote sensing/TOALS02crop_Prov4Cund2.TIF')

Podemos observar en la siguiente grafica que el primer componente, que en este caso corresponde a la banda PC1, es el que mayor grado de variación contiene y al combinarlo con PC2 y PC3, podemos llegar a obtener un poco mas de información sobre lo que contiene la imagen de landsat 7.

set.seed(1)
sr <- spatSample(TOALS02crop1, 10000)
ACP <- prcomp(na.omit(sr), scale = TRUE)
ACP
## Standard deviations (1, .., p=6):
## [1] 2.2424207 0.6691936 0.5958256 0.3347635 0.2029121 0.1244241
## 
## Rotation (n x k) = (6 x 6):
##             PC1        PC2        PC3         PC4          PC5         PC6
## TOA_1 0.4241340  0.3613298 -0.2336144  0.15808599 -0.456231611 -0.63390608
## TOA_2 0.4284333  0.2835549 -0.1848761 -0.05229875  0.830775558 -0.09454724
## TOA_3 0.4290365  0.3418544 -0.1555180  0.10351645 -0.278878308  0.76576016
## TOA_4 0.3574311 -0.7809173 -0.4763900  0.18667745 -0.009441656  0.02293748
## TOA_5 0.4177229 -0.1829132  0.2899267 -0.82821854 -0.144497546 -0.03416585
## TOA_6 0.3875997 -0.1699524  0.7589825  0.49067666  0.054063784 -0.03379110
screeplot(ACP)

# 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(TOALS02crop1, ACP, fun=pca_predict2)
plot(pci)

#Realizamos el raster con los 3 primeros componentes
pca <- rasterPCA(TOALS02crop)
PCA_123 = stack(pca$map$PC1, pca$map$PC2, pca$map$PC3)
writeRaster(PCA_123,"ACP_Cund4prov2001.tif", drivername="Gtiff", overwrite=TRUE)

A continuación, comparamos un ráster de 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). Lo que se alcanza a observar en la imagen de componentes principales, es un mayor detalle y discriminación de lo que se encuentra en la imagen falso color. En tonos fucsia resalta las nubes, en tonos verdes pareciera las zonas con menor vegetación y mayor contenido de agua o humedad, en tonos negros lo que podría corresponder a zonas urbanas y en diferentes tonos azules la vegetación. Dependiendo de la intensidad del tono azul, pareciera que es la densidad de la vegetación, entre más oscuro más denso.

par(mfrow = c(1,2))
plotRGB(TOALS02crop1,4,3,2, stretch = "lin", main = "Land7_2001 B543")
plotRGB(pca$map,1,2,3, stretch = "lin",  main = "ACP (PC -123)")

3.7 Transformación Tasseled Cap:

En dicha transformación se forman capas o bandas nuevas que tienen un significado fisico con respecto a la imagen. En la siguiente imagen observamos brightness, que corresponde a la banda TC1 y esta relacionada con el suelo descubierto o zona urbana. Lo primero que se resalta es parte de la nube que capturo la imagen, pero hacia al lado izquierdo si se puede diferenciar una parte mas oscura que corresponde a zona urbana. La banda greenness que corresponde a TC2 y esta relacionada con la vegetación, revela una gran cantidad de vegetación en la zona de estudio. Lo que se observa en color amarillo o verde mas claro, corresponde a suelo descubierto, zona urbana y la nube de la imagen.La banda wetness (TC3), muestra la humedad y en la imagen alcanza a capturar la humedad de la nube y por eso no se diferencia tan claramente como en las dos bandas anteriores.

lsat_tc <- tasseledCap(TOALS02crop[[c(1:6)]], sat = "Landsat7ETM")
lsat_tc
## class      : RasterBrick 
## dimensions : 3227, 2996, 9668092, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 511752.7, 601632.7, 464391.4, 561201.4  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +ellps=intl +units=m +no_defs 
## source     : C:/Users/nyluq/AppData/Local/Temp/RtmpkFEkuE/raster/r_tmp_2021-04-13_165708_35060_55345.grd 
## names      :   brightness,    greenness,      wetness 
## min values :  0.003841463, -1.916818053, -2.000451413 
## max values :    3.5705306,    0.4046172,    0.8196520
plot(lsat_tc)

3.8. Perfil espectral

El perfil espectral nos permite ver a partir de datos de referencia, en este caso polígonos de clasificación de paisaje, la información espectral que encontramos de cada clase en la imagen en cada una de las bandas. Cada tipo de superficie 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/Cund_Suelos_2001_4provin_2.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 tengamos por clase y del tipo de muestra que se realice, que para nuestro caso es regular.

#Se toma la muestra de los poligonos de referencia para extraer los valores de reflactancia de la imagen
ptCund4prov <- spsample(Cund4prov, 100, type='regular')
ptCund4prov$class <- over(ptCund4prov,Cund4prov)$PAISAJE
ptCund4prov <- vect(ptCund4prov)
xy <- as.matrix(geom(ptCund4prov)[,c('x','y')])
dfm <- extract(TOALS02crop, xy)
head(dfm)
##          blue     green       red       NIR       MIR      SWIR
## [1,] 0.055750 0.0788500 0.0922975 0.2602125 0.2655200 0.1542000
## [2,] 0.022750 0.0321825 0.0263525 0.2693975 0.1516425 0.0536875
## [3,] 0.085285 0.1192200 0.1675375 0.2862275 0.4349475 0.3559950
## [4,] 0.020715 0.0383150 0.0258850 0.3407600 0.1617350 0.0616625
## [5,] 0.089685 0.1072300 0.1160300 0.2429425 0.2249850 0.1515325
## [6,] 0.062295 0.0824800 0.1162775 0.2461050 0.3523925 0.2166250
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

De esta manera obtenemos el perfil espectral, para la muestra que realizamos de la clase paisaje que contiene polígonos de cuerpo de agua, montaña, pantano, planicie, valle, zona urbana, tomo solo cinco clases. El perfil que corresponde a cuerpo de agua se comporta muy similar a lo que se encuentra en la literatura, que la reflectancia es prácticamente 0 en las bandas RGB. Las clases de Montaña, planicie y valle tienen un comportamiento similar entre ellas, pero la banda que mejor las discrimina es la banda del infrarrojo cercano (4) y tal vez la banda roja (3). La zona urbana es casi lineal, se refleja de manera muy similar en todas las bandas. Al ver comportamiento de la reflectancia en la banda NIR (4), tal vez explica porque es mayor su variación y es por que obtenemos en esta banda una sutil diferenciación entre las clases, pero muy marcada, que muestra que el promedio de valores para cada clase está bien definido.

# 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,0.4), xlim = c(1,6), 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 (2002)", font.main = 2)
# Legend
legend("topleft", rownames(pe),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")

Podemos 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).

ptCund4prov1 <- spsample(Cund4prov, 100, type='random')
ptCund4prov1$class <- over(ptCund4prov1,Cund4prov)$PAISAJE
ptCund4prov1 <- vect(ptCund4prov1)
xy <- as.matrix(geom(ptCund4prov1)[,c('x','y')])
dfm <- extract(TOALS02crop, 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)

En el árbol para la muestra de datos obtenida, la banda de mayor relevancia es la banda blue. Cuando la reflectancia de la banda blue es menor a 0.027 encontramos montaña. Si blue es mayor a 0.027, adicionalmente la banda MIR es menor a 0.14 encontramos planicie. Sin embargo, la cantidad de datos y clasificaciones puede ser poca y tal vez por eso la mayoría de terminaciones muestra a la clase montaña.

3.9. Indice de vegetación 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). Los resultados obtenidos de este índice 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). Para las cuatro provincias de Cundinamarca, en su mayoría es vegetación. 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.2, siendo la vegetación más densa y sana mayor a 0.5.

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}
ndvi <- vi(LS02crop1, 4, 3)
plot(ndvi, col=rev(terrain.colors(10)), main = "NDVI")

##3.10. Indice de vegetación ARVI

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 la vegetación.

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(LS02crop1, 4, 3, 1)
plot(arvi, col=rev(terrain.colors(10)), main = "ARVI")

La matriz de correlación nos confirma que existe una fuerte relación entre estos dos índices y que ambos se mueven hacia la misma dirección, la tendencia es lineal, aunque se ven unos datos por fuera de la línea de tendencia. Para este caso parece que cualquiera de los índices nos da la misma información.

indices <- c(ndvi,arvi)
names(indices) <- c('ndvi', 'arvi')
pairs(indices)

3.11 Histograma NDVI y ARVI

par(mfrow=c(1,2))
hist(ndvi,
     main = "Distribucion de valores NDVI",
     xlab = "NDVI",
     ylab= "Frecuencia",
     col = "wheat",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n', maxcell=5000000)
## Warning: [hist] 10% of the raster cells were used. (of which 59% was NA)
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

hist(arvi,
     main = "Distribucion de valores ARVI",
     xlab = "ARVI",
     ylab= "Frecuencia",
     col = "wheat",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n', maxcell=5000000)
## Warning: [hist] 10% of the raster cells were used. (of which 59% was NA)
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

En los histogramas encontramos similitud entre los dos índices, aunque para NDVI tiene un rango más amplio y agrupa datos entre -0.5 y -0.2, en ARVI el mínimo valor es desde -0.2.

4. Conclusiones

El índice de vegetación NDVI y el índice ARVI se comportaron de manera similar. Para este caso en el que queremos ver como se encontraba la vegetación de acuerdo con los índices y a los datos de referencia, se puede analizar con uno de los dos índices. En la imagen obtenida del 2001, la mayor cobertura es vegetativa y lo encontramos dentro de la clase montaña y valle.

Los datos de referencia de suelos para las cuatro provincias de Cundinamarca no tienen mucho detalle con respecto a la discriminación de la vegetación, pero la información permitió realizar los perfiles espectrales y confirmar si lo que nos muestra los índices de vegetación está de acuerdo con lo que se obtuvo de los rangos del espectro para cada clase.

Obtener la información descriptiva y estadística de la imagen, es información valiosa que aporta al análisis de los datos.

Referencias

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.

Meneses, C., 2011. NDVI as indicator of degradation.Food and Agriculture Organization (FAO), Unasylva 238, Vol. 62.

Abdou, B., Morin, D., Bonn, F. & Huete, A. (1996). A review of vegetation indices. Remote Sensing Reviews. 13. 95-120. 10.1080/02757259509532298.

IGAC. 2001. Datos abiertos subdirección de agrología. Mapas de Suelos del Territorio Colombiano a escala 1:100.000.Departamento: Cundinamarca. https://geoportal.igac.gov.co/contenido/datos-abiertos-agrologia.

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

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.

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/.

Lizarazo, I, 2021. A tutorial on land cover classification from satellite imagery using R. Available at https://rpubs.com/ials2un/cart_landsat8.