Análisis Multitemporal: Humedal Juan Amarillo

Base de datos

library(readxl)

datos <- read_xlsx("datos_fin.xlsx")
datos 

El área del Humedal Juan Amarillo para el año 2012 es de 222,76 hectáreas (Fajardo López, 2012 - LINK), se puede observar que las áreas halladas con ArcGis varían entre 167.15 Ha y 250.05 Ha, cálculos que está sujetos a errores humanos ya que la delimitación de las mismas se realizaron de forma manual. Temperaturas promedio varían de 13.8 °C a 14.4 °C con unas temperaturas mínimas de 7.9 °C y unas máximas de 21.9 °C

Variación temporal

library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(ggpubr)
library(tidyr)

ggarrange(
  datos %>% 
  ggplot(data = ., aes(x = Fecha, y = Area_Ha)) +
  geom_point(size = 3, pch = 19) +
  geom_line() +
  ylim(0, 260) +
  labs(x = "Año", y = "Área (Ha)",
       title = "Variación temporal del área") +
  theme_bw() ,
  
  datos %>% 
  ggplot(data = ., aes(x = Fecha, y = poblacion/1e+06)) +
  geom_point(size = 3, pch = 19) +
  geom_line() +
  geom_text(aes(label = round(poblacion/1e+06, digits = 4)), nudge_y = 45000/1e+06) +
  labs(x = "Año", y = "Personas (millones)",
       title = "Población Bogotá") +
  theme_bw(),
  
  
  datos %>%
  select(Fecha, tem_pro, tem_min, tem_max) %>% 
  gather(key = "medida", value = "temp", -Fecha) %>% 
  ggplot(data = ., aes(x = Fecha, y = temp, color = medida)) +
  geom_point(size = 3) +
  geom_line() +
  scale_color_brewer(palette = "Set1") +
  labs(x = "Año", y = "°C",
       title = "T. máximas, mínimas y promedio",
       subtitle = "Estación El Dorado", color = "") +
  theme_bw() +
  theme(legend.position = "bottom"),
  
  datos %>%
  select(Fecha, dias_lluvia, dias_tormen) %>% 
  gather(key = "medida", value = "dias", -Fecha) %>% 
  ggplot(data = ., aes(x = Fecha, y = dias, color = medida)) +
  geom_point(size = 3) +
  geom_line() +
  scale_color_brewer(palette = "Set1") +
  labs(x = "Año", y = "°C",
       title = "Días de lluvia y tormenta",
       subtitle = "Estación El Dorado", color = "") +
  theme_bw() +
  theme(legend.position = "bottom"),
  
  ncol = 2, nrow = 2
)

Como ya se mencionó, el cálculo del área se realizó de forma manual, por lo tanto, puede presentar algunos sesgos en comparación con la literatura consultada, de igual forma, se puede apreciar que el área el humedal va disminuyendo a medida que pasan los años, fenómeno que se viene presentando de manera significativa desde los años 50 (Camargo, 2016). Se puede apreciar un aumento significativo de la población según bases de datos tomadas del DANE (Departamento Administrativo Nacional de Estadística) con relación a los años, valores que concuerdan con el Análisis demográfico y proyecciones poblacionales de Bogotá realizado en marzo del año 2018 por la secretaría distrital de planeación LINK Temperaturas máximas y mínimas no presentan una variación significativa, a excepción el año 2019, esto se puede presentar gracias a que aún no se posee los datos de los mese restantes, por lo que el promedio puede variar un poco más con respecto a los otros años. El comportamiento entre los días de lluvias y tormentas puede presentar muchas variaciones para cada año, ya que la fuente no posee los datos para todos los días

Asociación de variables

g1 <- datos %>%
  ggplot(data = ., aes(x = poblacion/1e+06, y = Area_Ha)) +
  geom_point(size = 4) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", size = 1.5) +
  labs(x = "Personas (millones)", y = "Área (ha)",
       subtitle = "Población vs Área") +
  theme_bw()

g2 <- datos %>%
  ggplot(data = ., aes(x = tem_pro, y = Area_Ha)) +
  geom_point(size = 4) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", size = 1.5) +
  labs(x = "°C", y = "Área (ha)",
       subtitle = "Temperatura promedio vs Área") +
  theme_bw()

