Practica 5 Percepcion remota

Supervised Classification

setwd("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P5") # Cambiar el directorio de trabajo

getwd()

[1] "C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P5"

library(raster)

Reference data

nlcd <- brick("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P5/rsdata/rs/nlcd-L1.tif")

names(nlcd) <- c("nlcd2001","nlcd2011")

nlcd

class      : RasterBrick 
dimensions : 1230, 1877, 2308710, 2  (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 
source     : C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P5/rsdata/rs/nlcd-L1.tif 
names      : nlcd2001, nlcd2011 
min values :        1,        1 
max values :        9,        9 

head(nlcd)

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

# Los nombres de clase y colores para plotear

nlcdclass <- c("Water", "Developed", "Barren", "Forest", "Shrubland", "Herbaceous", "Planted/Cultivated", "Wetlands")

classdf <- data.frame(classvalue1 = c(1,2,3,4,5,7,8,9), classnames1 = nlcdclass)
classdf

# Códigos hexadecimales de colores

classcolor <- c("#5475A8", "#B50000", "#D2CDC0", "#38814E", "#AF963C", "#D1D182", "#FBF65D", "#C8E6F8")
classcolor

[1] "#5475A8" "#B50000" "#D2CDC0" "#38814E" "#AF963C" "#D1D182" "#FBF65D" "#C8E6F8"

# Ahora ratificamos (RAT = "Tabla de atributos de ráster") el ncld2011 (definimos RasterLayer como una variable categórica). Esto es útil para trazar.

nlcd2011 <- nlcd[[2]]
nlcd2011

class      : RasterLayer 
band       : 2  (of  2  bands)
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     : C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P5/rsdata/rs/nlcd-L1.tif 
names      : nlcd2011 
values     : 1, 9  (min, max)

head(nlcd2011)

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

# Con ratify() se puede ver la "Tabla de atributos de ráster" (RAT) como un marco de datos en un nuevo espacio denominado attributes

nlcd2011 <- ratify(nlcd2011)

# Podemos apreciar que ahora el objeto nlcd2011 tiene el apartado "attributes"

nlcd2011

class      : RasterLayer 
band       : 2  (of  2  bands)
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     : C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P5/rsdata/rs/nlcd-L1.tif 
names      : nlcd2011 
values     : 1, 9  (min, max)
attributes :

NA

# Lo que buscamos es asociar los niveles de "nlcd2011" con las categorias (landcover), para lograrlo creamos un objeto "rat". Introducimos en el objeto rat los niveles (valores que pueden tener los pixeles) de nlcd2011

rat <- levels(nlcd2011)[[1]]
rat

# Adjudicamos a cada ID (nivel de nlcd2011) una categoria (landcover) con el objeto nlcdclass

rat$landcover <- nlcdclass
rat

levels(nlcd2011) <- rat
levels(nlcd2011)

[[1]]

NA

Generate sample sites

# Cargue las ubicaciones de los sitios de entrenamiento
# Configure el generador de números aleatorios para reproducir los resultados

set.seed(99)

# Muestreo

samp2011 <- sampleStratified(nlcd2011, size = 200, na.rm = TRUE, sp = TRUE)
samp2011

class       : SpatialPointsDataFrame 
features    : 1600 
extent      : -121.9257, -121.4225, 37.85415, 38.18536  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 2
names       :    cell, nlcd2011 
min values  :     413,        1 
max values  : 2307837,        9 

table(samp2011$nlcd2011)

  1   2   3   4   5   7   8   9 
200 200 200 200 200 200 200 200 

# Tracemos los sitios de capacitación sobre nlcd2011 RasterLayer para visualizar la distribución de las ubicaciones de muestreo.

library(rasterVis)
plt <- levelplot(nlcd2011, col.regions = classcolor, main = 'Distribution of Training Sites')
print(plt + layer(sp.points(samp2011, pch = 3, cex = 0.5, col = 1)))

Extract values for sites

# Cargamos los datos Landsat

landsat5 <- stack("C:/Users/USS/Desktop/Pedro/2019 - 2020/Percepcion remota/P5/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 

# Extrae los valores de capa para las ubicaciones

sampvals <- extract(landsat5, samp2011, df = TRUE)
sampvals

# sampvals ya no tiene la información espacial. Para mantener la información espacial, utilice el argumento "sp = TRUE" en la función "extract".

# Quitamos la columna ID

sampvals <- sampvals[, -1]
sampvals

# Se combina la información de la clase con los valores extraídos mediante un data.frame

sampdata <- data.frame(classvalue = samp2011@data$nlcd2011, sampvals)
sampdata

Train the classifier

library(rpart)

# Formamos el modelo. La formula que usaremos es tipo ~ ., la cual expresa que intentaremos clasificar classvalue usando a todas las demás variables como predictoras.

cart <- rpart(as.factor(classvalue)~., data=sampdata, method = 'class', minsplit = 5)
cart

n= 1600 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

  1) root 1600 1400 1 (0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12)  
    2) SWIR1< 0.07371388 207   28 1 (0.86 0.0097 0.034 0 0 0 0 0.092) *
    3) SWIR1>=0.07371388 1393 1193 4 (0.015 0.14 0.14 0.14 0.14 0.14 0.14 0.13)  
      6) blue>=0.1180971 595  425 2 (0.0067 0.29 0.28 0.012 0.035 0.27 0.067 0.039)  
       12) SWIR1< 0.1966332 146   57 2 (0.021 0.61 0.25 0.014 0.021 0.0068 0.0068 0.075) *
       13) SWIR1>=0.1966332 449  287 7 (0.0022 0.18 0.29 0.011 0.04 0.36 0.087 0.027)  
         26) blue>=0.1506717 54   22 3 (0 0.3 0.59 0 0 0.093 0.019 0) *
         27) blue< 0.1506717 395  238 7 (0.0025 0.16 0.25 0.013 0.046 0.4 0.096 0.03)  
           54) SWIR1< 0.2871708 302  213 7 (0.0033 0.2 0.28 0.017 0.056 0.29 0.11 0.036)  
            108) green>=0.1254639 182  110 3 (0.0055 0.26 0.4 0.011 0.016 0.19 0.088 0.027) *
            109) green< 0.1254639 120   66 7 (0 0.092 0.12 0.025 0.12 0.45 0.15 0.05) *
           55) SWIR1>=0.2871708 93   25 7 (0 0.065 0.14 0 0.011 0.73 0.043 0.011) *
      7) blue< 0.1180971 798  605 4 (0.021 0.035 0.033 0.24 0.22 0.046 0.2 0.2)  
       14) NIR< 0.2540509 666  475 4 (0.012 0.03 0.033 0.29 0.26 0.047 0.15 0.18)  
         28) NIR>=0.2062028 302  180 4 (0 0.033 0.017 0.4 0.13 0.073 0.17 0.17) *
         29) NIR< 0.2062028 364  232 5 (0.022 0.027 0.047 0.19 0.36 0.025 0.14 0.19)  
           58) blue< 0.09249049 51   17 5 (0 0 0 0.31 0.67 0 0 0.02) *
           59) blue>=0.09249049 313  215 5 (0.026 0.032 0.054 0.17 0.31 0.029 0.16 0.21)  
            118) SWIR2>=0.08096182 232  141 5 (0.013 0.026 0.06 0.17 0.39 0.039 0.19 0.11) *
            119) SWIR2< 0.08096182 81   40 9 (0.062 0.049 0.037 0.17 0.086 0 0.086 0.51) *
       15) NIR>=0.2540509 132   74 8 (0.068 0.061 0.03 0.015 0.053 0.045 0.44 0.29)  
         30) SWIR1>=0.1437143 112   55 8 (0.062 0.062 0.036 0.018 0.062 0.054 0.51 0.2) *
         31) SWIR1< 0.1437143 20    4 9 (0.1 0.05 0 0 0 0 0.05 0.8) *

