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

Genere sitios de muestra basados en un ráster de referencia

Extraer valores de celda de datos de Landsat para los sitios de muestra

Entrenar al clasificador usando muestras de entrenamiento

Clasifique los datos de Landsat usando el modelo entrenado

Evaluar la precisión del modelo. Datos de referencia La National Land Cover Database 2011 (NLCD 2011) es un producto de cobertura terrestre para los EE. UU. NLCD es una base de datos de cobertura del suelo basada en Landsat de 30 m que abarca 4 épocas (1992, 2001, 2006 y 2011). NLCD 2011 se basa principalmente en una clasificación de árbol de decisión de datos Landsat de alrededor de 2011.

PRIMERO CARGO LOS DATOS DE REFERENCIA -RASTER

library(raster)
## Loading required package: sp
nlcd <- brick('data/rs/nlcd-L1.tif')
names(nlcd) <- c("nlcd2001", "nlcd2011")

#nlcd raster de referencia; con nombres de 2001 y 2011
# The class names and colors for plotting - A continuación se generan los nombres de clase para el ploter

nlcdclass <- c("Agua", "Areas desarrolladas", "Esteriles", "Bosques", "Arbustos", "Herbaceas", "Plantado/Cultivado", "Humedales")
classdf <- data.frame(classvalue1 = c(1,2,3,4,5,7,8,9), classnames1 = nlcdclass)
# Hex codes of colors # Códigos hexadecimales de colores
classcolor <- c("#5475A8", "#B50000", "#D2CDC0", "#38814E", "#AF963C", "#D1D182", "#FBF65D", "#C8E6F8")
# Now we ratify RAT = "Tabla de atributos de ráster"el ncld2011 
# ratify: define RasterLayer como una variable categórica , el RasterLayer está vinculado a otros valores a través de una "Tabla de atributos de ráster" (RAT). Por lo tanto, los valores de las celdas son un índice, mientras que los valores reales de interés están en la RAT. El RAT es un data.frame. La primera columna en la RAT ("ID") tiene los valores de celda únicos de la capa; esta columna normalmente no se debe cambiar. Las otras columnas pueden ser de cualquier tipo básico (factor, carácter, entero, numérico o lógico). Las funciones documentadas aquí están principalmente disponibles de tal manera que los archivos con una RAT pueden leerse y procesarse;

nlcd2011 <- nlcd[[2]]
nlcd2011 <- ratify(nlcd2011)
rat <- levels(nlcd2011)[[1]]
#
rat$landcover <- nlcdclass
levels(nlcd2011) <- rat

Generar sitios de muestra

En este ejemplo, generaremos los sitios de muestra de capacitación y validación utilizando la referencia de NLCD RasterLayer. Alternativamente, puede usar sitios predefinidos que puede haber recopilado de otras fuentes. Generaremos los sitios de muestra siguiendo un muestreo aleatorio estratificado para asegurar muestras de cada clase LULC.

# Load the training sites locations
# Set the random number generator to reproduce the results
set.seed(99)
# Sampling
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
 #Number of samples in each class
table(samp2011$nlcd2011)
## 
##   1   2   3   4   5   7   8   9 
## 200 200 200 200 200 200 200 200
library("rasterVis")
## Loading required package: lattice
## Loading required package: latticeExtra

Se ve que hay dos variables en samp2011. La cell columna contiene números de celdas de nlcd2011 muestra. la columna nlcd2011ontiene los valores de clase (1-9). Aquí nlcdtiene valores enteros entre 1-9. A menudo encontrará que los nombres de clase se proporcionan como etiquetas de cadena (por ejemplo, agua, cultivo, vegetación). Necesitará ‘volver a etiquetar’ los nombres de clase a enteros o factores si solo se proporcionan etiquetas de cadena antes de usarlas como variable de respuesta en la clasificación. Hay varios enfoques que podrían usarse para convertir estas clases a códigos enteros. Podemos hacer una función que reclasifique las cadenas de caracteres que representan las clases de cobertura del suelo en números enteros basados en los niveles de factores existentes.

nlcd2011Tracemos los sitios de entrenamiento sobre RasterLayer para visualizar la distribución de las ubicaciones de muestreo.

library("lattice")
library("latticeExtra")
plt <- levelplot(nlcd2011, col.regions = classcolor, main = 'Distribucion sitios entrenamiento')
print(plt + layer(sp.points(samp2011, pch = 3, cex = 0.5, col = 1)))