g3 <- datos %>%
  ggplot(data = ., aes(x = tem_min, y = Area_Ha)) +
  geom_point(size = 4) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", size = 1.5) +
  labs(x = "°C", y = "Área (ha)",
       subtitle = "Temperatura mínima vs Área") +
  theme_bw()

g4 <- datos %>%
  ggplot(data = ., aes(x = tem_max, y = Area_Ha)) +
  geom_point(size = 4) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", size = 1.5) +
  labs(x = "°C", y = "Área (ha)",
       subtitle = "Temperatura máxima vs Área") +
  theme_bw()

g5 <- datos %>%
  ggplot(data = ., aes(x = dias_lluvia, y = Area_Ha)) +
  geom_point(size = 4) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", size = 1.5) +
  labs(x = "Días", y = "Área (ha)",
       subtitle = "Días con lluvia vs Área") +
  theme_bw()

g6 <- datos %>%
  ggplot(data = ., aes(x = dias_tormen, y = Area_Ha)) +
  geom_point(size = 4) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", size = 1.5) +
  labs(x = "Días", y = "Área (ha)",
       subtitle = "Días con tormenta vs Área") +
  theme_bw()

ggarrange(g1, g2, g3, g4, g5, g6, ncol = 3, nrow = 2)

Las anteiores gráficas están sujetas a un tamaño de muestra de 6 (los seis años analizados: 2014, 2015, 2016,2017,2018, 2019) las cúales pueden mostrar relaciones similares a la realidad aumentando el tamaño de la muestra.

Regresión lineal

modelo <- lm(Area_Ha ~ tem_pro, data = datos)
summary(modelo)

## 
## Call:
## lm(formula = Area_Ha ~ tem_pro, data = datos)
## 
## Residuals:
##       1       2       3       4       5       6 
##  10.043  32.103 -37.759 -18.013  -5.747  19.375 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1137.31     744.97   1.527    0.202
## tem_pro       -65.20      52.95  -1.231    0.286
## 
## Residual standard error: 28.68 on 4 degrees of freedom
## Multiple R-squared:  0.2749, Adjusted R-squared:  0.09357 
## F-statistic: 1.516 on 1 and 4 DF,  p-value: 0.2856

El modelo, aunque no es signifivativo (p=0.2856) permite inferir que por cada grado que aumenta la temperatura anualmente, el humedal disminuye aproximadamente 65.20 hectáreas. Esta relación es netamente matemática (no se presenta en la vida real) por lo tanto, esta sujeta a validaciones con un tamaño de muestra superior.

El modelo ajustado se puede expresar de la siguiente manera:

\[área = 1137.31 -\ 65.20*temperatura\]

Lectura de información espacial en R del año 2014

• Lectura de bandas para el color natural (Bandas 4,3,2)

# Cargando bibliotecas
library(sp)
library(raster)
library(landsat8)
library(tidyverse)

# Función brick() multibanda
# Función raster() una banda
tif_nat14 <- brick(x = "FINAL/imagenes/2014/inflarojo_verdadero/color_verdadero2014.tif")
plot(tif_nat14)

• Lectura de bandas para el color infrarrojo (Bandas 5,4,3)

tif_inf14 <- brick(x = "FINAL/imagenes/2014/inflarojo_verdadero/infrarrojo2014.tif")
plot(tif_inf14)

• Lectura del shape de área realizado en el ArcGis

# Lectura de .shp
shape_ja <- shapefile("FINAL/Humedal_2014_01_01.shp")
plot(shape_ja)

La lectura del formato .TIF muestra las bandas que componen cada imagen por separado, lo cual facilita la observación de cada una.

Banda 4 (Color natural)

# Cuatro gráficos en uno
par(mfrow = c(2, 2))

# Una banda
plot(tif_nat14[[1]],
     main = "Banda 4")

plot(tif_nat14[[1]],
     main = "Banda 4\nDos colores",
     col = terrain.colors((2)))