# Ploteamos (model.class)
# Trazar el árbol de clasificación formado.

par(xpd = FALSE, mar = c(0,1,1,0))
plot(cart, uniform = TRUE, main="Classification Tree", )
text(cart, cex = 0.8)

Classify

# Ahora que tenemos nuestro modelo de clasificación entrenado, podemos usarlo para hacer predicciones, es decir, para clasificar todas las celdas en el raster stack Landsat 5.

# Ahora se predice los datos del subconjunto en función del modelo; la predicción para toda el área lleva más tiempo

pr2011 <- predict(landsat5, cart, type='class')
pr2011

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     : 1, 9  (min, max)
attributes :

NA

# Ahora trace el resultado de la clasificación usando rasterVis. Observar que establecerá los nombres de clase para los valores de clase

pr2011 <- ratify(pr2011)
rat <- levels(pr2011)[[1]]
rat$legend <- classdf$classnames
levels(pr2011) <- rat
levelplot(pr2011, maxpixels = 1e6,
          col.regions = classcolor,
          scales=list(draw=FALSE),
          main = "Decision Tree classification of Landsat 5")

Question 1:Plot “nlcd2011” and “pr2011” side-by-side and comment about the accuracy of the prediction (e.g. mixing between cultivated crops, pasture, grassland and shrubs).

