Práctica 4 percepción remota

Unsupervised Classification

library(raster)

landsat5 <- stack("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P4/rsdata/rs/centralvalley-2011LT5.tif")

names(landsat5) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')

landsat5

class      : RasterStack 
dimensions : 1230, 1877, 2308710, 6  (nrow, ncol, ncell, nlayers)
resolution : 0.0002694946, 0.0002694946  (x, y)
extent     : -121.9258, -121.42, 37.85402, 38.1855  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : blue, green, red, NIR, SWIR1, SWIR2 

Question 1: Make a 3-band False Color Composite plot of “landsat5”.

landsatRGB <- landsat5[[c(3,2,1)]] # Creamos un objeto con las capas RGB de landsat

# Representamos True color y False color para comparar

par(mfrow = c(1,2)) 
plotRGB(landsatRGB, r=3, g=2, b=1, axes=TRUE, stretch="lin", main="Landsat5 True Color Composite")
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Landsat5 False Color Composite")

# Realizaremos una clasificación no supervisada en un subconjunto espacial de la capa ndvi

ndvi <- (landsat5[['NIR']] - landsat5[['red']]) / (landsat5[['NIR']] + landsat5[['red']])

ndvi

class      : RasterLayer 
dimensions : 1230, 1877, 2308710  (nrow, ncol, ncell)
resolution : 0.0002694946, 0.0002694946  (x, y)
extent     : -121.9258, -121.42, 37.85402, 38.1855  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -0.4932148, 0.8103668  (min, max)

kmeans classification

# Creamos una extensión para recortar la capa ndvi

e <- extent(-121.807, -121.725, 38.004, 38.072)

# Se recorta la imagen Landsat con la extensión

ndvi <- crop(ndvi, e)

ndvi # Si comparamos con el antiguo NVDI se ve que ha sido recortado (menores dimensiones)

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)

# convertir el ráster a vector / matriz

nr <- getValues(ndvi)
str(nr)

 num [1:76608] 0.245 0.236 0.272 0.277 0.277 ...

# Es importante configurar el generador de semillas porque "kmeans" inicia los centros en ubicaciones aleatorias

set.seed (99)

# 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 devuelve un objeto de la 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"

# Se usa el objeto ndvi para establecer los valores del clúster en un nuevo ráster

knr <- setValues(ndvi, kmncluster$cluster)

# También se puede hacer de esta forma

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)

head(knr)

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1  4 4 3 3 3 3 3 4 4  4  4  4  4  4  4  4  4  4  4  4
2  4 4 3 3 3 3 3 4 4  4  4  4  4  4  4  4  4  4  4  4
3  4 3 3 3 3 3 3 4 4  4  4  4  4  4  4  4  4  4  4 10
4  3 3 3 3 3 3 3 4 4  4  4  3  4  4  4  4  4  4  4  4
5  3 3 4 4 4 3 3 4 4  4  4  3  4  4  4  4  4  4  4  4
6  4 4 4 4 4 4 4 4 4  4  4  4  4  4  4  3  4  4  4  4
7  4 4 4 4 4 4 4 4 4  4  4  4  4  4  4  3  3  3  3  4
8  4 4 4 4 4 4 4 4 4  4  3  3  3  7  7  7  7  7  7  7
9  3 3 3 4 4 4 4 4 3  3  3  3  7  9  9  9  9  7  7  7
10 3 3 3 3 3 3 3 4 3  3  3  7  7  9  9  9  9  9  7  7

# Se define un vector de color para los 10 clústers

mycolor <- c("#fef65b","#ff0000", "#daa520","#0000ff","#0000ff","#00ff00","#cbbeb5",
             "#c3ff5b", "#ff7373", "#00ff00", "#808080")

# Se representa el NDVI y la clasificación supervisada para comprobar los grupos creados

