Lectura de imágenes Infrarojo (3 bandas)

# Bibliotecas
library(raster)
library(sp)
library(tidyverse)

# Importando imágenes 
inf14 <- brick("imagenes/satelitales/Nuevas_completas/infrarrojo2014.tif")
inf15 <- brick("imagenes/satelitales/Nuevas_completas/infrarrojo2015.tif")
inf16 <- brick("imagenes/satelitales/Nuevas_completas/inflarrojo2016.tif")
inf18 <- brick("imagenes/satelitales/Nuevas_completas/infrarrojo2018.tif")

Ejemplo: número de bandas de inf14

nlayers(inf14)

## [1] 3

Ejemplo: información de inf14

inf14

## class      : RasterBrick 
## dimensions : 262, 473, 123926, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 593655, 607845, 517695, 525555  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /home/edimer/Documentos/U/Otros/Mariana/Avances/Final/imagenes/satelitales/Nuevas_completas/infrarrojo2014.tif 
## names      : infrarrojo2014.1, infrarrojo2014.2, infrarrojo2014.3 
## min values :                2,             -253,             -232 
## max values :             6472,             6848,             6896

Ejemplo: coordenadas de inf14

crs(inf14)

## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0

Ejemplo: dimensiones de inf14

dim(inf14)

## [1] 262 473   3

Visualización banda 5

par(mfrow = c(2, 2))
plot(inf14[[1]], main = "Banda 5 - 2014")
plot(inf15[[1]], main = "Banda 5 - 2015")
plot(inf16[[1]], main = "Banda 5 - 2016")
plot(inf18[[1]], main = "Banda 5 - 2018")

Visualización 2014

# Cambiando nombres a las capas
names(inf14) <- c("Banda5", "Banda 4", "Banda 3")

# Asignando colores
library(RColorBrewer)
colores <- brewer.pal(n = 7, name = "RdYlGn")

# Gráfico
spplot(inf14, main = "2014", col.regions = colores, cuts = 7)

Visualización 2015

# Cambiando nombres a las capas
names(inf15) <- c("Banda5", "Banda 4", "Banda 3")

# Gráfico
spplot(inf15, main = "2015", col.regions = colores, cuts = 7)

Visualización 2016

# Cambiando nombres a las capas
names(inf16) <- c("Banda5", "Banda 4", "Banda 3")

# Gráfico
spplot(inf16, main = "2016", col.regions = colores, cuts = 7)

Visualización 2018

# Cambiando nombres a las capas
names(inf18) <- c("Banda5", "Banda 4", "Banda 3")

# Gráfico
spplot(inf18, main = "2018", col.regions = colores, cuts = 7)

NDVI

Cálculo de NDVI

\[NDVI = \frac{Infrarojo\ cercano\ -\ Rojo}{Infrarojo\ cercano\ +\ Rojo} = \frac{Banda5\ -\ Banda4}{Banda5\ +\ Banda4}\]

# Calculando NDVI
ndvi14 <- (inf14[[1]] - inf14[[2]])/(inf14[[1]] + inf14[[2]])
ndvi15 <- (inf15[[1]] - inf15[[2]])/(inf15[[1]] + inf15[[2]])
ndvi16 <- (inf16[[1]] - inf16[[2]])/(inf16[[1]] + inf16[[2]])
ndvi18 <- (inf18[[1]] - inf18[[2]])/(inf18[[1]] + inf18[[2]])

Distribución de NDVI

par(mfrow = c(2, 2))
hist(ndvi14, main = "NDVI 2014")
hist(ndvi15, main = "NDVI 2015")
hist(ndvi16, main = "NDVI 2016")
hist(ndvi18, main = "NDVI 2018")

Eliminando valores menores a -1 y mayores a 1

# Eliminando valores > 1 y < -1
ndvi14 <- calc(ndvi14, function(x){x[x < -1] <- NA; return(x)})
ndvi14 <- calc(ndvi14, function(x){x[x > 1] <- NA; return(x)})
ndvi15 <- calc(ndvi15, function(x){x[x < -1] <- NA; return(x)})
ndvi15 <- calc(ndvi15, function(x){x[x > 1] <- NA; return(x)})
ndvi16 <- calc(ndvi16, function(x){x[x < -1] <- NA; return(x)})
ndvi16 <- calc(ndvi16, function(x){x[x > 1] <- NA; return(x)})
ndvi18 <- calc(ndvi18, function(x){x[x < -1] <- NA; return(x)})
ndvi18 <- calc(ndvi18, function(x){x[x > 1] <- NA; return(x)})

