Aquí exploramos la clasificación supervisada. Existen varios algoritmos de clasificación supervisada, y la elección del algoritmo puede afectar a los resultados. Aquí exploramos dos algoritmos relacionados (CART y RandomForest).

En la clasificación supervisada, tenemos conocimiento previo sobre algunos de los tipos de cubierta terrestre a través de, por ejemplo, el trabajo de campo, los datos espaciales de referencia o la interpretación de imágenes de alta resolución (como las disponibles en los mapas de Google). Se identifican sitios específicos en la zona de estudio que representan ejemplos homogéneos de estos tipos de cubierta terrestre conocidos. Esas zonas se denominan comúnmente lugares de capacitación porque las propiedades espectrales de esos lugares se utilizan para capacitar al algoritmo de clasificación.

En los siguientes ejemplos se utiliza un clasificador de árboles de clasificación y regresión (CART) (Breiman et al. 1984) (lectura adicional para predecir las clases de cubierta terrestre en el área de estudio.

Realizaremos los siguientes pasos: - Generar sitios de muestra basados en un raster de referencia - Extraer los valores de las células de los datos del Landsat para los lugares de muestra - Entrenar el clasificador usando muestras de entrenamiento - Clasificar los datos del Landsat usando el modelo entrenado - Evaluar la precisión del modelo

Datos de referencia

La Base de Datos Nacional de Cubierta Terrestre 2011 (NLCD 2011) es un producto de cubierta terrestre para los Estados Unidos. La NLCD es una base de datos sobre la cubierta terrestre basada en Landsat de 30 m que abarca 4 épocas (1992, 2001, 2006 y 2011). La NLCD 2011 se basa principalmente en una clasificación de árbol de decisión de los datos de Landsat de alrededor de 2011.

Puede encontrar los nombres de las clases en el NLCD 2011 (aquí)[https://www.mrlc.gov/nlcd11_leg.php]. Tiene dos pares de valores y nombres de clase que corresponden a los niveles de uso de la tierra y el sistema de clasificación de la cubierta terrestre. Estos niveles suelen representar el nivel de complejidad, siendo el nivel I el más simple con amplias categorías de uso de la tierra y cubierta terrestre. Lea este informe de Anderson et al para aprender más sobre este sistema de clasificación de uso de la tierra y cobertura de la tierra.

library(raster)
## Loading required package: sp
nlcd <- brick('data/rs/nlcd-L1.tif')
names(nlcd) <- c("nlcd2001", "nlcd2011")
# Los nombres de las clases y los colores para el trazado
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)
# Códigos hexagonales de colores
classcolor <- c("#5475A8", "#B50000", "#D2CDC0", "#38814E", "#AF963C", "#D1D182", "#FBF65D", "#C8E6F8")
# Ahora ratificamos (RAT = "Raster Attribute Table") el ncld2011 (definir RasterLayer como una variable categórica). Esto es útil para el trazado.
nlcd2011 <- nlcd[[2]]
nlcd2011 <- ratify(nlcd2011)
rat <- levels(nlcd2011)[[1]]
#
rat$landcover <- nlcdclass
levels(nlcd2011) <- rat

Hicimos muchas cosas aquí. Da un paso atrás y lee más sobre la ratificación.

Nota No hay ninguna clase con valor 6.

Generar sitios de muestra Como discutimos en la clase, los datos de entrenamiento y/o validación pueden provenir de una variedad de fuentes. En este ejemplo generaremos los sitios de muestra de capacitación y validación utilizando la referencia RasterLayer del NLCD. Alternativamente, se pueden utilizar sitios predefinidos que pueden haber sido recogidos de otras fuentes. Generaremos los sitios de muestra siguiendo un muestreo aleatorio estratificado para asegurar las muestras de cada clase de LULC.

# Cargar las ubicaciones de los lugares de entrenamiento
# Establecer 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

Puedes ver que hay dos variables en samp2011. La columna de la celda contiene los números de la celda de nlcd2011 muestreada. La columna de nlcd2011 contiene los valores de la clase (1-9). Dejaremos la columna de la célula más tarde.

Aquí nlcd tiene valores enteros entre 1-9. A menudo encontrará que los nombres de clase se proporcionan como etiquetas de cadenas (por ejemplo, agua, cultivo, vegetación). Necesitará “re-etiquetar” los nombres de clase a números enteros o factores si sólo se suministran etiquetas de cadena antes de usarlos como variable de respuesta en la clasificación. Hay varios enfoques que podrían utilizarse para convertir estas clases en códigos enteros. Podemos hacer una función que reclasifique las cadenas de caracteres que representan las clases de cobertura terrestre en números enteros basados en los niveles de factores existentes.

Vamos a trazar los sitios de entrenamiento sobre el nlcd2011 RasterLayer para visualizar la distribución de los lugares de muestreo.

En caso de ser necesario ejecutar install.packages(“rasterVis”)

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

rasterVis ofrece un trazado más avanzado (enrejado/enrejado) de los objetos Raster*. Por favor, instale el paquete si no está disponible para su máquina.

Extraer los valores de los sitios

Aquí están nuestros datos del Landsat.

landsat5 <- stack('data/rs/centralvalley-2011LT5.tif')
names(landsat5) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')

Una vez que tengamos los sitios, podemos extraer los valores de las células de landsat5 RasterStack. Estos valores de banda serán las variables predictoras y los “valores de clase” de nlcd2011 serán la variable de respuesta.

# Extraer los valores de las capas para las ubicaciones
sampvals <- extract(landsat5, samp2011, df = TRUE)
# Sampvals ya no tiene la información espacial. Para mantener la información espacial se utiliza el argumento `sp=TRUE` en la función `extracto`.
# Suelta la columna ID
sampvals <- sampvals[, -1]
# Combinar la información de la clase con los valores extraídos
sampdata <- data.frame(classvalue = samp2011@data$nlcd2011, sampvals)

Entrena el clasificador Ahora entrenaremos el algoritmo de clasificación usando el conjunto de datos de training2011.

library(rpart)
# Train the model
cart <- rpart(as.factor(classvalue)~., data=sampdata, method = 'class', minsplit = 5)
# print(model.class)
# Plot the trained classification tree
plot(cart, uniform=TRUE, main="Classification Tree")
text(cart, cex = 0.8)

En el árbol de clasificación se imprimen valores de clase en los nodos de la hoja. Los nombres de las cubiertas de tierra correspondientes se encuentran en el marco de datos de la clasificación.

Véase ?rpart.control para establecer diferentes parámetros para la construcción del modelo.

Puede imprimir/trazar más sobre el modelo de carro creado en el ejemplo anterior. Por ejemplo, puede utilizar plotcp(cart) para conocer la complejidad de los costes (argumento cp en rpart).

Clasificar

Ahora que tenemos nuestro modelo de clasificación entrenado (cart), podemos usarlo para hacer predicciones, es decir, para clasificar todas las células del landsat5 RasterStack.

Importante Los nombres en el objeto Raster a clasificar deben coincidir exactamente con los esperados por el modelo. Este será el caso si el mismo objeto Raster fue utilizado (vía extracto) para obtener los valores que se ajustan al modelo.

# Ahora predice los datos del subconjunto basado en el modelo; la predicción para toda el área toma 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 :
##        ID value
##  from:  1     1
##   to :  8     9

Ahora traza el resultado de la clasificación usando rasterVis. See establecerá los nombres de las clases para los valores de las clases.

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")

Pregunta 1:Traza nlcd2011 y pr2011 uno al lado del otro y comenta la exactitud de la predicción (por ejemplo, mezclando entre cultivos, pastos, praderas y arbustos).

Es posible que tenga que seleccionar más muestras y utilizar variables de predicción adicionales. La elección del clasificador también desempeña un papel importante.

Evaluación del modelo Ahora evaluemos la precisión del modelo para tener una idea de cuán preciso puede ser el mapa clasificado. Dos medidas muy utilizadas en la teledetección son la “exactitud global” y el “kappa”. La evaluación de la exactitud se puede realizar utilizando las muestras independientes (validación2011).

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

En caso de ser necesario ejecutar install.packages(“dismo”)

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 computando 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')
    # Crear un marco de datos usando la referencia y la predicción
    x[[k]] <- cbind(test$classvalue, as.integer(pclass))
}