plot(tif_nat14[[1]],
     main = "Banda 4\nCuatro colores",
     col = terrain.colors((4)))

plot(tif_nat14[[1]],
     main = "Banda 4\nSeis colores",
     col = terrain.colors((6)))

A medida que se va aumentando el número de colores en cada banda, se puede apreciar una diferenciación más clara de cada color, lo que permite la observación de direntes cobertura vegetales.

Banda 4,3,2 (Color natural)

# 3 gráficos en uno
par(mfrow = c(2, 2))

# Todas las bandas
plot(tif_nat14[[1]],
     main = "Banda 4")

plot(tif_nat14[[2]],
     main = "Banda 3")

plot(tif_nat14[[3]],
     main = "Banda 2")

Se puede apreciar la lectura de las bandas que se combinan para dar paso al el color natural de la imagen, la cúal involucra tres bandas visibles (2,3,4), a cada banda se le asigna su verdadero color y cuando se combinan, genera una imágen con el color real de la escena. LINK

Banda 5 de 2014 (Infrarrojo)

# Cuatro gráficos en uno
par(mfrow = c(2, 2))

# Una banda
plot(tif_inf14[[1]],
     main = "Banda 5")

plot(tif_inf14[[1]],
     main = "Banda 5\nDos colores",
     col = terrain.colors((2)))

plot(tif_inf14[[1]],
     main = "Banda 5\nCuatro colores",
     col = terrain.colors((4)))

plot(tif_inf14[[1]],
     main = "Banda 5\nSeis colores",
     col = terrain.colors((6)))

A medida que se van sumando colores a la banda, se puede apreciar de una manera más clara la diferenciación entre coberturas vegetales, se puede además ver que se puede apreciar mejor la diferencia entre coberturas en el área urbana cuando se hace el paso de cuatro colores a seis colores.

Bandas 5, 4 y 3 (Color infrarrojo)

# 3 gráficos en uno
par(mfrow = c(2, 2))

# Todas las bandas
plot(tif_inf14[[1]],
     main = "Banda 5")

plot(tif_inf14[[2]],
     main = "Banda 4")

plot(tif_inf14[[3]],
     main = "Banda 3")

Al asignarle el color rojo a la banda 4 (infrarroja cercana) esta banda va a tener ciertas similitudes con la combinación RGB 4,3,2 . Sin embargo, al dar más peso a la región infrarroja (bandas 4 y 5) se ve realzada la diferencia de humedad en suelos y vegetales. Generalmente cuanto mayor es la humedad del suelo más oscuro aparecerá éste. LINK

Distribución de niveles digitales (ND)

• Bandas 4, 3 y 2:

par(mfrow = c(1, 3))
# tres gráficos en uno
hist(tif_nat14[[1]],
     xlab = "Frecuencia", ylab = "Nivel Digital (ND)",
     main = "Distribución de ND\nBanda 4")

hist(tif_nat14[[2]],
     xlab = "Frecuencia", ylab = "Nivel Digital (ND)",
     main = "Distribución de ND\nBanda 3")


hist(tif_nat14[[3]],
     xlab = "Frecuencia", ylab = "Nivel Digital (ND)",
     main = "Distribución de ND\nBanda 2")

LINK

Cada ND (Nivel digital) representa un área de la tierra en un lugar específico. El valor de ND asignado al píxel esta en función de la luz reflejada o emitida de la porción de la superficie terrestre LINK,en el histograma de la banda dos, algunos son valores negativos, aunque son pocos, pueden indicar daños en la imagen o en los pixeles, el histograma de los niveles digitales puede resultar muy útil a la hora de visualzar valores extremos.

• Bandas 5, 4 y 3:

par(mfrow = c(1, 3))
# tres gráficos en uno
hist(tif_inf14[[1]],
     xlab = "Frecuencia", ylab = "Nivel Digital (ND)",
     main = "Distribución de ND\nBanda 5")

hist(tif_inf14[[2]],
     xlab = "Frecuencia", ylab = "Nivel Digital (ND)",
     main = "Distribución de ND\nBanda 4")

hist(tif_inf14[[3]],
     xlab = "Frecuencia", ylab = "Nivel Digital (ND)",
     main = "Distribución de ND\nBanda 3")

