library(raster)
## Loading required package: sp
landsat5 <- stack('data/rs/centralvalley-2011LT5.tif')
names(landsat5) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
library(sp)

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

plotRGB(landsat5, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Landsat5 Falso Color")

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 en particular). Esto puede parecer extraño, pero puede ser útil cuando no tenemos mucho conocimiento previo de un área de estudio. O si tiene un amplio conocimiento de la distribución de las clases de interés de la cobertura del suelo, pero no tiene datos específicos sobre el terreno.

El algoritmo agrupa píxeles con características espectrales similares en grupos.

Obtenga más información sobre K-means y otros algoritmos supervisados sin supervisión aquí .

Realizaremos una clasificación no supervisada en un subconjunto espacial de la ndvicapa. Aquí hay otra forma de calcular ndvi. En este caso, no usamos una función separada, sino una notación algebraica directa

ndvi <- (landsat5[['NIR']] - landsat5[['red']]) / (landsat5[['NIR']] + landsat5[['red']])
plot(ndvi, main="Landsat5 NDVI")

aremos un kmeansagrupamiento de los ndvidatos. Primero usamos croppara hacer un subconjunto espacial de ndvi, para permitir un procesamiento más rápido (puede seleccionar cualquiera extentusando la drawExtent() función).

kmeans clasificación

# Extent to crop ndvi layer; hacer el extent para cortar la img
e <- extent(-121.807, -121.725, 38.004, 38.072)
# crop landsat by the extent- aplicar el extent hecho para hacer recorte de la img
ndvi2 <- crop(ndvi, e)
ndvi2
## 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)
plot(ndvi2)

ndvi2
## 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)

convert the raster to vecor/matrix - convertir raster en matriz

nr <- getValues(ndvi2) Se puede consultar la estructura de la columna animales con str(nr) ## num [1:76608] 0.245 0.236 0.272 0.277 0.277 …

nr <- getValues(ndvi2)
str(nr)
##  num [1:76608] 0.245 0.236 0.272 0.277 0.277 ...

Tenga en cuenta que getValues convirtió el ndviRasterLayer en una matriz (matriz). Ahora realizaremos el kmeans agrupamiento en la matriz e inspeccionaremos la salida.

nr esigual a la matriz que hicimos de ndvi con getvalues

primero configurar el generador de semillas: La función set.seed se usa para establecer la semilla aleatoria para todas las funciones de aleatorización. Si está usando R para crear una aleatorización que desea reproducir, primero debe usar set.seed .

set.seed(99)

para kmeans metodo, # Queremos crear 10 grupos, permitir 500 iteraciones, comenzar con 5 conjuntos aleatorios usando el método “Lloyd”

kmncluster <- kmeans(na.omit(nr), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")
# kmeans returns an object of class "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"

kmeansdevuelve un objeto con 9 elementos. La longitud del clusterelemento dentro kmnclusteres 76608 que es igual a la longitud de nr (matrizcreada) creado desde ndvi. Los valores de celda de kmncluster\(cluster rango entre 1 a 10 correspondientes al número de entrada del clúster que proporcionamos en la kmeansfunción. kmncluster\)clusterindica la etiqueta del clúster para el píxel correspondiente.

Necesitamos convertir los kmncluster$cluster valores nuevamente a RasterLayer de la misma dimensión que ndvi2.

# Use the ndvi object to set the cluster values to a new raster - set values lo que hace es Asignar valores (nuevos) a un objeto Raster

knr <- setValues(ndvi2, kmncluster$cluster)
# You can also do it like this
knr <- raster(ndvi2)
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 knres un RasterLayer, pero no sabemos qué grupo (1-10) pertenece a qué clase de cobertura terrestre (y si pertenece a una clase que reconoceríamos). Puede averiguarlo trazándolos lado a lado con capas de referencia y usando un color único para cada grupo.

# Define a color vector for 10 clusters (learn more about setting the color later)
mycolor <- c("#fef65b","#ff0000", "#daa520","#0000ff","#0000ff","#00ff00","#cbbeb5",
             "#c3ff5b", "#ff7373", "#00ff00", "#808080")


plot(ndvi2, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')

plot(knr, main = 'Unsupervised classification', col = mycolor )

Si bien para otros fines suele ser mejor definir más clases (y posiblemente fusionar clases más adelante), una clasificación simple como esta 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 mycolor. Obtenga más información sobre la selección de colores en R aquí y aquí .

Pregunta 2 : Trace el RGB de 3 bandas de landsat5 para el subconjunto (extensión e) y el resultado de la agrupación de 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.

#respuesta parte 1 RGB del extent
#aplicar el extent hecho para hacer recorte de la img RGB de 3 bandas de `` landsat5

landsat5b <- crop(landsat5, e)
plotRGB(landsat5b, r=4, g=2, b=1, axes=TRUE, stretch="lin", main="Landsat5B 421")

#respuesta parte 2 KMEANS 421 - LANDSAT5B Extent
#aplicar el extent hecho para hacer recorte de la img RGB de 3 bandas de `` landsat5

#primero convierto raster en matriz
nrejem <- getValues(landsat5b)
str(nrejem)
##  num [1:76608, 1:6] 0.136 0.139 0.139 0.147 0.147 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:6] "blue" "green" "red" "NIR" ...
#luego hago el set.seed y aleatorizacion
set.seed(99)
kmnclusterejem <- kmeans(na.omit(nrejem), centers = 5, iter.max = 500, nstart = 5, algorithm="Lloyd")
str(kmnclusterejem)
## List of 9
##  $ cluster     : int [1:76608] 2 2 2 2 2 2 2 2 1 1 ...
##  $ centers     : num [1:5, 1:6] 0.128 0.149 0.11 0.111 0.104 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:6] "blue" "green" "red" "NIR" ...
##  $ totss       : num 2016
##  $ withinss    : num [1:5] 43 81.7 19 57.7 32.3
##  $ tot.withinss: num 234
##  $ betweenss   : num 1783
##  $ size        : int [1:5] 13704 8026 25866 19202 9810
##  $ iter        : int 62
##  $ ifault      : NULL
##  - attr(*, "class")= chr "kmeans"
knrejem <- raster(landsat5b)
values(knrejem) <- kmnclusterejem$cluster
knrejem
## 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, 5  (min, max)

finalmente para hacer o definir los valores se utiliza

mycolor <- c('green', 'yellow', 'blue', 'burlywood', 'cyan')

plot(knrejem, col = mycolor, main = 'CLASIFICACION NO SUPERVISADA')

legend("topright",
       legend = c("cultivo 1", "cultivo 2", "Agua", "urbanizados", "pastos humedos"),
       fill = c("green", "yellow", "blue", "cyan", "burlywood"),
       border = FALSE,
       bty = "n") # turn off legend border

library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(ggplot2)