Clasificación no supervisada
Clasificación no supervisada
En este capítulo exploramos la clasificación no supervisada. Existen varios algoritmos de clasificación no supervisada, y la elección del algoritmo puede afectar a los resultados. Exploraremos sólo un algoritmo (k-means) para ilustrar el principio general.
Para este ejemplo, seguiremos el esquema de clasificación de la Base de Datos Nacional de Cubierta Terrestre 2011 (NLCD 2011) para un subconjunto de las regiones del Valle Central. Utilizamos la imagen compuesta sin nubes de Landsat 5 con 6 bandas.
library(raster)## Loading required package: splandsat5 <- stack('data/rs/centralvalley-2011LT5.tif')
names(landsat5) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')//Pregunta 1: Haz una trama compuesta de 3 bandas de colores falsos de landsat5.
En la clasificación no supervisada, utilizamos los datos de reflectancia, pero no proporcionamos ningún dato de respuesta (es decir, no identificamos ningún píxel como perteneciente a una clase particular). Esto puede parecer extraño, pero puede ser útil cuando no tenemos mucho conocimiento previo de un área de estudio. O si se tiene un amplio conocimiento de la distribución de las clases de cobertura terrestre de interés, pero no hay datos terrestres específicos.
El algoritmo agrupa píxeles con características espectrales similares en grupos.
Aprenda más sobre los K-means y otros algoritmos no supervisados aquí.
Realizaremos una clasificación no supervisada en un subconjunto espacial de la capa ndvi. Aquí hay otra forma de calcular la NDVI. En este caso no usamos una función separada, sino que usamos una notación algebraica directa.
ndvi <- (landsat5[['NIR']] - landsat5[['red']]) / (landsat5[['NIR']] + landsat5[['red']])Haremos un agrupamiento de los medios de los datos de la NDVI. Primero usamos crop para hacer un subconjunto espacial de la ndvi, para permitir un procesamiento más rápido (puede seleccionar cualquier extensión usando la función drawExtent()).
# Extensión del cultivo de la capa ndvi
e <- extent(-121.807, -121.725, 38.004, 38.072)
# crop landsat by the extent
ndvi <- crop(ndvi, e)
ndvi## class : RasterLayer
## dimensions : 252, 304, 76608 (nrow, ncol, ncell)
## resolution : 0.0002694946, 0.0002694946 (x, y)
## extent : -121.807, -121.725, 38.00413, 38.07204 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : -0.3360085, 0.7756007 (min, max)nr <- getValues(ndvi)
str(nr)## num [1:76608] 0.245 0.236 0.272 0.277 0.277 ...Tengan en cuenta que getValues convirtió la ndvi RasterLayer en una matriz. Ahora realizaremos la agrupación de los medios en la matriz e inspeccionaremos la salida.
set.seed(99)#Semilla
# Queremos crear 10 cúmulos, permitir 500 iteraciones, comenzar con 5 conjuntos aleatorios usando el método de "Lloyd"
kmncluster <- kmeans(na.omit(nr), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")
# kmeans devuelve un objeto de clase "kmeans"
str(kmncluster)## List of 9
## $ cluster : int [1:76608] 4 4 3 3 3 3 3 4 4 4 ...
## $ centers : num [1:10, 1] 0.55425 0.00498 0.29997 0.20892 -0.20902 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:10] "1" "2" "3" "4" ...
## .. ..$ : NULL
## $ totss : num 6459
## $ withinss : num [1:10] 5.69 6.13 4.91 4.9 5.75 ...
## $ tot.withinss: num 55.8
## $ betweenss : num 6403
## $ size : int [1:10] 8932 4550 7156 6807 11672 8624 8736 5040 9893 5198
## $ iter : int 108
## $ ifault : NULL
## - attr(*, "class")= chr "kmeans"kmeans devuelve un objeto con 9 elementos. La longitud del elemento del cluster dentro del kmncluster es 76608 que es igual a la longitud de nr creado a partir del ndvi. Los valores de las celdas del kmncluster\(cluster van de 1 a 10 correspondientes al numero de entrada del cluster que proveimos en la funcion de kmeans. kmncluster\)cluster indica la etiqueta del cluster para el pixel correspondiente. Necesitamos convertir los valores de kmncluster$cluster de nuevo a RasterLayer de la misma dimensión que la ndvi.
# Usar el objeto ndvi para fijar los valores del cúmulo a un nuevo raster
knr <- setValues(ndvi, kmncluster$cluster)
# También puedes hacerlo así
knr <- raster(ndvi)
values(knr) <- kmncluster$cluster
knr## class : RasterLayer
## dimensions : 252, 304, 76608 (nrow, ncol, ncell)
## resolution : 0.0002694946, 0.0002694946 (x, y)
## extent : -121.807, -121.725, 38.00413, 38.07204 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 1, 10 (min, max)Podemos ver que knr es un RasterLayer pero no sabemos qué cúmulo (1-10) pertenece a qué clase de cubierta terrestre (y si pertenece a una clase que reconoceríamos). Puedes averiguarlo trazándolos lado a lado con una capa de referencia y usando un color único para cada cúmulo.
# Definir un vector de color para 10 cúmulos (aprender más sobre la configuración del color más tarde)
mycolor <- c("#fef65b","#ff0000", "#daa520","#0000ff","#0000ff","#00ff00","#cbbeb5",
"#c3ff5b", "#ff7373", "#00ff00", "#808080")
par(mfrow = c(1,2))
plot(ndvi, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')
plot(knr, main = 'Unsupervised classification', col = mycolor )
Mientras que para otros fines suele ser mejor definir más clases (y posiblemente fusionarlas más tarde), una clasificación simple como ésta podría ser útil, por ejemplo, fusionar los grupos 4 y 5 para construir una máscara de agua para el año 2011.
Puedes cambiar los colores en mi micrón. Aprende más sobre la selección de colores en R aquí y aquí.
Pregunta 2: Traza la banda 3 RGB de landsat5 para el subconjunto (extensión e) y el resultado de los medios de los clusters uno al lado del otro y haz una tabla de etiquetas de cobertura de tierra para los clusters. Por ejemplo, los grupos 4 y 5 son agua.