#Extraer valores para sitios 
#Aquí están nuestros datos de Landsat.
library(raster)
landsat5 <- stack('data/rs/centralvalley-2011LT5.tif')
names(landsat5) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
#Una vez que tenemos los sitios, podemos extraer los valores de las celdas de landsat5 RasterStack. Estos valores de banda serán las variables predictoras y los "valores de clase" nlcd2011serán la variable de respuesta.
# Extract the layer values for the locations
sampvals <- extract(landsat5, samp2011, df = TRUE)
# sampvals no longer has the spatial information. To keep the spatial information you use `sp=TRUE` argument in the `extract` function.
# drop the ID column

sampvals <- sampvals[, -1]

# combine the class information with extracted values

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

Entrenar al clasificador Ahora entrenaremos el algoritmo de clasificación utilizando el training2011 conjunto de datos.

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)

n the classification tree plot classvalues are printed at the leaf nodes. You can find the corresponding land use land cover names from the classdf data.frame.

See ?rpart.control to set different parameters for building the model.

You can print/plot more about the cart model created in the previous example. E.g. you can use plotcp(cart) to learn about the cost-complexity (cp argument in rpart).

Clasificar Ahora que tenemos nuestro modelo de clasificación entrenado ( cart), podemos usarlo para hacer predicciones, es decir, para clasificar todas las celdas en landsat5 RasterStack.

Importante Los nombres en el objeto Ráster a clasificar deben coincidir exactamente con los esperados por el modelo. Este será el caso si se utilizó el mismo objeto Ráster (a través de extract) para obtener los valores que se ajustan al modelo.

# Now predict the subset data based on the model; prediction for entire area takes longer time-Ahora prediga los datos del subconjunto según el 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 :
##        ID value
##  from:  1     1
##   to :  8     9

Ahora trace el resultado de la clasificación usando rasterVis. Ver establecerá el classnamespara el classvalues.

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

levelplot(nlcd2011, maxpixels = 1e6,
          col.regions = classcolor,
          scales=list(draw=FALSE),
          main = "nlcd2011")

La National Land Cover Database 2011 (NLCD 2011) es un producto de cobertura terrestre para los EE. UU. NLCD es una base de datos de cobertura del suelo basada en Landsat de 30 m que abarca 4 épocas (1992, 2001, 2006 y 2011)

levelplot(nlcd2011, maxpixels = 1e6,
          col.regions = classcolor,
          scales=list(draw=FALSE),
          main = "nlcd2011")

Evaluación del modelo

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 (típicamente 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
## 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')
    # create a data.frame using the reference and prediction
    x[[k]] <- cbind(test$classvalue, as.integer(pclass))
}

Ahora combine los cinco elementos de la lista en un solo data.frame, usando do.cally calcule una matriz de confusión.

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              Agua Areas desarrolladas Esteriles Bosques Arbustos
##   Agua                 175                   6         0       3        0
##   Areas desarrolladas    2                  90        51       8       10
##   Esteriles              7                  39        82       4       19
##   Bosques                0                   2         1     106       57
##   Arbustos               0                   3         5      59      102
##   Herbaceas              0                   9        36      10       27
##   Plantado/Cultivado     0                   7        11      34       42
##   Humedales             18                  10         6      36       29
##                      predicted
## observed              Herbaceas Plantado/Cultivado Humedales
##   Agua                        0                  7         9
##   Areas desarrolladas        22                 11         6
##   Esteriles                  38                  5         6
##   Bosques                     1                  6        27
##   Arbustos                   12                 12         7
##   Herbaceas                 109                  8         1
##   Plantado/Cultivado         19                 69        18
##   Humedales                   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
## [1] 1600
# number of correctly classified cases per class
diag <- diag(conmat)
# Overall Accuracy
OA <- sum(diag) / n
OA
## [1] 0.4975
## [1] 0.4975

Kappa:

# observed (true) cases per class
rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
# predicted cases per class
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
kappa
## [1] 0.4257143
## [1] 0.4257143

Productor y precisión del usuario

# Producer accuracy
PA <- diag / colsums
# User accuracy
UA <- diag / rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
outAcc
##                     producerAccuracy userAccuracy
## Agua                       0.8663366        0.875
## Areas desarrolladas        0.5421687        0.450
## Esteriles                  0.4270833        0.410
## Bosques                    0.4076923        0.530
## Arbustos                   0.3566434        0.510
## Herbaceas                  0.5291262        0.545
## Plantado/Cultivado         0.4569536        0.345
## Humedales                  0.4598540        0.315