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 1

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.

CONTENIDOS PARTE 1

1.INTRODUCCIÓN 2.OBJETIVO 3.ÁREA ESTUDIO 4.DATOS 5.MÉTODOS.

1. INTRODUCCIÓN

Este es un cuaderno de R Markdown en el cual se realiza el procesamiento de un subset espacial de una escena LANDSAT 8, con PATH ROW 8 54, con identificador del dato LC08_008054_20190811, para generar dos versiones del mapa de coberturas la tierra, en un área de prueba ubicada en el departamento del Magdalena; los cuales se generan a partir de clasificacion supervisada y no supervisada; esto, con el fin de realizar un primer acercamiento del uso y analisis de los algoritmos K-means y random forest; respectivamente, utilizados al realizar las clasificaciones.

Dado que corresponde a un ejercicio academico, al mismo tiempo se realiza la exploración, y caracterización, de los datos obtenidos del subset de la escena Landsat 8 utilizada para el presente ejercicio.

Para la validación de la exactitud temática; en el caso de la clasificación supervisada se utiliza un archivo formato shapefile; que contiene la clasificacion de coberturas de la tierra, del área de estudio a escala 1:25.000, la cual se encuentra de acuerdo a la Leyenda Nacional de Coberturas de la Tierra, Metodología Corine Land Cover, adaptada para Colombia (IDEAM, 2010), en el nivel de clase 3, con algunas modificacines en las clases; la cual fue generada a partir de interpretación visual de imagenes Sentinel 2, por el Grupo de Interpretación de Coberturas de la Tierra escala semidetallada y detallada; del Instituto Geográfico Agustin Codazzi.

2. OBJETIVO

Esta publicación tiene como fin, generar dos versiones del mapa de coberturas de la tierra, en un área de prueba (area de estudio seleccionada), ubicada, en el departamento del Magdalena; los cuales se realizan a partir de clasificación supervisada y no supervisada en una imagen LANDSAT 8.

Así mismo, ilustrar al lector acerca de como realizar la exploración, caracterización del conjunto de datos de la imagen y uno de los usos obtenidos en la escena seleccionada.

3. ÁREA DE ESTUDIO

La zona de estudio, cuenta con un área de 72.291,43 ha, abarcando parte de seis municipios del departamento del Magdalena. A continuación se puede observar la ubicación de coordenada de referencia de la zona, en el municipio de Pijiño del Carmen (lng=-74.15, lat=9.59).

m <- leaflet()
m <- addTiles(m)
m <- addMarkers(m, lng=-74.15, lat=9.59, popup="COORDENADA REFERENCIA ÁREA ESTUDIO")
m

Los municipios que abarca el área corresponden a Santa Ana, Pijiño del Carmen, San Zénon, San Sebastian de Buena Vista, Guamal y El Banco; en el departamento del Magdalena. A continuación se observa el shapefile de coberturas utilizado, que define el área de estudio.

cobertura3<- shapefile("C:/informe1/shape_info1")
plot(cobertura3, axes=F, main="Area de Estudio")

En la Figura 1, a continuación, se observan los municipios que se encuentran en el área de estudio seleccionada, dentro del deparatmento del Magdalena, la codificación de municipios de acuerdo a DANE corresponde a (Tabla1 ):

MUNICIPIO	COD MPIO
El Banco	245
Guamal	318
Pijiño del Carmen	545
San Sebastian de Buena Vista	692
San Zenon	703
Santa Ana	707.

A continuación. municipios del área de estudio (Figura 1).

knitr::include_graphics("C:/informe1/mpios.jpg")

4.DATOS

4.1.LANDSAT 8

La mision Landsat 8 fue lanzada en febrero de 2013; el satelite se encuentra equipado con dos sensores: OLI (Operational Land Imager) que cuenta con un radiometro de barrido multicanal; y TIRS (Thermal Infrared Sensor), el cual es un radiometro infrarrojo de dos canales.

Las caracteristicas de las imagenes se listan a continuación; esta información se extractó textual del “Data Users HandBook - Landsat 8 (L8)” (USGS, 2019), Version 5.0.

L8 adquiere habitualmentealrededor de 700 escenas por día, Los datos OLI y TIRS para cada escena, se fusionan para crear un solo producto que contiene los datos de ambos sensores. Los datos de ambos sensores se corrigen radiométricamente y se registran conjuntamente en una proyección cartográfica, con correcciones por desplazamiento del terreno que dan como resultado una imagen digital ortorectificada estándar llamada Producto de Nivel 1 con corrección de terreno (L1T).

