Imagen LandSat 8, Sibate 2015 - Clasificacion supervisada

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

TITULO: CLASIFICACION SUPERVISADA

Aqui exploramos la clasificacion supervisada. Existen varios algoritmos de clasificacion supervisados, y la eleccion del algoritmo puede afectar los resultados. Aqui exploramos dos algoritmos relacionados (CART y RandomForest).

En la clasificacion supervisada, tenemos conocimiento previo sobre algunos de los tipos de cobertura del suelo mediante, por ejemplo, trabajo de campo, datos espaciales de referencia o interpretación de imagenes de alta resolucion (como las disponibles en los mapas de Google). Se identifican sitios especificos en el area de estudio que representan ejemplos homogeneos de estos tipos de cobertura terrestre conocidos. Estas areas se conocen comunmente como sitios de entrenamiento porque las propiedades espectrales de estos sitios se usan para entrenar el algoritmo de clasificacion.

Los siguientes ejemplos utilizan un clasificador de arboles de clasificacion y regresion (CART) (Breiman et al. 1984) ( lecturas adicionales para predecir las clases de cobertura del suelo en el area de estudio).

Realizaremos los siguientes pasos:

Genere sitios de muestra basados en un raster 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 precision del modelo.

TITULO: 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.

Puede encontrar los nombres de clase en NCLD 2011 (aquí) -[ https://www.mrlc.gov/nlcd11_leg.php ]-. Tiene dos pares de valores de clase y nombres que corresponden a los niveles de uso del suelo y el sistema de clasificacion de la cobertura del suelo. Estos niveles generalmente representan el nivel de complejidad, siendo el nivel I el mas simple con amplias categorias de cobertura del suelo. Lea este informe de Anderson et al. Para obtener mas informacion sobre este sistema de clasificacion de uso y cobertura de la tierra.

library(raster)

Loading required package: sp

nlcd = brick('data/rs/nlcd-L1.tif')
names(nlcd)  =  c("nlcd2001", "nlcd2011")


# The class names and colors for plotting
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)
# Hex codes of colors
classcolor  =  c("#5475A8", "#B50000", "#D2CDC0", "#38814E", "#AF963C", "#D1D182", "#FBF65D", "#C8E6F8")
# Now we ratify (RAT = "Raster Attribute Table") the ncld2011 (define RasterLayer as a categorical variable). This is helpful for plotting.
nlcd2011 = nlcd[[2]]
nlcd2011 = ratify(nlcd2011)
rat  =  levels(nlcd2011)[[1]]
#
rat$landcover = nlcdclass
levels(nlcd2011)  =  rat

Hicimos muchas cosas aqui. Da un paso atras y lee mas sobre ratify.

#####Nota No hay clase con valor 6.

TITULO: GENERAR SITIOS DE MUESTRA

Como discutimos en la clase, los datos de capacitacion y / o validacion pueden provenir de una variedad de fuentes. En este ejemplo, generaremos los sitios de muestra de capacitacion y validacion 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 

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

##
##   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 -cell- columna contiene numeros de celdas de -nlcd2011- muestra. -nlcd2011- La columna contiene los valores de clase (1-9). Dejaremos la -cell- columna mas tarde.

Aqui -nlcd- tiene valores enteros entre 1-9. A menudo encontrara que los nombres de clase se proporcionan como etiquetas de cadena (por ejemplo, agua, cultivo, vegetacion). Necesitara ‘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 clasificacion. Hay varios enfoques que podrian usarse para convertir estas clases a codigos enteros. Podemos hacer una funcion que reclasifique las cadenas de caracteres que representan las clases de cobertura del suelo en numeros enteros basados en los niveles de factores existentes.

-nlcd2011- Tracemos los sitios de entrenamiento sobre RasterLayer para visualizar la distribucion de las ubicaciones de muestreo.

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 mas avanzado (enrejado / enrejado) de objetos Raster *. Instale el paquete si no esta disponible para su maquina.

Aqui estan nuestros datos de Landsat.

landsat5 = stack('C:/Users/David-PC/Documents/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 seran las variables predictoras y los “valores de clase” -nlcd2011- seran 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)

TITULO: ENTRENAR AL CLASIFICADOR

Ahora entrenaremos el algoritmo de clasificaciOn 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)

