Introducción

En este ejercicio se realizará la clasificación no supervisada de un mosaico de imágenes multiespectrales de alta resolución capturadas con un vehículo aéreo no tripulado (UAV) sobre un cultivo experimental de papa ubicado en el municipio de Subachoque, Cundinamarca.

Datos

Se utilizará un mosaico de imágenes áereas multiespectrales obtenidas con el sensor Micasense RedEdge, la información de las bandas del sensor se muestra a continuación.

rededge <- read.csv("./datos/rededge.csv", sep = ",")

knitr::kable(rededge)

X..Banda	Nombre	Centro.de.Banda..nm.	Ancho.Banda.FWHM..nm.
Banda 1	Azul	475	20
Banda 2	Verde	560	20
Banda 3	Rojo	668	10
Banda 4	Infrarrojo cercano	840	40
Banda 5	Red edge	717	10
## Método
El ejercici	o consta de las sigui	entes etapas:

Realizar la clasificación no supervisada mediante los algoritmos K-Means y CLARA.
Comparar el desempeño de los dos algoritmos de clasificación.
Evaluar el número de grupos que mejor resultados produce en la clasificación.

A continuación cargamos el mosaico de la zona de interés:

library(raster)

## Loading required package: sp

# Load the multi-band GeoTIFF file with brick function
rst <- brick("./datos/M_20180526_SPEC_3116_clip_10cm.tif")
rst[rst==0] <- NA
# Change band names      
names(rst) <- paste("b",1:5,sep="")

Composición falso color del mosaico.

plotRGB(rst, r=4, g=3, b=2, scale=10000, stretch="lin", main="RGB composite (b4,b3,b2) Micasense RedEdge")

Ahora empieza el proceso de clasificación no supervisada.

library(cluster)

# Extract all values from the raster into a data frame
rstDF <- values(rst)

# Check NA's in the data
idx <- complete.cases(rstDF)

# Initiate the raster datasets that will hold all clustering solutions 
# from 2 groups/clusters up to 12
rstKM <- raster(rst[[1]])
rstCLARA <- raster(rst[[1]])


for(nClust in 2:12){
  
  cat("-> Clustering data for nClust =",nClust,"......")
  
  # Perform K-means clustering
  km <- kmeans(rstDF[idx,], centers = nClust, iter.max = 50)
  
  # Perform CLARA's clustering (using manhattan distance)
  cla <- clara(rstDF[idx, ], k = nClust, metric = "manhattan")
  
  # Create a temporary integer vector for holding cluster numbers
  kmClust <- vector(mode = "integer", length = ncell(rst))
  claClust <- vector(mode = "integer", length = ncell(rst))
  
  # Generate the temporary clustering vector for K-means (keeps track of NA's)
  kmClust[!idx] <- NA
  kmClust[idx] <- km$cluster
  
  # Generate the temporary clustering vector for CLARA (keeps track of NA's too ;-)
  claClust[!idx] <- NA
  claClust[idx] <- cla$clustering
  
  # Create a temporary raster for holding the new clustering solution
  # K-means
  tmpRstKM <- raster(rst[[1]])
  # CLARA
  tmpRstCLARA <- raster(rst[[1]])

  # Set raster values with the cluster vector
  # K-means
  values(tmpRstKM) <- kmClust
  # CLARA
  values(tmpRstCLARA) <- claClust
  
  # Stack the temporary rasters onto the final ones
  if(nClust==2){
    rstKM    <- tmpRstKM
    rstCLARA <- tmpRstCLARA
  }else{
    rstKM    <- stack(rstKM, tmpRstKM)
    rstCLARA <- stack(rstCLARA, tmpRstCLARA)
  }
  
  cat(" done!\n\n")
}

## -> Clustering data for nClust = 2 ...... done!
## 
## -> Clustering data for nClust = 3 ...... done!
## 
## -> Clustering data for nClust = 4 ...... done!
## 
## -> Clustering data for nClust = 5 ...... done!
## 
## -> Clustering data for nClust = 6 ......

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 116837450)

##  done!
## 
## -> Clustering data for nClust = 7 ......

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 116837450)

##  done!
## 
## -> Clustering data for nClust = 8 ...... done!
## 
## -> Clustering data for nClust = 9 ...... done!
## 
## -> Clustering data for nClust = 10 ...... done!
## 
## -> Clustering data for nClust = 11 ......

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 116837450)

##  done!
## 
## -> Clustering data for nClust = 12 ......

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 116837450)

##  done!

# Write the clustering solutions for each algorithm
writeRaster(rstKM,"./datos/M_20180526_SPEC_3116_clip_10cm_KMeans_nc2_12-1.tif", overwrite=TRUE)
writeRaster(rstCLARA,"./datos/M_20180526_SPEC_3116_clip_10cminst_CLARA_nc2_12-1.tif", overwrite=TRUE)

Evaluación del desempeño de los algoritmos

library(clusterCrit)


# Start a data frame that will store all silhouette values
# for k-means and CLARA   
clustPerfSI <- data.frame(nClust = 2:12, SI_KM = NA, SI_CLARA = NA)


