Imagen LandSat 8, Sibate 2015 - Clasificacion no supervisada

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

En este capitulo exploramos la clasificacion no supervisada. Existen varios algoritmos de clasificación no supervisados, y la eleccion del algoritmo puede afectar los resultados. Exploraremos un solo algoritmo (k-means) para ilustrar el principio general

Para este ejemplo, seguiremos el esquema de clasificacion de la -Base de datos nacional de cobertura de la tierra 2011 (NLCD 2011)- para un subconjunto de las regiones del Valle Central. Utilizamos imagenes compuestas sin nubes de Landsat 5 con 6 bandas.

Sibate, Imagen LandSat 8 2015, utilizando 7 bandas

library(raster)
landsat8 = stack(b2,b3,b4,b5,b6,b7)
names(landsat8) = c ('blue', 'green', 'red', 'NIR', 'SWIR1','SWIR2')
plot(landsat8)

Pregunta 1 : Haga una trama compuesta de 3 bandas de falso color de `` landsat8 ’’.

landsat8FCC_imagen = stack(b5, b6, b4)
plotRGB(landsat8FCC_imagen, axes=TRUE, stretch="lin", main="Landsat composicion de color falsa - Agua / Tierra")

# recortando a sibate

sibatextend = shapefile('C:/Users/David-PC/Desktop/ESCRITORIO/SIBATE_RURAL/rural.shp')

ext = extent(sibatextend)

# crop landsat by the extent

landsatfcc = crop(landsat8FCC_imagen, ext)