En el diagrama de arbol de clasificacion, los valores de clase se imprimen en los nodos de hoja. Puede encontrar los nombres de cobertura de suelo correspondientes al uso del suelo en -classdf- data.frame.

Consulte -?rpart.control- para establecer diferentes parametros para construir el modelo.

Puede imprimir / trazar mas sobre el -cart- modelo creado en el ejemplo anterior. Por ejemplo, puede usar -plotcp(cart)- para aprender sobre la complejidad de costos ( -cp- argumento en -rpart-).

TITULO: CLASIFICAR

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

Importante Los nombres en el objeto Raster a clasificar deben coincidir exactamente con los esperados por el modelo. Este sera el caso si se utilizo el mismo objeto Raster (a traves de extracto) 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
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 :

## 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 clasificacion usando -rasterVis-. Ver establecera el -classnames- para 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")

Pregunta 1 : Trace “nlcd2011” y “pr2011” lado a lado y comente sobre la precision de la prediccion (p. Ej., Mezcla entre cultivos, pastos, pastizales y arbustos).

Es posible que deba seleccionar mas muestras y usar variables predictoras adicionales. La eleccion del clasificador tambien juega un papel importante.

library(gridExtra)
plot2011nlcd = levelplot(nlcd2011,
                       maxpixels = 1e6, 
                       xlab= NULL, ylab= NULL,
                       main='National Land Cover Database 2011',colorkey = FALSE,
                       scales=list(draw=FALSE), col.regions = classcolor,
                       labels = TRUE)


plot2011pr = levelplot(pr2011,
                       maxpixels = 1e6, 
                       xlab= NULL,nylab= NULL,
                       main='Prediccion 2011', colorkey = FALSE,
                       scales=list(draw=FALSE), col.regions = classcolor,
                       labels = TRUE)

grid.arrange(plot2011nlcd, plot2011pr, ncol=2)

NA
NA
NA

El modelo de prediccion es mucho mas detallado, lo que infiere que, probablemnete, es un modelo mucho mas acorde a la realidad de los datos

TITULO: EVALUACION DEL MODELO

Ahora, evaluemos la precision del modelo para tener una idea de la precision del mapa clasificado. Dos medidas ampliamente utilizadas en la teledeteccion son “precision general” y “kappa”. Puede realizar la evaluacion de precision utilizando las muestras independientes ( -validation2011-).

Para evaluar cualquier modelo, puede usar la validacion cruzada k-fold. En esta tecnica, los datos utilizados para ajustar el modelo se dividen en -k- grupos (tipicamente 5 grupos). A su vez, uno de los grupos se usara para la prueba del modelo, mientras que el resto de los datos se usara 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 confusion 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 confusion.

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

##                     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 sobre la clasificación erronea entre diferentes clases.

Pregunta 3 : ¿Se te ocurren formas de mejorar la precision?

Calcule la precision general y la estadistica “Kappa”.

Precision 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 precision del usuario

# Producer accuracy
PA = diag / colsums

# User accuracy
UA = diag / rowsums
outAcc = data.frame(producerAccuracy = PA, userAccuracy = UA)
outAcc

##                    producerAccuracy userAccuracy
## Water                     0.8663366        0.875
## Developed                 0.5421687        0.450
## Barren                    0.4270833        0.410
## Forest                    0.4076923        0.530
## Shrubland                 0.3566434        0.510
## Herbaceous                0.5291262        0.545
## Planted/Cultivated        0.4569536        0.345
## Wetlands                  0.4598540        0.315

Pregunta 4 : Realice la clasificacion utilizando clasificadores de bosque aleatorio del paquete “randomForest”

Pregunta 5 : Grafique los resultados del clasificador rpart y Random Forest uno al lado del otro.

Pregunta 6 (opcional) : repita los pasos para el año 2001 utilizando Random Forest . Use la imagen compuesta sin nubes -data/centralvalley-2001LE7.tif-. Estos son datos de Landsat 7 . Utilice como datos de referencia la Base de Datos Nacional de Cobertura de Tierras 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. Investigue el efecto del tamaño de la muestra en la clasificacion. Repita los pasos con diferentes subconjuntos, por ejemplo, un tamaño de muestra de 150, 100, 50 por clase, y compare los resultados. Use las mismas muestras de reserva para la evaluacion del modelo.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.