par(mfrow = c(2, 2))
plot(ndvi14, main = "NDVI 2014", col = colores, cuts = 6)
plot(ndvi15, main = "NDVI 2015", col = colores, cuts = 6)
plot(ndvi16, main = "NDVI 2016", col = colores, cuts = 6)
plot(ndvi18, main = "NDVI 2018", col = colores, cuts = 6)

Graficando valores NDVI mayores a 0.4 (vegetación)

ndvi14.2 <- calc(ndvi14, function(x){x[x < 0.4] <- NA; return(x)})
ndvi15.2 <- calc(ndvi15, function(x){x[x < 0.4] <- NA; return(x)})
ndvi16.2 <- calc(ndvi16, function(x){x[x < 0.4] <- NA; return(x)})
ndvi18.2 <- calc(ndvi18, function(x){x[x < 0.4] <- NA; return(x)})

par(mfrow = c(2, 2))
plot(ndvi14.2, main = "NDVI >0.4 2014", col = colores, cuts = 6)
plot(ndvi15.2, main = "NDVI >0.4 2015", col = colores, cuts = 6)
plot(ndvi16.2, main = "NDVI >0.4 2016", col = colores, cuts = 6)
plot(ndvi18.2, main = "NDVI >0.4 2018", col = colores, cuts = 6)

Graficando valores NDVI mayores a 0.4 (vegetación)

par(mfrow = c(2, 2))
plot(ndvi14.2, main = "NDVI >0.4 2014")
plot(ndvi15.2, main = "NDVI >0.4 2015")
plot(ndvi16.2, main = "NDVI >0.4 2016")
plot(ndvi18.2, main = "NDVI >0.4 2018")

NDVI en área de estudio

Corte del área de interés

# Importando shape
shape_ja <- shapefile("imagenes/Areas/Humedal_2014_01_01.shp")

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

# Bandas a cortar
ndvi14_z <- crop(ndvi14, extent(shape_ja2))
ndvi14_z <- mask(ndvi14_z, shape_ja2)

ndvi15_z <- crop(ndvi15, extent(shape_ja2))
ndvi15_z <- mask(ndvi15_z, shape_ja2)

ndvi16_z <- crop(ndvi16, extent(shape_ja2))
ndvi16_z <- mask(ndvi16_z, shape_ja2)

ndvi18_z <- crop(ndvi18, extent(shape_ja2))
ndvi18_z <- mask(ndvi18_z, shape_ja2)

NDVI Área de estudio (1)

# Gráfico
par(mfrow = c(2, 2))
plot(ndvi14_z , main = "NDVI 2014", col = colores, cuts = 6)
plot(ndvi15_z , main = "NDVI 2015", col = colores, cuts = 6)
plot(ndvi16_z , main = "NDVI 2016", col = colores, cuts = 6)
plot(ndvi18_z , main = "NDVI 2018", col = colores, cuts = 6)

NDVI Área de estudio (2)

# Gráfico
par(mfrow = c(2, 2))
plot(ndvi14_z , main = "NDVI 2014", col = c("dodgerblue", "firebrick", "forestgreen"))
plot(ndvi15_z , main = "NDVI 2015", col = c("dodgerblue", "firebrick", "forestgreen"))
plot(ndvi16_z , main = "NDVI 2016", col = c("dodgerblue", "firebrick", "forestgreen"))
plot(ndvi18_z , main = "NDVI 2018", col = c("dodgerblue", "firebrick", "forestgreen"))

Espejo de agua

# Calculando valores inferiores a -0.2
# 2014
df_14 <- as.data.frame(ndvi14_z)
df_14 <- na.omit(df_14)
prop.table(table(df_14$layer <= 0))
## 
##      FALSE       TRUE 
## 0.96250455 0.03749545

# 2015
df_15 <- as.data.frame(ndvi15_z)
df_15 <- na.omit(df_15)
prop.table(table(df_15$layer <= 0))
## 
##      FALSE       TRUE 
## 0.96978522 0.03021478