par(mfrow = c(1,2))
plot(ndvi, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')
plot(knr, main = 'Unsupervised classification', col = mycolor )

Question 2: Plot 3-band RGB of “landsat5” for the subset (extent “e”) and result of “kmeans” clustering side-by-side and make a table of land-use land-cover labels for the clusters. E.g. cluster 4 and 5 are water.

# Se recorta la imagen Landsat5 con la extensión "e"

Landsat5_crop <- crop(landsat5, e)

Landsat5_crop

class      : RasterBrick 
dimensions : 252, 304, 76608, 6  (nrow, ncol, ncell, nlayers)
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      :        blue,       green,         red,         NIR,       SWIR1,       SWIR2 
min values : 0.084590562, 0.055765331, 0.040395390, 0.034050025, 0.010290771, 0.006807715 
max values :   0.3684652,   0.5407562,   0.5847761,   0.5988820,   0.5135725,   0.5523322 

# Comenzamos con la banda 1 (Azul) de Landsat5

L_Blue <- Landsat5_crop[[1]]

L_Blue

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      : blue 
values     : 0.08459056, 0.3684652  (min, max)

# convertir el ráster a vector / matriz

matriz_Blue <- getValues(L_Blue)
str(matriz_Blue)

 num [1:76608] 0.136 0.139 0.139 0.147 0.147 ...

# Establecemos semillas 

set.seed (99)

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

Cluster_Blue <- kmeans(na.omit(matriz_Blue), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")

# kmeans devuelve un objeto de la clase "kmeans"

str(Cluster_Blue)

List of 9
 $ cluster     : int [1:76608] 9 10 10 10 10 10 9 9 9 9 ...
 $ centers     : num [1:10, 1] 0.111 0.192 0.105 0.164 0.262 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "1" "2" "3" "4" ...
  .. ..$ : NULL
 $ totss       : num 30.8
 $ withinss    : num [1:10] 0.0889 0.1003 0.0555 0.0878 0.1034 ...
 $ tot.withinss: num 0.831
 $ betweenss   : num 30
 $ size        : int [1:10] 21214 771 13839 1845 200 99 17248 9229 7297 4866
 $ iter        : int 136
 $ ifault      : NULL
 - attr(*, "class")= chr "kmeans"

# Se usa el objeto L_Blue recortado para establecer los valores del clúster en un nuevo ráster

Lay_Blue <- setValues(L_Blue, Cluster_Blue$cluster)

Lay_Blue

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      : blue 
values     : 1, 10  (min, max)

head(Lay_Blue)

    1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
1   9 10 10 10 10 10  9  9  9  9  9  9  9  9  9  9 10 10 10 10
2  10 10 10 10 10 10  9  9  9  9  9  9 10 10 10 10 10 10 10 10
3  10  9  9  9  9  9  9  9  9  9  9  9 10 10 10 10 10 10 10 10
4  10 10  9  9  9  9  9  9  9  9  7  9 10 10 10 10 10 10 10 10
5  10 10 10  9  9  9  9  9  9  9  9 10 10 10 10 10  9  9  9  9
6   4  4 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
7  10  4 10  4  4  4  4  4  4  4 10 10 10 10 10 10 10 10 10 10
8  10 10 10 10 10  4  4  4  4  4 10 10 10 10 10  9  9  9  9  9
9  10 10 10 10 10 10  4  4 10 10 10 10  9  9  9  9  9  9  9  9
10 10 10 10 10 10 10 10 10 10 10 10 10  9  9  9  9  9  9  9  9

# Banda 2 (Verde)

L_Green <- Landsat5_crop[[2]]

# convertir el ráster a vector / matriz

matriz_Green <- getValues(L_Green)
str(matriz_Green)

 num [1:76608] 0.139 0.14 0.143 0.144 0.144 ...

# Establecemos semillas 

set.seed (99)

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

Cluster_Green <- kmeans(na.omit(matriz_Green), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")

# kmeans devuelve un objeto de la clase "kmeans"

str(Cluster_Green)

List of 9
 $ cluster     : int [1:76608] 6 6 6 6 6 6 6 6 6 6 ...
 $ centers     : num [1:10, 1] 0.095 0.199 0.1595 0.1073 0.0846 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "1" "2" "3" "4" ...
  .. ..$ : NULL
 $ totss       : num 60.6
 $ withinss    : num [1:10] 0.193 0.269 0.238 0.22 0.128 ...
 $ tot.withinss: num 1.95
 $ betweenss   : num 58.6
 $ size        : int [1:10] 19215 1059 3648 19585 12313 7599 9887 194 89 3019
 $ iter        : int 65
 $ ifault      : NULL
 - attr(*, "class")= chr "kmeans"

# Se usa el objeto L_Green recortado para establecer los valores del clúster en un nuevo ráster

Lay_Green <- setValues(L_Green, Cluster_Green$cluster)

Lay_Green

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      : green 
values     : 1, 10  (min, max)

head(Lay_Green)

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1  6 6 6 6 6 6 6 6 6  6  6  6  6  6  6  6  6  3  3  3
2  6 6 6 6 6 6 6 6 6  6  6  6  6  6  6  6  3  3  3  3
3  6 6 6 6 6 6 6 6 6  6  6  6  3  3  3  3  3  3  3  3
4  6 6 6 6 6 6 6 7 7  7  7  6  3  3  3  3  3  6  6  6
5  3 6 3 6 6 6 7 7 6  6  6  6  6  6  6  3  3  6  6  6
6  3 3 3 3 3 3 3 3 3  3  3  3  3  3  3  3  3  3  3  3
7  3 3 3 3 3 3 3 3 3  3  3  3  3  3  3  3  3  3  3  3
8  3 3 3 3 3 3 3 3 3  3  3  3  3  6  6  6  6  6  6  6
9  3 3 3 3 3 3 3 3 3  3  3  3  6  6  6  6  6  6  6  6
10 3 3 3 3 3 3 3 3 3  3  3  3  6  7  7  7  6  6  6  6

# Banda 3 (Roja)

L_Red <- Landsat5_crop[[3]]

# convertir el ráster a vector / matriz

matriz_Red <- getValues(L_Red)
str(matriz_Red)

 num [1:76608] 0.15 0.153 0.157 0.159 0.159 ...

# Establecemos semillas 

set.seed (99)

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

Cluster_Red <- kmeans(na.omit(matriz_Red), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")

# kmeans devuelve un objeto de la clase "kmeans"

str(Cluster_Red)

List of 9
 $ cluster     : int [1:76608] 5 5 8 8 8 8 8 8 5 5 ...
 $ centers     : num [1:10, 1] 0.2012 0.0848 0.0587 0.4989 0.1419 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "1" "2" "3" "4" ...
  .. ..$ : NULL
 $ totss       : num 114
 $ withinss    : num [1:10] 0.415 0.255 0.231 0.214 0.285 ...
 $ tot.withinss: num 2.87
 $ betweenss   : num 111
 $ size        : int [1:10] 1938 17063 7546 89 6824 9985 242 4517 19849 8555
 $ iter        : int 67
 $ ifault      : NULL
 - attr(*, "class")= chr "kmeans"

# Se usa el objeto La_Red recortado para establecer los valores del clúster en un nuevo ráster

Lay_Red <- setValues(L_Red, Cluster_Red$cluster)

Lay_Red

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      : red 
values     : 1, 10  (min, max)

# Comprobamos que cada capa tiene grupos distintos

head(Lay_Blue)

    1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
1   9 10 10 10 10 10  9  9  9  9  9  9  9  9  9  9 10 10 10 10
2  10 10 10 10 10 10  9  9  9  9  9  9 10 10 10 10 10 10 10 10
3  10  9  9  9  9  9  9  9  9  9  9  9 10 10 10 10 10 10 10 10
4  10 10  9  9  9  9  9  9  9  9  7  9 10 10 10 10 10 10 10 10
5  10 10 10  9  9  9  9  9  9  9  9 10 10 10 10 10  9  9  9  9
6   4  4 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
7  10  4 10  4  4  4  4  4  4  4 10 10 10 10 10 10 10 10 10 10
8  10 10 10 10 10  4  4  4  4  4 10 10 10 10 10  9  9  9  9  9
9  10 10 10 10 10 10  4  4 10 10 10 10  9  9  9  9  9  9  9  9
10 10 10 10 10 10 10 10 10 10 10 10 10  9  9  9  9  9  9  9  9

head(Lay_Green)

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1  6 6 6 6 6 6 6 6 6  6  6  6  6  6  6  6  6  3  3  3
2  6 6 6 6 6 6 6 6 6  6  6  6  6  6  6  6  3  3  3  3
3  6 6 6 6 6 6 6 6 6  6  6  6  3  3  3  3  3  3  3  3
4  6 6 6 6 6 6 6 7 7  7  7  6  3  3  3  3  3  6  6  6
5  3 6 3 6 6 6 7 7 6  6  6  6  6  6  6  3  3  6  6  6
6  3 3 3 3 3 3 3 3 3  3  3  3  3  3  3  3  3  3  3  3
7  3 3 3 3 3 3 3 3 3  3  3  3  3  3  3  3  3  3  3  3
8  3 3 3 3 3 3 3 3 3  3  3  3  3  6  6  6  6  6  6  6
9  3 3 3 3 3 3 3 3 3  3  3  3  6  6  6  6  6  6  6  6
10 3 3 3 3 3 3 3 3 3  3  3  3  6  7  7  7  6  6  6  6

head(Lay_Red)

   1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
1  5 5 8 8 8 8 8 8 5  5  5  5  5  5  5  5  8  8  8  8
2  5 8 8 8 8 8 8 5 8  8  8  5  8  8  8  8  8  1  1  1
3  8 5 5 5 5 5 5 5 8  8  5  5  8  8  8  8  8  8  8  8
4  8 8 8 5 5 5 5 5 5  5  5  5  8  8  8  8  8  8  8  8
5  8 8 8 8 8 8 5 5 5  5  8  8  8  8  8  8  8  8  8  8
6  1 1 1 1 1 1 8 8 8  8  8  8  8  8  8  8  8  8  8  8
7  1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1
8  8 1 1 1 1 1 1 1 1  1  1  8  8  8  8  8  5  5  5  8
9  8 8 8 1 1 1 1 1 8  8  8  8  5  5  5  5  5  5  5  8
10 8 8 8 8 8 8 8 1 8  8  8  8  5 10 10 10  5  5  5  5

# Creamos un raster con las 3 bandas RGB

LandsatRGB <- stack(Lay_Blue, Lay_Green, Lay_Red)
LandsatRGB

class      : RasterStack 
dimensions : 252, 304, 76608, 3  (nrow, ncol, ncell, nlayers)
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 
names      : blue, green, red 
min values :    1,     1,   1 
max values :   10,    10,  10 

# Se representa el Landsat5_crop y la clasificación supervisada para comprobar los grupos creados

plotRGB(LandsatRGB, r = 3, g = 2, b = 1, axes = TRUE, stretch = "lin", main = "Landsat5 True Color Composite")