Ahora combina los cinco elementos de la lista en un solo marco de datos, usando do.call y calcula una matriz de confusión.

y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
conmat <- table(y)
# cambiar el nombre de las clases
colnames(conmat) <- classdf$classnames
rownames(conmat) <- classdf$classnames
conmat
##                     predicted
## observed             Water Developed Barren Forest Shrubland Herbaceous
##   Water                175         6      0      3         0          0
##   Developed              2        90     51      8        10         22
##   Barren                 7        39     82      4        19         38
##   Forest                 0         2      1    106        57          1
##   Shrubland              0         3      5     59       102         12
##   Herbaceous             0         9     36     10        27        109
##   Planted/Cultivated     0         7     11     34        42         19
##   Wetlands              18        10      6     36        29          5
##                     predicted
## observed             Planted/Cultivated Wetlands
##   Water                               7        9
##   Developed                          11        6
##   Barren                              5        6
##   Forest                              6       27
##   Shrubland                          12        7
##   Herbaceous                          8        1
##   Planted/Cultivated                 69       18
##   Wetlands                           33       63

Pregunta 2: Comente la clasificación errónea entre las diferentes clases.

Pregunta 3: ¿Puede pensar en formas de mejorar la precisión.

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

Exactitud general:

# Número de casos
n <- sum(conmat)
n
## [1] 1600
# Número de casos correctamente clasificados por clase
diag <- diag(conmat)
# Precisión general
OA <- sum(diag) / n
OA
## [1] 0.4975

Kappa:

# Casos observados (verdaderos) por clase
rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
# Casos pronosticados por clase
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
kappa
## [1] 0.4257143

Precisión del productor y del usuario

# Precisión del productor
PA <- diag / colsums
# Precisión del usuario
UA <- diag / rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
outAcc

Pregunta 4: Realice la clasificación usando los clasificadores de Random Forest del paquete “RandomForest”.

Pregunta 5: Traza los resultados de rpart y del clasificador aleatorio de bosques uno al lado del otro.

Pregunta 6 (opcional): Repita los pasos para el año 2001 usando Random Forest. Utilice los datos de la imagen compuesta sin nubes/centralvalley-2001LE7.tif. Estos son los datos de Landsat 7 . Utilice como datos de referencia la National Land Cover Database 2001 (NLCD 2001) para el subconjunto del Valle Central de California.*

Pregunta 7 (opcional): Hemos entrenado a los clasificadores utilizando 200 muestras para cada clase. Investigar el efecto del tamaño de la muestra en la clasificación. Repita los pasos con diferentes subconjuntos, por ejemplo, un tamaño de muestra de 150, 100, 50 por clase, y compare los resultados. Utilizar las mismas muestras de reserva para la evaluación del modelo.