# 2016
df_16 <- as.data.frame(ndvi16_z)
df_16 <- na.omit(df_16)
prop.table(table(df_16$layer <= 0))
## 
##      FALSE       TRUE 
## 0.96869312 0.03130688

# 2018
df_18 <- as.data.frame(ndvi18_z)
df_18 <- na.omit(df_18)
prop.table(table(df_18$layer <= 0))
## 
##       FALSE        TRUE 
## 0.999271933 0.000728067

Resultado: el 3.74% (0.03749545) del área de estudio podría ser considerada como espejo de agua en el año 2014. Estos valores se validan posteriormente a través de clasificación no supervisada.

Cobertura vegetal (NDVI >0.4)

# Calculando valores inferiores a -0.2
# 2014
prop.table(table(df_14$layer >= 0.4))
## 
##      FALSE       TRUE 
## 0.09683291 0.90316709

# 2015
prop.table(table(df_15$layer >= 0.4))
## 
##     FALSE      TRUE 
## 0.1066618 0.8933382

# 2016
prop.table(table(df_16$layer >= 0.4))
## 
##     FALSE      TRUE 
## 0.1048416 0.8951584

# 2018
prop.table(table(df_18$layer >= 0.4))
## 
##     FALSE      TRUE 
## 0.1248635 0.8751365

Resultado: el 90.31% (0.90316709) del área de estudio podría ser considerada como cobertura vegetal (NDVI mayor a 0.4) en el año 2014. El área cubierta por vegetación parece variar a través del tiempo.

Clasificación no supervisada NDVI 2014

Algoritmos

Algoritmo k-means de 2 a 6 grupos con un máximo de 50 iteraciones e iniciando siempre con 2 grupos.
Algoritmo CLARA (Clustering Large Applications) de 2 a 6 grupos con tamaños de muestra fijos de 2500 observaciones (pixeles) y 5 muestras.
Comparación a través del índice silhouette con valores entre -1 y 1. Valores cercanos a 1 indican que la observación fue correctamente asignada al clúster. Dado que el método es costoso computacionalmente, se obtienen muestras aleatorias estratificadas de tamaño 5000, para calcular el índice silhouette.
Clasificación basada en la banda 1 (B5) de Landsat-8.
Las imágenes con los grupos (clúster) se exporta en formato .tif
Los códigos fueron adaptados de los siguientes documentos:

# Bibliotecas para clustering
library(cluster)
library(clusterCrit)

Clusterización con k-means y CLARA

# Extracción de valores del raster
bandasDF14 <- values(ndvi14)

# Recordando índices de celdas con valores
idx <- complete.cases(bandasDF14)

# Asignando bandas para clusterización
bandasKM <- raster(ndvi14[[1]])
bandasCLARA <- raster(ndvi14[[1]])

for(cluster in 2:6){
  cat("-> Proceso de clusterización para k=",cluster,"......")
  
  # Realización del agrupamiento o cluser de  K-means
  km <- kmeans(bandasDF14[idx], centers = cluster, iter.max = 50, nstart = 2)
  
  # Realización del agrupamiento o cluser de CLARA utilizando distancia manhattan
  cla <- clara(bandasDF14[idx], k = cluster, metric = "manhattan", samples = 5,
               sampsize = 2500)
  
  # Crear un vector que contenga los numeros de los clusters
  kmClust <- vector(mode = "integer", length = ncell(ndvi14))
  claClust <- vector(mode = "integer", length = ncell(ndvi14))
  
  # Genera el vector de los clusters de Kmeans, el cual realiza un seguimiento de los NA
  kmClust[!idx] <- NA
  kmClust[idx] <- km$cluster
  
  # Genera el vector temporal de los clusters de Clara y realiza seguimiento de los NA
  claClust[!idx] <- NA
  claClust[idx] <- cla$clustering
  
  # Crea el raster temporal que contiene los nuevos agrupamientos o clusters solucion
  # Para K-means
  tmpbandasKM <- raster(ndvi14[[1]])
  # Para CLARA
  tmpbandasCLARA <- raster(ndvi14[[1]])
  
  # Establece los  valores del raster con el vector cluster 
  # Para  K-means
  values(tmpbandasKM) <- kmClust
  # Para CLARA
  values(tmpbandasCLARA) <- claClust
  
  # Stack del raster temporal en los finales
  if(cluster == 2){
    bandasKM    <- tmpbandasKM
    bandasCLARA <- tmpbandasCLARA
  }else{
    bandasKM    <- stack(bandasKM, tmpbandasKM)
    bandasCLARA <- stack(bandasCLARA, tmpbandasCLARA)
  }
  
  cat(" done!\n\n")
}