for(i in 1:nlayers(rstKM)){ # Iterate through each layer
  
  cat("-> Evaluating clustering performance for nClust =",(2:12)[i],"......")
  
  # Extract random cell samples stratified by cluster
  cellIdx_RstKM <- sampleStratified(rstKM[[i]], size = 2000)
  cellIdx_rstCLARA <- sampleStratified(rstCLARA[[i]], size = 2000)
  
  # Get cell values from the Stratified Random Sample from the raster 
  # data frame object (rstDF)
  rstDFStRS_KM <- rstDF[cellIdx_RstKM[,1], ]
  rstDFStRS_CLARA <- rstDF[cellIdx_rstCLARA[,1], ]
  
  # Make sure all columns are numeric (intCriteria function is picky on this)
  rstDFStRS_KM[] <- sapply(rstDFStRS_KM, as.numeric)
  rstDFStRS_CLARA[] <- sapply(rstDFStRS_CLARA, as.numeric)
  
  # Compute the sample-based Silhouette index for: 
  #    
  # K-means
  clCritKM <- intCriteria(traj = rstDFStRS_KM, 
                          part = as.integer(cellIdx_RstKM[,2]), 
                          crit = "Silhouette")
  # and CLARA
  clCritCLARA <- intCriteria(traj = rstDFStRS_CLARA, 
                             part = as.integer(cellIdx_rstCLARA[,2]), 
                             crit = "Silhouette")

  # Write the silhouette index value to clustPerfSI data frame holding 
  # all results
  clustPerfSI[i, "SI_KM"]    <- clCritKM[[1]][1]
  clustPerfSI[i, "SI_CLARA"] <- clCritCLARA[[1]][1]
  
  cat(" done!\n\n")
  
}

## -> Evaluating clustering performance for nClust = 2 ...... done!
## 
## -> Evaluating clustering performance for nClust = 3 ...... done!
## 
## -> Evaluating clustering performance for nClust = 4 ...... done!
## 
## -> Evaluating clustering performance for nClust = 5 ...... done!
## 
## -> Evaluating clustering performance for nClust = 6 ...... done!
## 
## -> Evaluating clustering performance for nClust = 7 ...... done!
## 
## -> Evaluating clustering performance for nClust = 8 ...... done!
## 
## -> Evaluating clustering performance for nClust = 9 ...... done!
## 
## -> Evaluating clustering performance for nClust = 10 ...... done!
## 
## -> Evaluating clustering performance for nClust = 11 ...... done!
## 
## -> Evaluating clustering performance for nClust = 12 ...... done!

write.csv(clustPerfSI, file = "./datos/clustPerfSI.csv", row.names = FALSE)

Para la evaluación se tiene en cuenta un idice de silueta que permite comparar los resultados según cada cluster.

knitr::kable(clustPerfSI, digits = 3, align = "c", 
             col.names = c("#clusters","Avg. Silhouette (k-means)","Avg. Silhouette (CLARA)"))

#clusters	Avg. Silhouette (k-means)	Avg. Silhouette (CLARA)
2	0.412	0.407
3	0.365	0.341
4	0.302	0.276
5	0.304	0.262
6	0.319	0.321
7	0.302	0.261
8	0.290	0.255
9	0.277	0.256
10	0.275	0.246
11	0.262	0.252
12	0.234	0.225

Ahora se comparan los dos algoritmos.

plot(clustPerfSI[,1], clustPerfSI[,2], 
     xlim = c(1,13), ylim = range(clustPerfSI[,2:3]), type = "n", 
     ylab="Avg. Silhouette Index", xlab="# of clusters",
     main="Silhouette index by # of clusters")

# Plot Avg Silhouette values across # of clusters for K-means
lines(clustPerfSI[,1], clustPerfSI[,2], col="red")
# Plot Avg Silhouette values across # of clusters for CLARA
lines(clustPerfSI[,1], clustPerfSI[,3], col="blue")

# Grid lines
abline(v = 1:13, lty=2, col="light grey")
abline(h = seq(0.30,0.44,0.02), lty=2, col="light grey")

legend("topright", legend=c("K-means","CLARA"), col=c("red","blue"), lty=1, lwd=1)

Finalmente, podeos ver los resultados de los algoritmos de clasificación.

rstKM_df <-  as.data.frame(rstKM[[4]], xy = TRUE)
rstCLARA_df <-  as.data.frame(rstCLARA[[4]], xy = TRUE)

#unique(rstkm_df$layer.2.2.1)
#factor(rstKM_df$layer.2.2.1)

library(ggplot2)
my_col <- c("#00A600FF", "#ECB176FF", "#00FFAAFF","#FF00AAFF","#AAFF00FF")
#my_col2 <- c("#AAFF00FF",  "#00FFAAFF","#FF00AAFF","#00A600FF","#ECB176FF")

ggplot() +
 geom_raster(data = rstKM_df, aes(x = x, y = y,
                                      fill = factor(layer.2.2.1))) + 
    scale_fill_manual(values = my_col) +
    ggtitle("K-MEANS") +
    coord_quickmap()

## Warning: Removed 3054341 rows containing missing values (geom_raster).

ggplot() +
 geom_raster(data = rstCLARA_df, aes(x = x, y = y,
                                      fill = factor(layer.2.2.1))) + 
    scale_fill_manual(values = my_col) + 
    ggtitle("CLARA") +
    coord_quickmap()

## Warning: Removed 3054341 rows containing missing values (geom_raster).

Discusión

Las diferencias en el desempeño de los dos algoritmos son pequeñas como se evidencia en la tabla de ídices de silueta y el grafico de comparación. Los resultados en la clasificación son visualmente similares, en los dos resultados, por ejemplo, el agua es clasificada como si se tratara del suelo del cultivo.

En la zona sur occidental del mosaico hay un área afectada con tizón tardío que en ambos algoritmos es clasificada dentro de la misma clase que el suelo desnudo y el agua.

Estos resultados poco óptimos pueden tener una razón, que es la calibración radiométrica de las imágenes aéreas obtenidas, en este caso el mosaico empleado proviene de imágenes aéreas sin calibración radiométrica, se esperaría que una buena calibración radiométrica pudiera generar resultados mucho mejores.

1 - Clasificación No supervisada

Jorge Luis Rodriguez based on Joao Goncalves

06-07-2018

Introducción

Datos

Evaluación del desempeño de los algoritmos

Discusión