# Quitamos los ejes y la leyenda de colores para una mejor visualización. En gráficos lattice no funciona la funcion par().

library(gridExtra)
plot_nlcd2011 <- levelplot(nlcd2011,
                    maxpixels = 1e6, 
                    xlab= NULL,
                    ylab= NULL,
                    main='nlcd2011',
                    colorkey = FALSE,
                    scales=list(draw=FALSE),
                    col.regions = classcolor,
                    labels = FALSE)
plot_pr2011 <- levelplot(pr2011,
                    maxpixels = 1e6, 
                    xlab= NULL,
                    ylab= NULL,
                    main='pr2011 (predicción)',
                    colorkey = FALSE,
                    scales=list(draw=FALSE),
                    col.regions = classcolor,
                    labels = FALSE)
grid.arrange(plot_nlcd2011, plot_pr2011, ncol=2)

Model evaluation

# Ahora analicemos la precisión del modelo para tener una idea de la precisión del mapa clasificado.

# Para evaluar cualquier modelo, puede usar la validación cruzada k-fold. En esta técnica, los datos utilizados para ajustar el modelo se dividen en k grupos (generalmente 5 grupos). A su vez, uno de los grupos se usará para la prueba del modelo, mientras que el resto de los datos se usará para el entrenamiento del modelo (ajuste).

library(dismo)
set.seed(99)
j <- kfold(sampdata, k = 5, by=sampdata$classvalue)
table(j)

j
  1   2   3   4   5 
320 320 320 320 320 

# Ahora entrenamos y probamos el modelo cinco veces, cada vez calculando una matriz de confusión que almacenamos en una lista.

x <- list()
for (k in 1:5) {
    train <- sampdata[j!= k, ]
    test <- sampdata[j == k, ]
    cart <- rpart(as.factor(classvalue)~., data=train, method = 'class', minsplit = 5)
    pclass <- predict(cart, test, type='class')

# Creamos un data.frame usando la referencia y la predicción
# Cuando se desea unir una tabla a la derecha de otra se utiliza la función "cbind"
   
     x[[k]] <- cbind(test$classvalue, as.integer(pclass))
}

# Ahora combinamos los cinco elementos de la lista en un solo data.frame, usando la función "do.call" y calcule una matriz de confusión.
# Cuando se desea unir una tabla a debajo de otra, se utiliza la función "rbind"

y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
conmat <- table(y)
# change the name of the classes
colnames(conmat) <- classdf$classnames
rownames(conmat) <- classdf$classnames
conmat

                    predicted
observed             Water Developed Barren Forest Shrubland Herbaceous Planted/Cultivated Wetlands
  Water                175         6      0      3         0          0                  7        9
  Developed              2        90     51      8        10         22                 11        6
  Barren                 7        39     82      4        19         38                  5        6
  Forest                 0         2      1    106        57          1                  6       27
  Shrubland              0         3      5     59       102         12                 12        7
  Herbaceous             0         9     36     10        27        109                  8        1
  Planted/Cultivated     0         7     11     34        42         19                 69       18
  Wetlands              18        10      6     36        29          5                 33       63

Calcule la precisión general y la estadística “Kappa”.

Precisión general:

# number of cases
n <- sum(conmat)
n

[1] 1600

# Cantidad de casos correctamente clasificados por clase. La función "diag" nos da la diagonal de la tabla conmat, es decir, los casos en los que las clases de observed coinciden con las de predicted

diag <- diag(conmat)
diag

             Water          Developed             Barren             Forest          Shrubland 
               175                 90                 82                106                102 
        Herbaceous Planted/Cultivated           Wetlands 
               109                 69                 63 

# Precisión general

OA <- sum(diag) / n
OA

[1] 0.4975

Kappa:

# Casos observados (verdaderos) por clase. La función "apply" aplica una función a todos los elementos de una matriz. Con MARGIN = 1 sumamos los elementos de cada fila.

rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
p

             Water          Developed             Barren             Forest          Shrubland 
             0.125              0.125              0.125              0.125              0.125 
        Herbaceous Planted/Cultivated           Wetlands 
             0.125              0.125              0.125 

# Casos predichos por clase. Con MARGIN = 2 sumamos los elementos de cada columna.

colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
kappa

[1] 0.4257143

Producer and user accuracy

# Producer accuracy

PA <- diag / colsums
PA

             Water          Developed             Barren             Forest          Shrubland 
         0.8663366          0.5421687          0.4270833          0.4076923          0.3566434 
        Herbaceous Planted/Cultivated           Wetlands 
         0.5291262          0.4569536          0.4598540 

# User accuracy

UA <- diag / rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
outAcc