Exportando imágenes

writeRaster(bandasKM,"Kmeans_26.tif", overwrite=TRUE)
writeRaster(bandasCLARA,"CLARA_26.tif", overwrite=TRUE)

Importando imágenes

# Leyendo imágenes clasificadas
kmeans_26 <- brick("Kmeans_26.tif")
clara_26  <- brick("CLARA_26.tif")

Evaluando la clasificiación

# Extracción de valores del raster
bandasDF14 <- values(ndvi14)

# Dataframe para almacenar índices Silhouette
clustPerfSI <- data.frame(nClust = 2:6, SI_KM =NA, SI_CLARA=NA)

# Ciclo para calcular Silhouette a traves de cada capa
for(i in 1:nlayers(kmeans_26)){
  
  cat("-> Evaluación de la clasificación con k=",(2:6)[i],"......")
  
  # Extraccion de muestras aleatorias estratificadas por cluster
  cellIdx_bandasKM <- sampleStratified(kmeans_26[[i]], size = 5000)
  cellIdx_bandasCLARA <- sampleStratified(clara_26[[i]], size = 5000)
  
  # Obtiene los valores de celda de la muestra aleatoria estratificada del ráster
  bandasDFStRS_KM <- bandasDF14[cellIdx_bandasKM[,1]]
  bandasDFStRS_CLARA <- bandasDF14[cellIdx_bandasCLARA[,1]]
  
  # Para asegúrese que todas las columnas sean numéricas ya que la función intCriteria lo exige.
  bandasDFStRS_KM[] <- sapply(bandasDFStRS_KM, as.numeric)
  bandasDFStRS_CLARA[] <- sapply(bandasDFStRS_CLARA, as.numeric)
  
  # Calculo del índice de siluetas basado en muestras para: 
  #    
  # K-means
  clCritKM <- intCriteria(traj = as.matrix(bandasDFStRS_KM),
                          part = as.integer(cellIdx_bandasKM[,2]), 
                          crit = "Silhouette")
  # CLARA
  clCritCLARA <- intCriteria(traj = as.matrix(bandasDFStRS_CLARA), 
                             part = as.integer(cellIdx_bandasCLARA[,2]), 
                             crit = "Silhouette")
  
  # Escribir los valores del índice de silueta en el data frame clustPerfSI.
  # Resultados
  clustPerfSI[i, "SI_KM"]    <- clCritKM[[1]][1]
  clustPerfSI[i, "SI_CLARA"] <- clCritCLARA[[1]][1]
  
  cat(" done!\n\n")
  
}

## -> Evaluación de la clasificación con k= 2 ...... done!
## 
## -> Evaluación de la clasificación con k= 3 ...... done!
## 
## -> Evaluación de la clasificación con k= 4 ...... done!
## 
## -> Evaluación de la clasificación con k= 5 ...... done!
## 
## -> Evaluación de la clasificación con k= 6 ...... done!

Gráfico de índice de silhouette:

clustPerfSI %>% 
  gather(key = "algoritmo", value = "sill", -nClust) %>%
  mutate(algoritmo = factor(algoritmo, labels = c("CLARA", "Kmeans"))) %>% 
  ggplot(data = ., aes(x = factor(nClust), y = sill, color = algoritmo)) +
  geom_point() +
  geom_line(aes(group = algoritmo)) +
  labs(x = "Número de clúster", y = "Silhouette",
       subtitle = "Validación de algoritmos de clusterización",
       color = "") +
  theme_classic() +
  theme(legend.position = "bottom")

K-means con k = 2

plot(kmeans_26[[1]], col = c("forestgreen", "orangered"),
     main = "K-means\nk = 2")

Proporción (área) por clúster:

# Valores <= 1.5 (clúster 1) --> aproximado
round(prop.table(table(values(kmeans_26[[1]]) <= 1.5))*100, digits = 2)
## 
## FALSE  TRUE 
## 66.27 33.73

K-means con k = 3

plot(kmeans_26[[2]], col = c("orangered", "forestgreen", "black"),
     main = "K-means\nk = 3")

