Parte2_taller1
Julieth Perez
22/4/2020
DOS VERSIONES DEL MAPA DE LAS COBERTURAS DE LA TIERRA EN UN ÁREA DE PRUEBA SELECCIONADA, UBICADA EN EL DEPARTAMENTO DE MAGDALENA; GENERADOS MEDIANTE CLASIFICACION NO SUPERVISADA Y CLASIFICACION SUPERVISADA, UTILIZANDO LOS ALGORITMOS K-MEANS, Y RANDOM FOREST RESPECTIVAMENTE.-PARTE 2
El presente documento, se encuentra dividido en dos partes, en la primera parte, se hace una introducción al mismo, se presentan informacion de los datos; y los métodos a utilizar.En la segunda parte, se presenta el procesamiento de la clasificación no supervisada en el área de estudio, mediante el uso del algoritmo K-means; y el procesamiento realizado para generar la clasificación supervisada con el algoritmo Random Forest; los resultados, conclusiones y recomendaciones.
La primera parte de este taller se puede observar en: https://rpubs.com/juperezve/603314.
CONTENIDOS PARTE 2
6.RESULTADOS Y DISCUSIÓN DE RESULTADOS, 7.CONCLUSIONES, 8. RECOMENDACIONES.
6. RESULTADOS Y DISCUSIÓN DE RESULTADOS
6.1 CLASIFICACIÓN NO SUPERVISADA
6.1.1 NDVI
ndvi <- (Corte4[['NIR']] - Corte4[['red']]) / (Corte4[['NIR']] + Corte4[['red']])
plot(ndvi, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')
A continuación el histograma del NDVI generado (Figura 9); el mayor numero de datos se encuentran en los valores entre 0,55 y 0,65; con el histograma se puede observar que existe pocas áreas urbanizadas o completamente desprovistas de vegetación, dado que es bajo el numero de datos con valores de 0 o negativos.
hist(ndvi,
main = "DISTRIBUCIÓN VALORES NDVI",
xlab = "NDVI",
ylab= "Frecuencia",
col = "wheat",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-0.01,1, 0.05), labels = seq(-0.01,1, 0.05))
6.1.2 k MEANS
Se le asigna al algoritmo un numero de clases, y este agrupa los valores del ndvi, que se realizó en la primera parte, en 11 clusters.
tablavalores <- getValues(ndvi)
km <-kmeans(na.omit(tablavalores), 11)
rNA <- setValues(raster(ndvi), 0)#crear un raster en blanco con valores por defecto 0
for(i in 1:nlayers(ndvi)){ rNA[is.na(ndvi[[i]])] <- 1}# asignar valor de 1 a valores NA del raster rbrick
rNA <- getValues(rNA)#convertir rNA en un vector de valores
tablavalores <- as.data.frame(tablavalores)#convertir tablavalores en un data frame
tablavalores$clase[rNA==0] <- km$cluster #Si rNA es 0, asignar el valor del cluster k-mean
tablavalores$clase[rNA==1] <- NA#Si rNA es 1, asignar NA
clases <- raster(ndvi)
clases <- setValues(clases, tablavalores$clase)
#asignar los valores del raster clases a partir de tablavalores$claseA continuación se encuentra el resultado de la clasificación obtenido, en donde K-means genera los 11 cluster a partir de los resultados del NDVI (Figura 10).
#asignar los valores del raster clases a partir de tablavalores$clase
mycolor <- c("#fef65b","#ff0000", "#daa520","#0000ff","#0000ff","#00ff00","#cbbeb5",
"#c3ff5b", "#ff7373", "#00ff00", "#808080")
plot(clases, main = 'CLASIFICACIÓN NO SUPERVISADA', col = mycolor )
El resultado obtenido, puede ser una primera aproximación a la clasificación en un área donde no se cuenta con información o muestreo en campo de la zona. En la Figura 11, se muestra la comparación del NDV, con la clasificación generada.
mycolor <- c("#fef65b","#ff0000", "#daa520","#0000ff","#0000ff","#00ff00","#cbbeb5",
"#c3ff5b", "#ff7373", "#00ff00", "#808080")
par(mfrow = c(1,2))
plot(ndvi, col = rev(terrain.colors(11)), main = 'Landsat-NDVI')
plot(clases, main = 'CLASIF NO SUPERVISADA', col = mycolor )
La Figura 12, muestra un acercamiento en donde se sobrepone la capa de coberturas de la tierra escala 25.000, con la que se cuenta sobre toda el área de estudio; a la clasificación no supervisada generada por K-means; lo que da un indicio cualitativo del ajuste de la misma.
knitr::include_graphics("C:/informe1/ACERCAMIENTOIMGNOSUP.jpg")
6.2 CLASIFICACIÓN SUPERVISADA
6.2.1 RANDOM FOREST
En primer lugar se cargan las áreas de entrenamiento, las cuales se encuentran como poligonos en un shapefile; la función “head”, sirve para ver el encabezado de la capa, en donde se observa la clase, y el numero asignado a cada clase; en este caso se cuenta con muestras para 22 clases; esto es dado que el mapa de coberturas utilizado para el entrenamiento cuenta con estas 22 clases.
entrenamiento<-readOGR(".","Entrenamiento")## OGR data source with driver: ESRI Shapefile
## Source: "C:\informe1", layer: "Entrenamiento"
## with 487 features
## It has 4 fields
## Integer64 fields read as strings: CODIGOhead(entrenamiento@data)## CODIGO ID_NUM ID_STRING CLASE
## 0 231 6 6 PASTOS LIMPIOS
## 1 314 14 14 BOSQUE GALERIA
## 2 314 14 14 BOSQUE GALERIA
## 3 314 14 14 BOSQUE GALERIA
## 4 314 14 14 BOSQUE GALERIA
## 5 314 14 14 BOSQUE GALERIA
A continuación, se muestran las clases con las que cuenta el shapefile de áreas de entrenamiento y el numero de cada clase, asignado para generar la clasificación.
knitr::include_graphics("C:/informe1/NUMCLASES.jpg")
En la Figura 14, se muestra un acercamiento de la imagen L8, compocisión 543, con la capa de coberturas de la tierra escala 25.000, realizada a partir de clasificación visual de imagenes Sentinel 2; que se utiliza para para obtener las muestras de entrenamiento de la clasificaión supervisada.
knitr::include_graphics("C:/informe1/ACERCAMIENTOIMG1.jpg")
Una vez se cuenta con el raster sobre el cual se va a generar la clasificación, y las áreas de entrenamiento; estos ultimos se rasterizan, en la Figura 15 se muestra la distribución de los mismos en el área de estudio.
#rasterizar poligonos
clase <- rasterize(entrenamiento, Corte4, 'ID_NUM')
spplot(clase, scales=list(draw=T), main="POLÍGONOS ENTRENAMIENTO")
A continuación se pueden observar para cada clase, el numero de muestras o pixeles tomados para entrenar el algoritmo de clasificacion.Se observa que la clase con mayor numero de muestras es la 6, que corresponde a Pastos limpios; esto en atención a que es la cobertura con mayor porcentaje de área en la zona de estudio.
tab <- table(values(clase))
print(tab)##
## 1 2 3 4 5 6 7 8 9 10 11 12 13
## 453 242 78 337 523 39321 22718 2496 2800 213 264 141 84
## 14 15 16 17 18 19 20 21 22
## 5386 386 2696 6276 1404 537 589 51 104sum(tab)## [1] 87099
Los poligonos de entrenamiento se rasterizan, y con estos se extrae el corte de los mismos en el raster stack que contiene las 7 bandas a utilizar; es decir, se utiliza para el modelo, las áreas del raster que coinciden con las de entrenamiento y a estas se les asigna una clase (Figura 16).
#aplicar mascara sobre corte4 volver rastes el vector
mascara<- mask(Corte4, clase)
entrenamiento_r <- addLayer(mascara, clase)# entrenamiento_r es el raster de los poligonos
plot(entrenamiento_r)
Los atributos del raster de entrenamiento finalmente generado se muestran a continuación.
entrenamiento_r## class : RasterStack
## dimensions : 1334, 1154, 1539436, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 586275, 620895, 1027965, 1067985 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : blue, green, red, NIR, SWIR1, SWIR2, layer
## min values : 130, 378, 277, 1209, 855, 519, 1
## max values : 5326, 5679, 5805, 7151, 6602, 5144, 22- Validación cruzada
El siguiente código muestra como entrenar el algoritmo y a su vez generar la matriz de confusion; para este caso se entreno para generar 100 árboles de decisión.<div
names(entrenamiento_r)<- c("blue", "green", "red", "NIR", "SWIR1","SWIR2","clase")
tablavalores <- getValues(entrenamiento_r)
tablavalores <- na.omit(tablavalores)
tablavalores <- as.data.frame(tablavalores)
tablavalores$clase <- factor(tablavalores$clase, levels = c(1:22))
modelRF <- randomForest(x=tablavalores[ ,c(1:6)], y=tablavalores$clase, ntree=100, importance = TRUE, do.trace=FALSE)
colnames(modelRF$confusion) <- c("Nubes", "Tejido Urbano Discontinuo", "Vivienda Rural Dispersa", "Vias", "Cultivo permanente arboreo", "Pastos limpios", "Pastos arbolados", "Pastos Enmalezados", "Pastos inundables", "Mosaico de cultivos", "Mosaico de pastos y cultivos", "Mosaico de pastos y espacios naturales", "Mosaico cultivos y espacios nat", "Bosque de galeria", "Plantacion forestal", "Arbustal", "Vegetacion secundaria", "Tierras desnudas y degradadas", "Zonas quemadas", "Zonas pantanosas", "lagunas", "cuerpos agua art", "error de clase")
rownames(modelRF$confusion) <- c("Nubes", "Tejido Urbano Discontinuo", "Vivienda Rural Dispersa", "Vias", "Cultivo permanente arboreo", "Pastos limpios", "Pastos arbolados", "Pastos Enmalezados", "Pastos inundables", "Mosaico de cultivos", "Mosaico de pastos y cultivos", "Mosaico de pastos y espacios naturales", "Mosaico cultivos y espacios nat", "Bosque de galeria", "Plantacion forestal", "Arbustal", "Vegetacion secundaria", "Tierras desnudas y degradadas", "Zonas quemadas", "Zonas pantanosas", "lagunas", "cuerpos agua art")
#para mejorar la clasificación deben incluirse mas poligonos que representen adecuadamente las clases seleccionadas.#Exportar matriz de confusion a excel
#write.table(x = modelRF$confusion, file = "matriz2.txt", sep = ",", row.names = FALSE, col.names = TRUE)- Matriz de confusión
A continuación se muestra la matriz de confusión generada por el modelo, con el error obtenido para cada clase; los valores que se acercan a 1, son aquellos para los cuales el error es muy alto, es deir que el nivel de confiabilidad de la clasificación para esta clase es bajo; esto sucede principalmente con clases de poca área en las cuales no fue suficiente el muestreo realizado.
knitr::include_graphics("C:/informe1/matriz_final.jpg")
A continuación resultados del modelo. La figura muestra los parametros establecidos, y el resultado que arroja sobre la tasa de error estimado, que corresponde al 35,1%; el valor de exactitud del modelo corresponde a 100% menos la tasa de error estimada; por lo cual la exactitud del mismo es 64.9%
knitr::include_graphics("C:/informe1/res_entren_randomforest.jpg")
se muestran abajo, indicadores que posee el modelo;en primer lugar (Mean Decrease Accuracy), el índice de disminución media de la precisión, que muestra que bandas que mejoran la precisión de la clasificación; y, el indicador de Gini o disminución media de Gini (Mean Decrease Gini), que permite identificar que bandas presentan mayor homogeneidad en la clasificación.
varImpPlot(modelRF)
Con el modelo generado y matriz de confusión,el algoritmo genera la clasificación, que contiene las 22 clases establecidas inicialmente.
predLC <- predict(Corte4, model=modelRF, filename="predicccion.tif", na.rm=TRUE, overwrite=TRUE)
save(modelRF, file="prediction_modelRF.RData")
plot(predLC, main = 'CLASIF SUPERVISADA', col = mycolor )
- Índice Kappa
A continuación el calculo de precisión general del modelo y clasificación
#Precisión gral
#Numero muestras
#OA Overal Accuracy
n <- (87099)
# numero de muestras correctamente clasificadas
sum_diag <-(56619)
#Precisión
#OA <- sum(diag) / n
OA=sum_diag/n
OA## [1] 0.6500534A continuación, el resultado de la exactitud de usuario; es decir, la exactitud por cada clase, en donde se puede observar cuales clases arrojaron valores bajos de exactitud. Las clases que tienen mayores valores son aquellas en donde se tomaron mayor número de muestras de entrenamiento.
knitr::include_graphics("C:/informe1/exactitud usuario.jpg")
Finalmente, se calcula el Kappa global para medir la exactitud temática de la capa generada; el cual arroja un valor de 0,48; lo que indica que es necesario o bajar el nivel de la clasifcación, es decir hacer menos clases; o aumentar el numero de muestras por clase para mejorar la precisión; especialmente en las clases que arrojaron valores bajos en la matriz de confusión y exactitud de usario.
#kappa global
rowsums <- apply(modelRF$confusion, 1, sum)
p <- rowsums / n
colsums <- apply(modelRF$confusion, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
kappa## [1] 0.4848982A continuación, se observa la comparación de las dos clasificaciones realizadas, en donde se observa que es preferible realizar una supervisada con un buen numero de zonas de entrenamiento para mejorar la calidad de la clasificación.
par(mfrow = c(1,2))
plot(predLC, col =mycolor, main = 'CLASIF SUPERVISADA')
plot(clases, main = 'CLASIF NO SUPERVISADA', col = mycolor )
El mapa de coberturas generado con el algoritmo Random Forest, se muestra a continuación.
knitr::include_graphics("C:/informe1/mapa_clasfinal.jpg")
7. CONCLUSIONES
Con respecto a la clasifcación no supervisada se observa que el resultado genera una primera aproximación, alrealizar la separación en clusters; sin embargo, si no se conocen las caracteristicas del área, asignar un numero de clases para generar el algoritmo, genera una mayor incertidumbre en el resultado generado, pues esto se realiza a partir de supocisiones con poco fundamenteo.Sin embargo, esta información es importante para etapas exploratorias, o para definir un muestreo de campo permitiendo que este tipo de actividades se realicen con un mejor sustento.
Respecto a la clasificación supervisada realizada con el algoritmo Random Forest; se observa que al tener sitios de entrenamiento y conocimiento de las coberturas del área de estudio; se puede generar un mapa temático más ajustado la realidad; y de la misma forma, un producto en el que es posible hacer evaluación de su exactitud temática, y en caso de no encontrar un resultado aceptable, modificar los parametros y variables del modelo hasta tanto se generen resultados admisibles.
Así mismo, es importante resaltar que el resultado obtenido, depende también del nivel de detalle requerido; en este caso, la capa utilizada para generar los sitios de entrenamientos, es escala 1:25.000 , y fue generada con interpretación visual en imagenes de resolución espacial de 10 metros (Sentinel2A) con una escala de dibujo de 1:5.000; por lo cual, el resultado de la clasificacón automatica generada con el uso de imagenes landsat 8, de 30 metros de resolución espacial, utilizdas para generar salidad 1:50.000 o escalas mayores; obtendra un resultado, que aunque es bueno o aceptable; tendra sus limitaciones y es necesario tener en cuenta el alcance del mapa generado y el uso de la información que contiene.
En el caso del resultado obtenido de la validación del modelo y clasificacion supervisada realizada se observa que la tasa de error estimado (OOB), corresponde a 35,1 mientras que el valor de exactitud (1- oob) corresponde a 64.1.
Respecto al Mean Decrease Accuracy, (índice de disminución media de la precisión), se observa que las bandas que mejoran la precisión, corresponden al NIR, blue y SWIR1; encontrandose la relación alta entre el uso de los infrarojos en la clasificación.
El indicador de disminución media de Gini (Mean Decrease Gini), permite identificar que las bandas SWIR1, SWIR2 y NIR son aquellas que presentan mayor homogeneidad de los datos utilizados en la clasificación; lo cual concuerda con lo descrito en el parráfo anterior.
Las coberturas que arrojaron una mayor exactitud de clase fueron Nubeas, debido principalmente a que tiene una respuesta espectral muy caracteristica y fácil de diferencias; Pastos limpios, Arbustal, Pastos inundables, y vegetación secundaria.
El Kappa global corresponde a 0,48; el cual si bien no es bajo; no es el resultado esperado a un nivel de escala de reconocimiento o semidetallado por o cual se deben evaluar los resultados y tomar acciones que permitan mejorar el modelo para así obtener un resultados más satisfactorio.
8.RECOMENDACIONES
Hacer un PCA previo a la realización de la clasificación y realizar nuevamente el modelo con lo obtenido de los tres primeros componentes arrojados, podria eliminar sobrantes de información y ruido que altera la homogeneidad de la muestra.
Tomar mayor numero de muestras de entrenamiento en las clases que arrojaron una baja exactitud de usuario y un error alto en la matriz de confusión; es el caso del Tejido urbano discontinuo, Vivienda rural dispersa, vías, lagunas y los mosaicos.
Así mismo, debido a que algunas de las áreas existentes de las coberturas que obtuvieron una baja exactitud de usuario; es baja y no permite ampliar el numero de muestras de entrenamiento es preferible realizar en estas clases, una interpretación visual y generar una mascara para extraer estas zonas del raster a clasificar y una vez hecha la clasificación, incluir esta información. Esto con el fin de disminuir la desviación debido a datos atipicos.