Área de estudio

• Color natural:

# Conversión de coordenadas para poder realizar corte
shape_ja2 <- spTransform(shape_ja, CRS = crs(tif_nat14))

# Bandas a cortar
bandas_nat14 <- crop(tif_nat14, extent(shape_ja2))
corte_nat14 <- mask(bandas_nat14, shape_ja2)

# Gráfico
names(corte_nat14) <- paste("Banda", c(5, 3, 2), sep="")
plot(corte_nat14)

Se leyó el shape de área creado en ArcGis y se le introdujo las mismas coordenadas de la imagen para posteriormente recortarla y extraerla por máscara (procesos que se lleva a cabo en ArcGis con la función Extract by mask), cada banda se muestra por separado.

• Infrarojo:

# Bandas a cortar
bandas_inf14 <- crop(tif_inf14, extent(shape_ja2))
corte_inf14 <- mask(bandas_inf14, shape_ja2)

# Gráfico
names(corte_inf14) <- paste("Banda", c(5, 4, 3), sep="")
plot(corte_inf14)

Visualización Color Natural

plotRGB(tif_nat14, r=4, g=3, b=2, stretch="lin", main="RGB composición (b4,b3,b2) de Landsat-8")

plotRGB(corte_nat14, r=4, g=3, b=2, stretch="lin", main="RGB composición (b4,b3,b2) de Landsat-8")

Visualización Infrarojo

plotRGB(tif_inf14, r=5, g=4, b=3, stretch="lin", main="RGB composición (b8,b4,b2) de Landsat-8")

plotRGB(corte_inf14, r=5, g=4, b=3, stretch="lin", main="RGB composición (b8,b4,b2) de Landsat-8")

Se realizó la combinación de bandas para dar paso a las imágenes de infrarrojo y color natural, se procede a recortar el shape del área del humedal y se extrae por máscara.

NDVI

ndvi <- (tif_inf14[[1]] - tif_inf14[[2]])/(tif_inf14[[1]] + tif_inf14[[2]])
summary(ndvi)

##                layer
## Min.      -0.9926199
## 1st Qu.    0.1374676
## Median     0.2707559
## 3rd Qu.    0.5568026
## Max.       1.2758621
## NA's    1206.0000000

hist(ndvi)

Gráficos con NDVI

plot(ndvi, col = rev(terrain.colors(30)), main = 'NDVI desde landast8')