Proporción (área) por clúster:

# Valores <= 1.5 (clúster 1) --> aproximado
round(prop.table(table(values(kmeans_26[[2]]) <= 1.5))*100, digits = 2)
## 
## FALSE  TRUE 
## 74.23 25.77

# Valores > 1.5 y menores a 2.5 (clúster 2) --> aproximado
round(prop.table(table(values(kmeans_26[[2]]) > 1.5
                       & values(kmeans_26[[2]]) <= 2.5))*100, digits = 2)
## 
## FALSE  TRUE 
## 49.46 50.54

# Valores mayores a 2.5 (clúster 3) --> aproximado
round(prop.table(table(values(kmeans_26[[2]]) > 2.5))*100, digits = 2)
## 
## FALSE  TRUE 
## 76.31 23.69

CLARA con k = 2

plot(clara_26[[1]], col = c("forestgreen", "orangered"),
     main = "CLARA\nk = 2")

Proporción (área) por clúster:

# Valores <= 1.5 (clúster 1) --> aproximado
round(prop.table(table(values(clara_26[[1]]) <= 1.5))*100, digits = 2)
## 
## FALSE  TRUE 
## 65.82 34.18

CLARA con k = 3

plot(clara_26[[2]], col = c("orangered", "forestgreen", "black"),
     main = "CLARA\nk = 3")

Proporción (área) por clúster:

# Valores <= 1.5 (clúster 1) --> aproximado
round(prop.table(table(values(clara_26[[2]]) <= 1.5))*100, digits = 2)
## 
## FALSE  TRUE 
## 75.44 24.56

# Valores > 1.5 y menores a 2.5 (clúster 2) --> aproximado
round(prop.table(table(values(clara_26[[2]]) > 1.5
                       & values(clara_26[[2]]) <= 2.5))*100, digits = 2)
## 
## FALSE  TRUE 
## 73.26 26.74

# Valores mayores a 2.5 (clúster 3) --> aproximado
round(prop.table(table(values(clara_26[[2]]) > 2.5))*100, digits = 2)
## 
## FALSE  TRUE 
##  51.3  48.7

K-means con k = 4

plot(kmeans_26[[3]], col = c("blue", "forestgreen", "orangered", "black"),
     main = "K-means\nk = 4")

CLARA con k = 4

plot(clara_26[[3]], col = c("black", "forestgreen", "blue", "orangered"),
     main = "CLARA\nk = 4")

K-means con k = 5

plot(kmeans_26[[4]], col = c("blue", "black", "orangered", "forestgreen",
                             "yellow2"), main = "K-means\nk = 5")

CLARA con k = 5

plot(clara_26[[4]], col = c("black", "forestgreen", "blue", "yellow2",
                             "orangered"), main = "CLARA\nk = 5")

K-means con k = 6

plot(kmeans_26[[5]], col = c("forestgreen", "black", "red", "yellow2",
                             "navy", "magenta3"), main = "K-means\nk = 6")

CLARA con k = 6

plot(clara_26[[5]], col = c("yellow2", "black", "forestgreen", "red",
                             "magenta3", "navy"), main = "CLARA\nk = 6")

Clasificación de imágenes satelitales

Aprendizaje no supervisado

Mariana Flórez Restrepo, Edimer David Jaramillo

2019-08-15

Lectura de imágenes Infrarojo (3 bandas)

Visualización banda 5

Visualización 2014

Visualización 2015

Visualización 2016

Visualización 2018

NDVI

Cálculo de NDVI

Distribución de NDVI

Eliminando valores menores a -1 y mayores a 1

Graficando valores NDVI mayores a 0.4 (vegetación)

Graficando valores NDVI mayores a 0.4 (vegetación)

NDVI en área de estudio

Corte del área de interés

NDVI Área de estudio (1)

NDVI Área de estudio (2)

Espejo de agua

Cobertura vegetal (NDVI >0.4)

Clasificación no supervisada NDVI 2014

Algoritmos

Clusterización con k-means y CLARA

Exportando imágenes

Importando imágenes

Evaluando la clasificiación

K-means con k = 2

K-means con k = 3

CLARA con k = 2

CLARA con k = 3

K-means con k = 4

CLARA con k = 4

K-means con k = 5

CLARA con k = 5

K-means con k = 6

CLARA con k = 6