plotRGB(landsatfcc, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

NA
NA

En la clasificacion no supervisada, utilizamos los datos de reflectancia, pero no proporcionamos ningun dato de respuesta (es decir, no identificamos ningun pixel como perteneciente a una clase en particular). Esto puede parecer extraño, pero puede ser util cuando no tenemos mucho conocimiento previo de un area de estudio. O si tiene un amplio conocimiento de la distribucion de las clases de interes de la cobertura del suelo, pero no tiene datos especificos

El algoritmo agrupa pixeles con caracteristicas espectrales similares en grupos.

Obtenga mas informacion sobre K-means y otros algoritmos supervisados sin supervision -aqui- .

Realizaremos una clasificacion no supervisada en un subconjunto espacial de la -ndvi- capa. Aqui hay otra forma de calcular -ndvi-. En este caso, no usamos una funcion separada, sino una notacion algebraica directa.

ndviLandSat8=(landsat8[['NIR']]-landsat8[['red']])/(landsat8[['NIR']]+landsat8[['red']])
plot(ndviLandSat8)

NA
NA

Haremos un -kmeans- agrupamiento de los -ndvi- datos. Primero usamos -crop- para hacer un subconjunto espacial de -ndvi-, para permitir un procesamiento mas rapido (puede seleccionar cualquiera -extent- usando la -drawExtent()- función).

TITULO: CLASIFICACION KMEANS

# recortando a sibate

extent(landsat8)

class      : Extent 
xmin       : 446085 
xmax       : 673515 
ymin       : 362985 
ymax       : 595215 

landsat8crop2015 = crop(ndviLandSat8, ext)
landsat8crop2015

class      : RasterLayer 
dimensions : 617, 351, 216567  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 576795, 587325, 484005, 502515  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -0.1625031, 0.6681765  (min, max)

plot(landsat8crop2015)

# convertir el raster a vector / matriz
nr=getValues(landsat8crop2015)
str(nr)

 num [1:216567] 0.471 0.441 0.435 0.396 0.378 ...

Tenga en cuenta que -getValues- convirtio el -ndvi- RasterLayer en una matriz (matriz). Ahora realizaremos el -kmeans- agrupamiento en la matriz e inspeccionaremos la salida.

# Es importante configurar el generador de semillas porque `kmeans` inicia los centros en ubicaciones aleatorias, para eso utilizaremos set.seed

set.seed (99)

# Queremos crear 10 grupos, permitir 500 iteraciones, comenzar con 5 conjuntos aleatorios utilizando el método "Lloyd"

kmncluster = kmeans(na.omit(nr), centers = 10, iter.max = 500, nstart = 5, algorithm = "Lloyd")

# kmeans devuelve un objeto de la clase "kmeans"

str(kmncluster)

List of 9
 $ cluster     : int [1:216567] 10 4 4 8 8 8 8 4 10 8 ...
 $ centers     : num [1:10, 1] 0.181 0.285 0.548 0.425 0.236 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "1" "2" "3" "4" ...
  .. ..$ : NULL
 $ totss       : num 3096
 $ withinss    : num [1:10] 6.4 6 7.92 6.11 5.84 ...
 $ tot.withinss: num 57.6
 $ betweenss   : num 3039
 $ size        : int [1:10] 19149 32049 9132 28567 25928 34882 5743 33820 6570 20727
 $ iter        : int 103
 $ ifault      : NULL
 - attr(*, "class")= chr "kmeans"

-kmeans-devuelve un objeto con 9 elementos. La longitud del -cluster- elemento dentro -kmncluster- es 216567 que es igual a la longitud de -nr- creado desde -ndvi-. Los valores de celda de -kmncluster\(cluster- rango entre 1 a 10 correspondientes al numero de entrada del cluster que proporcionamos en la -kmeans- funcion. -kmncluster\)cluster- indica la etiqueta del cluster para el pixel correspondiente. Necesitamos convertir los -kmncluster$cluster- valores nuevamente a RasterLayer de la misma dimension que -ndvi-.

# Use el objeto ndvi (en este caso se llama "landsat8crop2015") para establecer los valores del cluster en un nuevo raster

knr = setValues(landsat8crop2015, kmncluster$cluster)

# Tambien puedes hacerlo asi

knr= raster(landsat8crop2015)
values(knr) = kmncluster$cluster
knr

class      : RasterLayer 
dimensions : 617, 351, 216567  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 576795, 587325, 484005, 502515  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +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 que grupo (1-10) pertenece a que clase de cobertura terrestre (y si pertenece a una clase que reconoceriamos). Puede averiguarlo trazandolos lado a lado con capas de referencia y usando un color unico para cada grupo.

#Defina un vector de color para 10 grupos (mas informacion sobre como configurar el color mas adelante)

mycolor = c ("#fef65b","#ff0000", "#daa520","#0000ff","#0000ff","#00ff00","#cbbeb5",
             "#c3ff5b", "#ff7373", "#00ff00", "#808080")
par(mfrow = c(1,2))
plot(landsat8crop2015, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')
plot(knr, main = 'Unsupervised classification', col = mycolor )

Si bien para otros fines suele ser mejor definir mas clases (y posiblemente fusionar clases mas adelante), una clasificacion simple como esta podria ser util, por ejemplo, fusionar los grupos 4 y 5 para construir una mascara de agua para el año 2015.

Puedes cambiar los colores en mi -mycolor-. Obtenga mas informacion sobre la seleccion de colores en R -aqui- y -aqui- .

Pregunta 2 : Grafique el RGB de 3 bandas de “landsat8” para el subconjunto (extensión “ext”) y el resultado de la agrupacion “kmeans” lado a lado y haga una tabla de cobertura del suelo para el uso del suelo Etiquetas para los grupos. Por ejemplo, los grupos 4 y 5 son agua.

landsat8p2 = stack(b6, b3, b2)
#plotRGB(landsat8p2, axes=TRUE, stretch="lin", main="Landsat composicion de color falsa ")

#crop landsat by the extent

landsatsibatep2 = crop(landsat8p2, ext)

plotRGB(landsatsibatep2, axes = TRUE, stretch = "lin", main = "Landsat FCC - SIBATE")

# convertir el raster a vector / matriz
nr_p2=getValues(landsatsibatep2)
str(nr_p2)

 int [1:216567, 1:3] 13348 12473 12269 13387 14576 14992 14842 12425 11290 9894 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:3] "LC08_L1TP_008057_20150104_20170415_01_T1_B6" "LC08_L1TP_008057_20150104_20170415_01_T1_B3" "LC08_L1TP_008057_20150104_20170415_01_T1_B2"

# Es importante configurar el generador de semillas porque `kmeans` inicia los centros en ubicaciones aleatorias, para eso utilizaremos set.seed

set.seed (99)

# Queremos crear 10 grupos, permitir 500 iteraciones, comenzar con 5 conjuntos aleatorios utilizando el método "Lloyd"

kmncluster_p2 = kmeans(na.omit(nr), centers = 10, iter.max = 500, nstart = 5, algorithm = "Lloyd")

# kmeans devuelve un objeto de la clase "kmeans"

str(kmncluster_p2)

List of 9
 $ cluster     : int [1:216567] 10 4 4 8 8 8 8 4 10 8 ...
 $ centers     : num [1:10, 1] 0.181 0.285 0.548 0.425 0.236 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "1" "2" "3" "4" ...
  .. ..$ : NULL
 $ totss       : num 3096
 $ withinss    : num [1:10] 6.4 6 7.92 6.11 5.84 ...
 $ tot.withinss: num 57.6
 $ betweenss   : num 3039
 $ size        : int [1:10] 19149 32049 9132 28567 25928 34882 5743 33820 6570 20727
 $ iter        : int 103
 $ ifault      : NULL
 - attr(*, "class")= chr "kmeans"

# Use el objeto ndvi (en este caso se llama "landsatsibatep2") para establecer los valores del cluster en un nuevo raster

knr_p2 = setValues(landsatsibatep2, kmncluster_p2$cluster)

# Tambien puedes hacerlo asi

knr_p2= raster(landsatsibatep2)
values(knr_p2) = kmncluster_p2$cluster
knr_p2

class      : RasterLayer 
dimensions : 617, 351, 216567  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 576795, 587325, 484005, 502515  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 1, 10  (min, max)

#Defina un vector de color para 10 grupos (mas informacion sobre como configurar el color mas adelante)

#mycolor_p2 = c (topo.colors(10))

plot(knr_p2, main = 'Unsupervised classification', col = rev(topo.colors(10))) 

NA
NA

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.