El intervalo de valores posibles oscila entre -1 y 1. Los valores negativos están asociados a zonas de agua y nieve. Valores positivos próximos a 0 representan zonas rocosas y desnudas que pueden adquirir algo de vegetación hasta llegar a valores próximos a 0.3. A partir de este valor encontramos presencia de vegetación densa. Cuanto mayor sea el valor más frondosa será la vegetación hasta adquirir valores próximos a 1. (Gisandbeers, 2016)[http://www.gisandbeers.com/calculo-del-indice-ndvi/] Se puede apreciar con claridad los cuerpos de agua que presentan un color rosado, el color verde denota vegetación y el amarillo área urbana. El histograma de los Niveles digitales muestra que hay un valor por encima de 1, por lo cúal se procede a la depuración de este valor

Quitando valores menores a -1 y mayores a 1 (valores erroneos).

veg1 <- calc(ndvi, function(x){x[x < -1] <- NA; return(x)})
veg1 <- calc(veg1, function(x){x[x > 1] <- NA; return(x)})
plot(veg1, main = 'Coberturas', col = rev(terrain.colors(8)))

Con la obteción de la gráfica de los niveles digitales, se pudo observar que había un valor mayor 1, se procedió a depurar los valores extremos para garantizar que los valor coincidan con los valores encontrados en la literaruta para el NDVI de -1 a 1

Quitando valores menores a 0.4 para mostrar vegetación.

veg2 <- calc(ndvi, function(x){x[x < 0.4] <- NA; return(x)})
veg2 <- calc(veg2, function(x){x[x > 1] <- NA; return(x)})
plot(veg2, main = 'Cobertura vegetal\nNDVI > 0.4')

Cuerpos de agua parecen estar entre -0.14 y -1.

cuts <- seq(minValue(veg1), maxValue(veg1), length.out=8)
cuts = round(cuts, digits = 2)
col.regions = brewer.pal(length(cuts)+2, "RdYlGn")
spplot(as(ndvi, 'SpatialGridDataFrame'),
       at=cuts, col.regions=col.regions,colorkey=list(labels=list(at=cuts),at=cuts),
       pretty=TRUE,scales=list(draw=T), main = "Cuerpos de agua\nNDVI entre -0.14 y -1")

Todo el proceso anterior se lleva a cabo con los años 2015,2016 y 2018

Los valores negativos indican la presencia de agua (color rosado y rojo), valores cercanos al 1 indican vegetación saludable (color verde y verde claro) y valores cercanos a cero indican áreas urbanizadas.

NDVI para el área de estudio:

Clasificación de imagen (algoritmo kmeans)

• Algoritmo k-means - idea intuitiva:

Agrupa las observaciones en K clusters distintos, donde el número K lo determina el analista antes de ejecutar del algoritmo. K-means clustering encuentra los K mejores clusters, entendiendo como mejor cluster aquel cuya varianza interna (intra-cluster variation) sea lo más pequeña posible. Se trata por lo tanto de un problema de optimización, en el que se reparten las observaciones en K clusters de forma que la suma de las varianzas internas de todos ellos sea lo menor posible. Amat Rodrigo, 2017

• Clasificación con k-means - 5 grupos:

• Se realizan 500 iteraciones.
• El algoritmo empieza con 3 grupos.

• El número de segmentos es igual a 5.

# Semilla para garantizar replicabilidad en el resultado
set.seed(1000)

# Obteniendo valores del ndvi para algoritmo kmeans
tablavalores <- getValues(ndvi)

# Algoritmo kmeans con 5 grupos
km <-kmeans(na.omit(tablavalores), centers = 5, iter.max = 500, nstart = 3)

rNA <- setValues(raster(ndvi), 0)#crear un raster en blanco con valores por defecto 0
for(i in 1:nlayers(ndvi)){ rNA[is.na(ndvi[[i]])] <- 1}# asignar valor de 1 a valores NA del raster
rNA <- getValues(rNA)#convertir rNA en un vector de valores
tablavalores <- as.data.frame(tablavalores)#convertir tablavalores en un data frame
tablavalores$clase[rNA==0] <- km$cluster #Si rNA es 0, asignar el valor del cluster k-mean 
tablavalores$clase[rNA==1] <- NA#Si rNA es 1, asignar NA
clases <- raster(ndvi)
clases <- setValues(clases, tablavalores$clase)#asignar los valores del raster clases a partir de  tablavalores$clase

colors<-c("black", "blue", "white","red","green")
plot(clases, legend=TRUE, col=colors)

En la imagen el azul parece indicar las carreteras, el color verde la vegetacíón y el color rojo las áreas urbanas. Aunque no es la mejor clasificación permite apreciar algunos segmentos.

• Clasificación con k-means - 6 grupos:

• Se realizan 500 iteraciones.
• El algoritmo empieza con 3 grupos.

• El número de segmentos es igual a 6.

# Semilla para garantizar replicabilidad en el resultado
set.seed(1000)

# Algoritmo kmeans con 5 grupos
km6 <-kmeans(na.omit(tablavalores), centers = 6, iter.max = 500, nstart = 3)

rNA <- setValues(raster(ndvi), 0)#crear un raster en blanco con valores por defecto 0
for(i in 1:nlayers(ndvi)){ rNA[is.na(ndvi[[i]])] <- 1}# asignar valor de 1 a valores NA del raster
rNA <- getValues(rNA)#convertir rNA en un vector de valores
tablavalores <- as.data.frame(tablavalores)#convertir tablavalores en un data frame
tablavalores$clase[rNA==0] <- km6$cluster #Si rNA es 0, asignar el valor del cluster k-mean 
tablavalores$clase[rNA==1] <- NA#Si rNA es 1, asignar NA
clases <- raster(ndvi)
clases <- setValues(clases, tablavalores$clase)#asignar los valores del raster clases a partir de  tablavalores$clase

colors<-c("orangered", "blue", "black","green","red", "magenta4")
plot(clases, legend=TRUE, col=colors)

Los seis grupos se pueden resumir para conocer sus valores (centroides) promedio:

km6$centers

##   tablavalores clase
## 1   0.39720287     1
## 2   0.54481958     2
## 3   0.23015482     3
## 4   0.79987596     5
## 5   0.64666709     2
## 6   0.09766934     4

• Clasificación con k-means - 7 grupos:

• Se realizan 500 iteraciones.
• El algoritmo empieza con 3 grupos.

• El número de segmentos es igual a 7.

# Semilla para garantizar replicabilidad en el resultado
set.seed(1000)

# Algoritmo kmeans con 5 grupos
km7 <-kmeans(na.omit(tablavalores), centers = 7, iter.max = 500, nstart = 3)

rNA <- setValues(raster(ndvi), 0)#crear un raster en blanco con valores por defecto 0
for(i in 1:nlayers(ndvi)){ rNA[is.na(ndvi[[i]])] <- 1}# asignar valor de 1 a valores NA del raster
rNA <- getValues(rNA)#convertir rNA en un vector de valores
tablavalores <- as.data.frame(tablavalores)#convertir tablavalores en un data frame
tablavalores$clase[rNA==0] <- km7$cluster #Si rNA es 0, asignar el valor del cluster k-mean 
tablavalores$clase[rNA==1] <- NA#Si rNA es 1, asignar NA
clases <- raster(ndvi)
clases <- setValues(clases, tablavalores$clase)#asignar los valores del raster clases a partir de  tablavalores$clase

colors<-c("orangered", "blue", "black","green","red", "magenta4", "yellow")
plot(clases, legend=TRUE, col=colors)

as.data.frame(km7$centers)

• Clasificación con k-means - 8 grupos:

• Se realizan 500 iteraciones.
• El algoritmo empieza con 3 grupos.

• El número de segmentos es igual a 8.

# Semilla para garantizar replicabilidad en el resultado
set.seed(1000)

# Algoritmo kmeans con 5 grupos
km8 <-kmeans(na.omit(tablavalores), centers = 8, iter.max = 500, nstart = 3)

rNA <- setValues(raster(ndvi), 0)#crear un raster en blanco con valores por defecto 0
for(i in 1:nlayers(ndvi)){ rNA[is.na(ndvi[[i]])] <- 1}# asignar valor de 1 a valores NA del raster
rNA <- getValues(rNA)#convertir rNA en un vector de valores
tablavalores <- as.data.frame(tablavalores)#convertir tablavalores en un data frame
tablavalores$clase[rNA==0] <- km8$cluster #Si rNA es 0, asignar el valor del cluster k-mean 
tablavalores$clase[rNA==1] <- NA#Si rNA es 1, asignar NA
clases <- raster(ndvi)
clases <- setValues(clases, tablavalores$clase)#asignar los valores del raster clases a partir de  tablavalores$clase

colors<-c("white", "blue", "black","forestgreen","red", "magenta4", "yellow", "cyan")
plot(clases, legend=TRUE, col=colors)

Resumen de los 8 grupos:

as.data.frame(km8$centers)

Clasificación de imagen (algoritmo CLARA)

• Algoritmo CLARA - idea intuitiva:

CLARA selecciona una muestra aleatoria de un tamaño determinado y le aplica el algoritmo de PAM (K-medoids) para encontrar los clusters óptimos acorde a esa muestra. Utilizando esos medoids se agrupan las observaciones de todo el set de datos. La calidad de los medoids resultantes se cuantifica con la suma total de las distancias entre cada observación del set de datos y su correspondiente medoid (suma total de distancias intra-clusters). AmatRodrigo,2017

En la clasificación del algoritmo CLARA con 6 grupos, aún no se puede apreciar bien la diferencia entre el cuerpo de agua y algunas áreas urbanas, por lo tanto una clasificación con 8 grupos como en el algoritmo K-means se podría apreciar mejor el cuerpo de agua.