El sensor OLI recopila datos de imagen para 9 bandas espectrales de onda corta en una franja de 190 km con una resolución espacial de 30 metros (m) para todas las bandas, excepto la banda Pancromatica que es de 15 m; El sensor TIRS, recopila datos de imagen para dos bandas térmicas con una resolución espacial de 100 m en una franja de 190 km. En la Figura 2 a continuación se observa la descripción de las bandas de los dos sensores.

knitr::include_graphics("C:/informe1/Bandas_L8.jpg")

Igualmente, a continuación (Figura 3), se observa el ancho de las bandas previamente descritas, de estos dos sensores, en el espectro electromagnetico.

knitr::include_graphics("C:/informe1/bandas_dibujo.jpg")

Fuente imagen: http://mappingandco.com/blog/disfrutando-del-landsat-8-1-parte-especificaciones-tecnicas/

Productos Nivel 1

Los datos de Nivel 1 disponibles para los usuarios son una imagen corregida radiométricamente y geométricamente. Se utilizan entradas tanto de los sensores como de la nave espacial, así como GCP y DEM. El resultado es un producto rectificado geométricamente libre de distorsiones relacionadas con el sensor y Tierra. La imagen también se corrige radiométricamente para eliminar las diferencias relativas del detector, el sesgo de corriente oscura y algunos artefactos. La imagen de Nivel 1 se presenta en unidades de DN, que se pueden reescalar fácilmente a radiancia espectral o reflectancia TOA.

Un producto completo de Nivel 1 consta de 13 archivos, incluidas las imágenes de 11 bandas, un archivo de metadatos específicos del producto y una imagen de Evaluación de calidad (QA). Los archivos de imagen son todas las imágenes GeoTIFF de 16 bits. Las bandas OLI son las bandas 1-9. Las bandas TIRS se designan como Bandas 10 y 11.

La imagen de control de calidad es una máscara de 16 bits, que marca nubes, datos de relleno y algunos tipos de cobertura del suelo.El archivo de metadatos de Nivel 1 (MTL) contiene parámetros de identificación para la escena, junto con la extensión espacial de la escena y los parámetros de procesamiento utilizados para generar el producto de Nivel 1. Este archivo es un archivo de texto legible.

Exploración de datos

En primer lugar se carga el raster seleccionado, que como se describio antes corresponde a una escena Landsat8;con PATH ROW 8 54, con identificador del dato LC08_008054_20190811; del 8 de noviembre de 2019; surface reflectance.

En la lista se observa que las bandas de la imagen, se encuentran en nivel L1T, con corrección geometrica y radiometrica presentandose en “SURFACE REFLECTANCE” se observan únicamente las bandas 1 a 7, debido a que, el Servicio Geologíco de los Estados Unidos - USGS, entrega este tipo de imagenes bajo demanda; y la reflectancia en superficie, no se encuentra disponible para las bandas termicas ni pancromatica.

Es de resaltar que SR (surface reflectance), mide la fracción de radiación solar reflejada desde la superficie de la tierra, hasta el sensor; por lo que esta imagen presenta la reflectancia de la superficie, la cual se genera a partir de entradas de nivel 1 que cumplen con la restricción de ángulo solar de <76 grados e incluyen las entradas de otros datos auxiliares necesarios para hacer la corrección y generar el producto en este nivel de procesamiento.

Es importante mencionar, que el USGS, cuenta con una interfaz bajo demanda de productos landsat (EROS Science Processing Architecture - ESPA). En la página de la ESPA (https://espa.cr.usgs.gov/); se pueden solicitar productos a granel o pedidos masivos en diferentes niveles de procesamiento.Para desplegar el raster seleccionado en primer lugar, se cargan las bandas a utilizar, que como se describio anteriormente corresponden a las bandas 2 a la 7.

Posteriormente se realiza un stack con las bandas; la imagen se encuentra en el sistema UTM ZONA 18N; cuenta con una resolución espacial de 30 metros.A continuación en la Figura 4, se observa la escena L8 completa, desplegadada en infrarojo (utilizando las bandas b5, b4, b3), combinación especifica para resaltar vegetación.

landsatRGB <- stack(b5, b4, b3)
plotRGB(landsatRGB, axes = TRUE, stretch = "hist", main = "ESCENA L8 FALSO COLOR COMPOSICIÓN 5,4,3")

A continuación se observa el summary de las bandas desplegadas, con sus valores mínimos, maximos, quantiles y medidas de tendencia central.

landsat_prop <- stack(b2, b3, b4, b5, b6, b7)
summary(landsat_prop)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.17% of all cells)

##         LC08_L1TP_008054_20190811_20190820_01_T1_sr_band2
## Min.                                             -1333.00
## 1st Qu.                                            206.00
## Median                                             405.00
## 3rd Qu.                                           1029.75
## Max.                                             10564.00
## NA's                                          17932851.00
##         LC08_L1TP_008054_20190811_20190820_01_T1_sr_band3
## Min.                                                 -234
## 1st Qu.                                               489
## Median                                                774
## 3rd Qu.                                              1207
## Max.                                                11348
## NA's                                             17932851
##         LC08_L1TP_008054_20190811_20190820_01_T1_sr_band4
## Min.                                                 -313
## 1st Qu.                                               335
## Median                                                682
## 3rd Qu.                                              1117
## Max.                                                11901
## NA's                                             17932851
##         LC08_L1TP_008054_20190811_20190820_01_T1_sr_band5
## Min.                                                  291
## 1st Qu.                                              2998
## Median                                               3474
## 3rd Qu.                                              3955
## Max.                                                12678
## NA's                                             17932851
##         LC08_L1TP_008054_20190811_20190820_01_T1_sr_band6
## Min.                                                  106
## 1st Qu.                                              1491
## Median                                               1939
## 3rd Qu.                                              2446
## Max.                                                 8568
## NA's                                             17932851
##         LC08_L1TP_008054_20190811_20190820_01_T1_sr_band7
## Min.                                                   59
## 1st Qu.                                               656
## Median                                               1001
## 3rd Qu.                                              1403
## Max.                                                 6467
## NA's                                             17932851

A continuación, se observan algunas de las principales caracteristicas, del metadato (LC08_L1TP_008054_20190811_20190820_01_T1_MTL) de la escena L8 utilizada, se observa que es colección 1, nivel L1T, tomada el 8 de agosto de 2019 y procesada el 20 de agosto del mismo año; entre otros.

Metadato Imagen
LANDSAT_SCENE_ID = “LC80080542019223LGN00”
LANDSAT_PRODUCT_ID = “LC08_L1TP_008054_20190811_20190820_01_T1”
COLLECTION_NUMBER = 01
COLLECTION_CATEGORY = “T1”
WRS_PATH = 8
WRS_ROW = 54
MAP_PROJECTION = “UTM”
DATUM = “WGS84”
ELLIPSOID = “WGS84”
UTM_ZONE = 18
GRID_CELL_SIZE_PANCHROMATIC = 15.00
GRID_CELL_SIZE_REFLECTIVE = 30.00
GRID_CELL_SIZE_THERMAL = 30.00
ORIENTATION = “NORTH_UP”
RESAMPLING_OPTION = “CUBIC_CONVOLUTION”.

Seguidamente, se genera un stack con las bandas 2 a 7 (se exluye la banda costera) y se genera el recorte del raster stack con el shapefile del área de estudio.

landsat_unido <- stack(b2, b3, b4, b5, b6, b7)
#Se genera el recorte con la función crop a continuación 
cobertura_corte<- shapefile("C:/informe1/shape_info2")
corte = crop(landsat_unido,cobertura_corte)
corte2 = mask(corte, cobertura_corte)
#Exporto el recorte generado 
#writeRaster(corte2, filename="prueba4_recorte", format="GTiff", overwrite=TRUE)

El raster stack del recorte realizado, corresponde al área en donde se realizara el mapa de coberturas mediante los dos métodos descritos (clasificación no supervisada y clasificación suervisada). En este raster hago un list con los nombres de las bandas en su interior que corresponden a las bandas 2 al 7.

Corte4 <- stack('C:/informe1/prueba4_recorte.tif')
names(Corte4) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')

A continuación, en la Figura 7, se observa el recorte realizado, en compocision 543 (infrarojo), resaltando la vegetación.

plotRGB(Corte4, r="NIR", g="red", b="green", axes=TRUE, stretch="lin", main="RECORTE COMPOCISIÓN INFRAROJO (543)")

Los valores resumidos de las bandas, una vez realizado el recorte se encuentran a continuación; se observa que los maximos, minimos y valores de tendencia central cambian con respecto a la escena Completa; esto es debido a que al hacer el corte, el numero de datos disminuye, por lo cual si solo se utiliza un subset de la imagen para el procesamiento; es necesario calcular las estadisticas despues de haber realizado el corte, para obtener una información mas ajustada a la realidad.

summary(Corte4)

## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (6.5% of all cells)

##           blue  green    red    NIR  SWIR1  SWIR2
## Min.       100    405    282    855    701    474
## 1st Qu.    340    785    648   3422   2266   1182
## Median     463    922    838   3652   2570   1406
## 3rd Qu.    818   1053    998   3903   2882   1657
## Max.      5180   5491   5606   6790   6123   4809
## NA's    811514 811514 811514 811514 811514 811514

5. MÉTODOS

5.1 Clasificación digital de imagenes

Corresponde a la clasificación realizada a partir de la respuesta espectral de pixeles en imagenes de satelite; esto se realiza para pasar de una variable continua a una discreta o categorica, consistente en clases. En este caso, a clases de cobertura de la tierra.

Las técnicas de clasificación por respuesta espectral se pueden separar en supervisadas y no supervisadas.

Las clasificaciónes realizadas a continuación, se hacen con base en la guía “Clasificación sin supervision”, (GFC) en https://rspatial.org/raster/rs/4-unsupclassification.html; la guía “Herramientas avanzadas de investigación - Análisis espacial con R” (Ballari, 2015) en https://rpubs.com/daniballari/analisis_espacial_indice, y la a guía “Clasificación supervisada”, (GFC) en https://rspatial.org/raster/rs/5-supclassification.html.

5.2 Clasificación digital no supervisada

En este tipo de clasificación es el propio algoritmo el qué decide qué características han de cumplir los elementos que pertenecerán a una determinada clase; y clasifica de acuerdo a los valores espectrales de los pixeles desconocidos en una imagen; posteriormente los agrupa en una serie de clases que comparten afinidad estadística y espacial de los valores espectrales.

5.3 Clasificación digital supervisada

La clasificación supervisada es aquella en que se le señala al algoritmo, cuantas clases generar, y se toman muestras o áreas de entrenamiento en la imagen, para cada clase de forma que se toma una muestra representativa de cada una;con las áreas de entrenamiento, el algoritmo agrupa los pixeles por su similitud espectral.

5.4 NDVI

El Indice de Vegetación dDiferencial Normalizado NDVI, separa la vegetación verde del suelo (Rouse et al; 1974) la formula de este corresponde a: NVDI=NIR-RED/NIR+RED;sus valores se encuentran entre -1 a 1; siendo 1 área cubierta de vegetación con alta densidad y -1 un área desprovista de vegetación. La clasificación se realiza con el raster generado del NDVI, dado que mejora la discriminación entre el suelo y la vegetación, reduciendo el efecto del relieve en la clasificación(Rouse et al; 1974) En la Figura 8, se muestra el resultado obtenido del NDVI, a simple vista, se puede notar que existen pocas áreas desprovistas de vegetación.

5.5. Algoritmos a utilizar

5.5.1 K-means

El algoritmo a utilizar para realizar la clasificación no supervisada es el K-means; definido por J, MacQueen, en 1967; permite realizar agrupamientos o clusterización en un conjunto de datos, a partir de la respuesta espectral de los mismos; en donde cada grupo esta representado por los promedios de los valores espectrales.

5.5.2 Random Forest

Es un algoritmo, de técnica Machine Learning, O aprendizaje automatico, el cual utiliza árboles de decisión como clasificadores base (Brieman, 2001). Para los agrupamientos Random Forests selecciona al azar un subconjunto de los atributos para luego volver a seleccionar el mejor corte entre estos. Posteriormente el proceso se repite en cada uno de los árboles; los árboles son usados en el resultado final a partir del promedio de los resultados de cada uno de los árboles (Brieman, 2001).

Validación cruzada

El algoritmo utilizado, genera una validación cruzada interna, en donde se calcua un un valor de exactitud y error de clasificación. Esta toma el raster que contiene las áreas que se extrajeron del raster stack a clasificar,a partir del recorte con los poligonos de entrenamiento que contienen las clases asignadas a estos polgígonos; y, utiliza un porcentaje para hacer la clasificación, vs otro porcentaje para realizar la validación del modelo.

Debido a que el algoritmo en cada árbol de decision generado,“aprende automaticamente”, este toma la decision de cuantas muestras tomar para el entrenamiento y cuantas para la validación del modelo.

Seguidamente se realiza el entrenamiento del algoritmo y se genera el modelo. En este caso, el algoritmo genera coo resultado una matriz de confusión, en donde se observa la validación; la diagonal de esta matriz muestra para cada clase, el numero de muestras que se encontraron bien clasificadas.

Indice Kappa

Mide el grado de concordancia de las evaluaciones nominales u ordinales realizadas por múltiples evaluadores cuando se evalúan las mismas muestras.Los valores de kappa varían de –1 a +1. Mientras más alto sea el valor de kappa, más fuerte será la concordancia.

Parte1_taller1

Julieth Perez

22/4/2020