Análisis de Imágenes Satelitales Landsat 8 con Herramientas R
Alexander Vergara
Abril, 2020
1.Introducción
A partir de los diecisiete Objetivos de Desarrollo Sostenible establecidos en la Agenda 2030 de las Naciones Unidas, derivado de labores conjuntas entre los 193 Estados Miembros, y con los cuales se busca generar instrumentos que permitan alcanzar progresos a nivel social, económico, ambiental; conseguir naciones, ciudades y territorios sostenibles, capacitados para afrontar de manera apropiada el cambio climático, la reducción de recursos minerales y el uso adecuado de la tierra, a partir de mecanismos de reducción de la desigualdad en todas sus dimensiones. En este sentido y con el propósito de planificar adecuadamente el territorio se hace indispensable contar con herramientas y métodos que permitan generar información geográfica y ambiental con la cual se logren prevenir riesgos de desastres, localización de espacios adecuados para la agricultura, ganadería y asentamientos poblacionales, entre otros. Así, desde el advenimiento de diversas tecnologías que han apoyado el reconocimiento y mapeo de grandes extensiones de territorio, se cuenta con las que se ocupan del monitoreo y observación y de la tierra. A este segmento corresponden los satélites de la agencia espacial norteamericana, la cual ha venido desarrollando desde los años 60 uno de proyectos más fructíferos de teledetección espacial dedicado a la observación de los recursos terrestres. el primer satélite fue el ERTS Eart Resource Technology Satellite, denominado posteriormente Landsat (Chuvieco, 2010). Esta serie de satélites son hoy día los que cuentan con la mayor información, dada su vida útil. De manera paralela al desarrollo tecnológico satelital, ha sido menester desarrollar algoritmos y herramientas que, a partir del procesamiento de los datos provenientes de estos sistemas satelitales, permiten generar información valiosa en cuanto a la dinámica de la cobertura terrestre. Muchos de estos, conocidos como software de procesamiento de imágenes, han sido concebidos en enfoques académicos resultado de proyectos de investigación, como podrían citarse Idrisi e Ilwis, de la universidad de Clark e ITC de Holanda, respectivamente (ITC, 2018). Otros por su parte han sido desarrollo de instituciones militares, como Grass, de las Fuerzas Armadas de los Estados Unidos, al igual que aquellos que han sido fruto de labores de investigación de agencias espaciales de diversas latitudes, como es el caso de SNAP de la Agencia Espacial Europea (ESA) (ESA, 2018). De igual manera, empresas privadas de desarrollo de software específico se han dedicado a la generación de robustos programas con enfoque profesional, algunos de ellos denominados de coste y dentro de los cuales vale citar ENVI, PCIGeomatica, y ERDASAl entre otros. Por su parte, ciertas asociaciones de especialistas en áreas específicas de la teledetección, las geociencias y el computo entre otras, han aunado esfuerzos en desarrollar herramientas que satisfacen ciertas necesidades, con la posibilidad de involucrar profesionales que adhieren a la mejora continua de prestaciones de estos, los cuales tiene como característica ser de licencia de uso público, General Public License (GNU). En este segmento se hallan para diversas disciplinas, una buena cantidad de desarrollos informáticos, como es caso del software R, el cual ha tenido una importante evolucionan, desde sus denominaciones previas, como S, cuyo origen se enfocaba al análisis estadístico, hoy días gracias a la licencia GNU, se ha popularizado en áreas de la medicina, economía, bioinformática, modelado y estadística, entre otras tantas De esta última aplicación, es de donde probablemente se valen desarrollos en lo que respecta al análisis y procesamiento de imágenes de sensores remotos. Actualmente se han conseguido importantes rutinas específicas, denominadas paquetes, las cuales, a partir de funciones concretas, permiten a especialistas en teledetección y ciencias de la tierra, obtener información valiosa de la cobertura terrestre por medio del procesamiento de imágenes de múltiples sensores, así como de herramientas específicas desarrolladas para trabajo en R, z como rLandsat8 (GitHub, 2017).
2. Datos y métodos
2.1 Área de estudio
El área de estudio comprende un sector del Magdalena medio colombiano que cubre parcialmente los departamentos de Antioquia, Boyacá, Cundinamarca y Caldas que son irrigados por el Rio de la Magdalena. Son de importancia para el análisis los cascos urbanos de Puerto Berrío, Puerto Nare, Puerto Salgar, Puerto Triunfo y la Victoria. De igual manera lo hacen el área urbana de Puerto Boyacá, especialmente ciertos sectores de explotación de hidrocarburos de este municipio.
Figura 1. Localización del área de estudio sector del Magdalena medio (Colombia), en color rojo footprint de la escena Landsat 8 seleccionada.
2.2. Datos
Se obtuvo una escena Landsat 8 LC08_L1TP_008056_20200322_20200326_01_T1 del sitio https://glovis.usgs.gov/app compuesta por 11 bandas discriminadas por los sensores Operational Land Imager (OLI) el cual consigue información en las bandas denominadas: aerosol costero, azul, verde, rojo, NIR, SWIR1, SWIR2, Pancromático y Cirrus. Por su parte el sensor térmico infrarrojo Termal Infrared Sensor (TIRS) complementó los datos para las bandas 10 y 11 del Infrarrojo térmico. Los datos proporcionados por el U.S. Geological Survey (USGS) a través de glovis, comprendieron dos archivos de compresión. Como parte de estos, se contó con archivos LC08_L1TP_008056_20200322_20200326_01_T1_MTL de metadato de la escena, insumo requerido en cierto fragmento del procesamiento. De este se citan algunos aspectos, considerando su importancia. Cobertura de nubes en escena y en tierra: 13.4%, tipo de dato (nivel del producto): L1TP, corrección geométrica apoyada con puntos de control terrestre (GCP): 184, error medio cuadrático de corrección geométrica (RMSE): 9.26 m, ángulo de elevación solar en grados del momento de su adquisición: 61.85°, valores de píxel: 65535 (16bits). Datos como fecha de adquisición del archivo, al igual que los denominados path y row, se omiten en este aparte en tanto hacen parte del nombre de cada una de las bandas.
El procesamiento y análisis estadístico de las imágenes que componen la ensena Landsat 8, se realizó por medio del empleo del software R de código abierto. Tanto el cargue de las imágenes, como los procesos siguientes, se llevaron a cabo por medio de librerías específicas de manejo de información espacial, al igual que código definido por el software. Procesos concretos como cálculo de reflectancia por banda, demandaron desarrollos enfocados exclusivamente en la extracción y manejo de datos de imágenes Landsat 8 por medio de R, para el caso la librería rLandsat8. Se adquirió información espacial en formato .shp correspondiente al mapa de coberturas de la tierra a escala 1:100.000 del periodo 2010 - 2012, elaborado por el Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM), a partir de la metodología Corine Land Cover adaptada para Colombia. El sitio de descarga de la información corresponde al Sistema de Información Ambiental de Colombia (SIAC) Se plantea un método alternativo de evaluación cualitativa y cuantitativa de los datos contenidos en una imagen multiespectral Landsat8. Para ello se hace uso de herramientas de desarrollo libre en el software R para la selección de bandas que permitan la mejor discriminación de coberturas terrestre con propósitos de mapeo.
2.4. Métodos
Esta propuesta considera un método compuesto por las siguientes etapas: (i) análisis de datos fuente imagen, (ii) procesamiento de datos imagen y (iii) productos preliminares
2.5. Análisis de los datos imagen
Dentro del procedimiento de análisis preliminar de los datos, se procedió con una serie de códigos de solicitud de información y procedimiento en cada uno de los segmentos de escritura en R, y de los cuales se describen y resumen los pasos efectuados a continuación.
• Descripción de propiedades de las imágenes, y con lo cual por medio de la creación de objetos RasterLayer por banda individual, se accede a sus datos haciendo su llamado. Estos datos con: clase, dimensiones de la escena, resolución espacial, extensión sistema de referencia, nombre de cada banda y valores máximo y mínimo de niveles digitales contenidos por banda.
2.6. Procesamiento de datos imagen
• Conversión de valores de nivel digital DN a unidades absolutas de reflectancia, conforme al enfoque de determinación de coberturas. No se aplica en este estudio a las bandas térmicas. (Young et al., 2017), centrándose sobre los datos adquiridos por el sensor (OLI) al techo de la atmosfera. (Department of the Interior U.S. Geological Survey, 2016)
\[ ρλ = (Mρ * Qcal + Aρ)/θSE \]
Donde: Ρλ`= Valor de reflectancia planetaria. Mρ = Factor multiplicativo de escalado especifico por banda, en metadato (REFLECTANCE_MULT_BAND_x, (x) = numero de la banda). Aρ = Es el factor aditivo de escalado especifico por banda, en metadato (REFLECTANCE_ADD_BAND_x , (x) = numero de la banda). Qcal = (DN) de cada cada una de las bandas de la imagen. Θse = Angulo de elevación solar del centro de la escena, metadato (SUN_ELEVATION).
• Creación de pilas de bandas (RasterStack) a partir de los objetos RasterLayer de banda única llamados por código. Alternativamente se pueden crear las pilas de bandas (RasterStack) llamando directamente la banda desde la carpeta que las almacena.
• Ploteo de bandas únicas y composición de mapas. Este punto hace parte de las etapas del reconocimiento de las bandas que componen laescena. Es posible hacerlo de manera independiente por banda desde un RasterStack al igual que la generación de composiciones en color verdadero y falso color para evidenciar determinados comportamientos del paisaje en las imágenes.
• Renombre de bandas de la escena cuando se requiere extraer de la pila ciertas bandas de interés y asociar con los nombres asignados en función del segmento del espectro electromagnético en el que cada una de estas obtiene información.
• Creación y almacenamiento de un subconjunto espacial cuando se requiere extraer una porción espacial de la escena. Para ello se define la extensión con valores de coordenadas que esta debe tener acorde al sistema de referencia espacial que la escena tenga.
• La relación entre bandas puede llevarse a cabo por medio de matrices numéricas como la varianza y covarianza. Generar un esquema de matriz de dispersión, permite una mayor comprensión, por medio de la aproximación visual respecto a relaciones entre bandas.
• Extracción de valores de píxeles de áreas geográficas conforme a coberturas terrestres que componen la escena para la generación de las firmas espectrales. Para no depender únicamente del punto de vista de interpretación visual, se permite emplear información de referencia. Para el caso, fuente IDEAM del mapa de coberturas de la tierra.
• Generación de firmas espectrales específicas de la escena a partir de zonas de muestreo definidas en el proceso de extracción de valores de píxeles por cobertura y banda.
• Procedimientos aritméticos concretos entre bandas para la generación de Índices por medio de funciones de R
Índice normalizado de vegetación: \[ NDVI = NIR-Red/NIR+ Red\]
Índice normalizado de agua: \[ NDWI = NIR-SWIR/NIR+ SWIR\]
Índice normalizado de construcción: \[ NDBI = NIR/SWIR + NIR\]
2.7 Productos preliminares
• Clasificación no supervisada por medio del método más común de agrupamiento k-means. Este método busca delimitar clases espectrales presentes en la imagen. No demanda de conocimiento previo del área de estudio, su enfoque es más interpretativo que de obtención de resultados. (Chuvieco, 2010)
• Clasificación supervisada por medio de los algoritmos comunes CART y Random-Forest. Estos métodos parten de cierto conocimiento del área de interés a interpretar y delimitar en la imagen. La comprensión de la zona de estudio puede ser derivada de en trabajos previos de campo o referencias cartográficas para la delimitación de las categorías que conforman la leyenda. (Robert A. Schowengerdt, 1997)
• Extracción de valores de lugares puntuales desde los datos de referencia para la clasificación supervisada a partir de un clasificador. Los dos factores que influyen en el rendimiento del entrenamiento son tamaño y estrategia de muestreo del conjunto de datos. (Chandra P. Giri, 2016)
3. Resultados.
Los siguientes chunks (bloques de código), permiten hacer una exploración y reconocimiento la has bandas azul, verde, rojo, y NIR obteniendo los datos de: clase, dimensiones de la escena, resolución espacial, extensión sistema de referencia, nombre de cada banda y valores máximo y mínimo de niveles digitales contenidos por banda.
Estas imágenes cargadas corresponden a cortes realizados en software externo ArcGis. A manera de manera comparación, se cargaron las imágenes a las que se les aplicó la corrección de reflectancia TOA por medio de herramientas de rLandsat8. Notese de estas la extensión y los valores maximos y minomos de píxel. Parte de la identificación esta con el texto ‘TOA’ junto al nombre de la banda.
## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:magrittr':
##
## extract## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/57310/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/57310/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.4-1## class : RasterLayer
## dimensions : 2388, 3116, 7441008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 513234.9, 606714.9, 633089.7, 704729.7 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/57310/Documents/taller1/CTLC08_L1TP_008056_20200322_20200326_01_T1_B2.tif
## names : CTLC08_L1TP_008056_20200322_20200326_01_T1_B2
## values : 8640, 62412 (min, max)## class : RasterLayer
## dimensions : 2388, 3116, 7441008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 513234.9, 606714.9, 633089.7, 704729.7 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/57310/Documents/taller1/CTLC08_L1TP_008056_20200322_20200326_01_T1_B3.tif
## names : CTLC08_L1TP_008056_20200322_20200326_01_T1_B3
## values : 7233, 63208 (min, max)## class : RasterLayer
## dimensions : 2388, 3116, 7441008 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 513234.9, 606714.9, 633089.7, 704729.7 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/57310/Documents/taller1/CTLC08_L1TP_008056_20200322_20200326_01_T1_B4.tif
## names : CTLC08_L1TP_008056_20200322_20200326_01_T1_B4
## values : 6370, 65535 (min, max)## class : RasterLayer
## dimensions : 7731, 7571, 58531401 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 479985, 707115, 523185, 755115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/57310/Documents/taller1/LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B2.tif
## names : LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B2
## values : 0.04372, 1.14824 (min, max)## class : RasterLayer
## dimensions : 7731, 7571, 58531401 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 479985, 707115, 523185, 755115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/57310/Documents/taller1/LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B3.tif
## names : LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B3
## values : 0.02374, 1.16874 (min, max)Dentro del reconocimiento de las bandas, se esploraron sus estadísticas para la banda 2 sin corrección de reflectancia TOA
## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0## [1] 7441008## [1] 30 30## [1] 1## [1] TRUE## [1] 2388 3116 1Comparativamente se crearon pilas de bandas (RasterStack) en los dos bloques siguientes, llamando las bandas cargadas previo, o directamente desde la carpeta que las almacena. Sin importar el procedimiento, el objeto es proceso su requerimiento posterior
## class : RasterStack
## dimensions : 2388, 3116, 7441008, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 513234.9, 606714.9, 633089.7, 704729.7 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : CTLC08_L1TP_008056_20200322_20200326_01_T1_B2, CTLC08_L1TP_008056_20200322_20200326_01_T1_B3, CTLC08_L1TP_008056_20200322_20200326_01_T1_B4
## min values : 8640, 7233, 6370
## max values : 62412, 63208, 65535Dado el objeto del estudio, solo se cargaron las banas 1 a 7 con corrección de reflectancia TOA. Poseriormente se renombrarán acorde a la laongitud de onda de operación.
## [1] "taller1/LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B1.tif"
## [2] "taller1/LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B2.tif"
## [3] "taller1/LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B3.tif"
## [4] "taller1/LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B4.tif"
## [5] "taller1/LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B5.tif"
## [6] "taller1/LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B6.tif"
## [7] "taller1/LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B7.tif"## class : RasterStack
## dimensions : 7731, 7571, 58531401, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 479985, 707115, 523185, 755115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : LC08_L1TP//_T1_TOA_B1, LC08_L1TP//_T1_TOA_B2, LC08_L1TP//_T1_TOA_B3, LC08_L1TP//_T1_TOA_B4, LC08_L1TP//_T1_TOA_B5, LC08_L1TP//_T1_TOA_B6, LC08_L1TP//_T1_TOA_B7
## min values : 0.05988, 0.04372, 0.02374, 0.01274, 0.00468, -0.00720, -0.00294
## max values : 1.07446, 1.14824, 1.16874, 1.21070, 1.21070, 1.21070, 1.21070Sin entrar netamente en procesamiento digital de imágenes, se empleó el despliegue de las bandas 2 y 3 (azul y verde) respectivamente,sin corrección de reflectancia TOA. El comando strech provoca un estiramiento lineal en los valores del ráster Con lo cual se proporcionó un intervalo adecuado de salida visual.
Cuando no se realizan el reescalado de ND a reflectancia, y dadas las cualidades de las imágenes se presentan dificultades para su representación a color con un estiramiento (stretch) de tipo lineal, así que se optó en este caso por un estiramiento del histograma ‘hist’, como lo denomina R. Esta alternativa de visualización esta también Presente también en los software de pago especializados en procesamiento de imágenes.
Precisamente en los siguientes chunks, ilustran los comportamientos similares tanto en la iamgens de corte realizado en software de pago, pero con valores ND. El siguiente chunk presenta imagen de reescalado sa reflectancia en su extensión original. A simple vista el ojo humano independiente de las escalas de ploteo, no percibe diferencias entre ellas.
Una de las composiciones mayormente empleadas es la conocida como falso color o infrarrojo color. Esta es el resultado de aplicar cañones de color rojo verde y azul, sobre las bandas correspondientes al NIR, rojo y verde respectivamente. Esta composición permite la cartografía de extensiones de cobertura vegetal, cuerpos de agua y ciudades. Aspecto por el cual se cita en el presente, dada su aplicabilidad en estudios por medio de análisis visual.
El renombre de bandas de la escena se aplicó para extraer de la pila de bandas (stack) ciertas bandas de interés para asociar con los nombres asignados en función del segmento del espectro electromagnético en el que cada una de estas obtiene información. De igual manera en el bloque siguiente, se hazo una verificacin de las bandas extraidas.
## [1] 7## [1] 3## [1] 3De los siguientes bloques en adelante, se emplearon solamente las imágenes que pasaron por el cálculo de reflectancia TOA, dada su importancia en el análisis de coberturas. No obstante, este articulo no trata procesamiento de multi-temporalidad, aspecto que si seria de importancia para aplicar valores de reflectancia TOA que usar ND de las imágenes.
## [1] "LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B1"
## [2] "LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B2"
## [3] "LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B3"
## [4] "LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B4"
## [5] "LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B5"
## [6] "LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B6"
## [7] "LC08_L1TP_008056_20200322_20200326_01_T1_TOA_B7"## [1] "ultra.blue" "blue" "green" "red" "NIR"
## [6] "SWIR1" "SWIR2"Como sucede en muchas situaciones, el área de interés se halla en solo una porción de una escena Landsat, en tal sentido se requirió generar un subconjunto espacial o recorte de dicha área. Para ello se definió la extensión acorde a los valores de coordenadas del sistema de referencia espacial de la escena.
## class : Extent
## xmin : 479985
## xmax : 707115
## ymin : 523185
## ymax : 755115## class : Extent
## xmin : 514005
## xmax : 577995
## ymin : 634695
## ymax : 672105A manera de verificación del área recortada se consultan de esta, el número y nombres asignados a las bandas.
## [1] 7## [1] "ultra.blue" "blue" "green" "red" "NIR"
## [6] "SWIR1" "SWIR2"## [1] "red" "green" "blue"## class : RasterBrick
## dimensions : 1247, 2133, 2659851, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 514005, 577995, 634695, 672105 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : red, green, blue
## min values : 0.03700, 0.05640, 0.08738
## max values : 1.04330, 0.99180, 0.97966Al igual que se presentó en los bloques anteriores, se generó para el área específica de corte hecho en R, una composición RGB con el propósito de evidenciar coherencia con los procesos intermedios generados y un control del flujo de procesos.
Para facilidad de comprensión, al igual que sucede con las matrices citadas, una fácil lectura de la relación entre bandas puede hacerse en diagonal, bajando del verde al SWIR1 por los histogramas y verificando del diagrama de dispersión el comportamiento. Puntualmente se puede deducir que una combinación de bandas rojo- NIR, verde-NIR, incluso verde-SWIR, permiten la identificación de agua, por su parte bajos valores en verde y altos en NIR, que son bastante parecidos entre rojo y NIR, permitirían identificar cobertura vegetal. Esto representa, independiente de sus valores numéricos, que las menores correlaciones se dan entre la banda verde -SWIR, rojo- NIR, NIR SWIR y rojo- NIR. Por su parte NIR y SWIR, ofrecen algo de información para definir una clase de suelos. Claramente verde y rojo presentan del esquema la más lata correlación = 0.99, y como tal no se ofrece nada adicional de entre ellas.
Para las firmas espectrales se puede contar son diferentes procedimientos, para este ejercicio se tomó como referencia la información de cobertura de la tierra del Sistema de Información Ambiental de Colombia SIAC, cuyo propósito era el de contar con informan con grados de confianza aceptables, no solo para esta etapa, sino para algunas posteriores, como es la clasificación supervisada.
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\57310\Documents\Shp\Export_Output.shp", layer: "Export_Output"
## with 27 features
## It has 2 fields
## Integer64 fields read as strings: Id## class : SpatialPolygonsDataFrame
## features : 27
## extent : 522360.6, 572548.8, 635080.7, 667002.1 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 2
## names : Id, Cobertura
## min values : 0, Agua
## max values : 0, PastosDe los poligonos cargados, se generaron 200 muestras de puntos.Este proceso puede realizarse de manera alternativa con las mismas herramientas de R, no obstante, al realizarlo de esta manera se interioriza más el concepto de la variación espacial de los niveles digitales en las distintas bandas y su comportamiento en diagrama de curvas espectrales.
## class : SpatialPoints
## features : 193
## extent : 522432.9, 572313.1, 635379.6, 666882.9 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0A partir del mapa de cobertura de la tierra del SIAC, se obtuvieron las seis coberturas para la generación de las firmas espectrales. Este puede ser un buen ejercicio cuando no se tiene mayor adiestramiento en la interpretación visual, dado que permite generar polígonos de muestra de las coberturas respecto de la información vectorial complementaria con menores grados de incertidumbre en esta etapa.
## class : RasterBrick
## dimensions : 1247, 2133, 2659851, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 514005, 577995, 634695, 672105 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : ultra.blue, blue, green, red, NIR, SWIR1, SWIR2
## min values : 0.11226, 0.08738, 0.05640, 0.03700, 0.05862, 0.01412, 0.00588
## max values : 0.93314, 0.97966, 0.99180, 1.04330, 1.13464, 0.79704, 0.58688En el siguiente apartado se citan varios de los diferentes índices normalizados más empleados, ya que para vegetación la literatura cita que son alrededor de 66. Algunos de ellos son: Índice de vegetación de diferencia normalizada (NDVI), Índice de vegetación mejorado (EVI), Índice de vegetación de diferencia (DVI). Algunos de estos llegan a tener valores de correlación altos, sin embargo, pueden tener aplicados ambientales específicas. De estos se tienen Vegetación ajustada al suelo SAVI. Índice de vegetación atmosféricamente resistente ARVI. Índice monitoreo global ambiental. GEMI. (Chuvieco, 2010)
Para este ejercicio se aplicó uno de los usos que ofrecen los NDVI que corresponde con la discriminación de coberturas con reflectancias diferentes dos bandas a analizar Índice de vegetación de diferencia normalizada: \[ NDVI = NIR-Red/NIR+ Red\]
Una de las ventajas de los índices, el mejorar la discriminación de las coberturas vegetales y estimar ciertas variables físicas como se presenta en el siguiente ploteo.
Como se citó previo, dada la baja correlación que presentan las bandas rojo- NIR= 0.53 del diagrama de dispersión, se presenta una buena distinción de la vegetación. Son precisamente estas dos bandas las que permiten identificar la salud de una cobertura vegetal, dado que en la medida que estas bandas tengan mayor variación entre sus reflectividades, implicará un buen estado de la vegetación. En este punto es conveniente remitirse al esquema de curvas espectrales para tener una mayor compresión, con lo cual lo citado, se evidencia en las curvas especificas de coberturas vegetales (bosques y pastos) entre las bandas 3 y 5 para el esquema (rojo, NIR) respectivamente.
Adicionalmente se generaron los histogramas de frecuencia para cada uno de estos índices, permitiendo estimar para procesos posteriores, en que valores se estaba presentando mayor frecuencia, aspecto importante no solo para la definición de umbrales, sino base para la clasificación, en especial la no supervisada, en tanto ésta no cuenta con información auxiliar.
Un aspecto interesante del NDVI es su variación entre márgenes conocidos, facilitando la interpretación. Cuando la vegetación es densa este estará en 0.5 y 0.7 como se puede apreciar en el histograma, justo el pico valores se halla entre 0.5 y 0.6. No así sucede para un umbral crítico de cobertura vegetal próximo a 0.1.
Con fines comparativos, se presentan diferentes ploteos de índices de agua. Se citan en particular tres métodos: ZU (2006), MC Feeters (1996) y Gao (1996).
Métodos ZU: \[ NDWI = Green-SWIR1/Green+SWIR1\]
Método MC Feeters:
\[ NDWI = Green-NIR/Green+NIR\]
De los índices de agua entre los métodos ZU y Mc Feeters, se aprecian el histograma qué tanto las frecuencias de ZU como las de Mc Feeters son bastante diferentes respecto a la concentración de datos. Por su parte ZU tiene un Rango de valores más próximos al pico, razón por la cual se evidencia más los cuerpos de agua respecto de Mc Feeters, este por su parte sólo presenta un pico más alto con valores entre -0.45 y -0.5. En síntesis, y comparando estos dos se presenta una mayor definición de los cuerpos de agua en ZU, lo que no implica que éste sea el mejor de los métodos, pero para este caso, pudiese ser más apropiado para generar cartografía temática de manera semi-automática. Sin embargo, la mejor selección e bandas a partir de sus valores de correlación, seria Mc Feeters, quien emplea verde-NIR, cuyo valor es conforme a la matriz de correlación = 0.58 frente a 0.61 que propuesta por ZU.
Método Gao:
\[ NDWI = NIR-SWIR/NIR+SWIR\]
Por su parte GAO define muy bien el cuerpo de agua principal (Rio Magdalena), sin embargo, algunos afluentes a este se logran difuminar un poco. Esto se debe a que como se evidencia en su histograma, hay un mayor rango de valores de la diferencia en torno a lo que debería ser el pico. Se podría decir que, para aplicar este método, deben cumplirse ciertas características de dimensiones y valores de reflectancia de los cuerpos de agua para que sean bien definidas. Su valor de correlación es bueno, como se citó previo: Dada su concentración de valores permite la discriminación de agua, pero el pico inferior apreciable en la matriz de dispersión es más reducido en comparación con los otros métodos citados, lo que implica que sean menos datos destinados a la delimitación de cuerpos de agua.
El siguiene corresponde
Índice de normalizado de área construida:
\[ NDBI = SWIR-NIR/SWIR+NIR\]
Para este índice solo se aplicó la formula propuesta, no se genero un histograma. No obstante, se permite apreciar en los tonos azueles, algunos sectores donde se hallan asentamientos poblacionales, sin embargo y dada las características de este tipo de cobertura, las cuales podría decirse que casi comparten curva espectral, arenas (áreas de minería a cielo abierto) y zonas de contrición, como se presentaron en el diagrama de curvas espectrales. Allí se aprecian leves diferencias entre estas coberturas para las bandas que se aplican en este método, NIR-SWIR. Razón por la cual se presentan confusiones entre estas coberturas al plotear este índice.
El siguiente corresponde a la umbralización. Con lo cual se buscó complementar la comparabilidad de la cobertura del suelo, vegetación y cuerpos de agua, siendo capital y de mayor interés dados los valores de confianza respecto a la vegetación, determinar la densidad de cobertura como lo presenta el siguiente ploteo.
Se mapeó el área que corresponde al pico entre 0.55 y 0.60 del histograma NDVI para ilustrar respecto al area cubierta.
Este ploteo muestra la derivación de 3 clases que se obtuvieron del NDVI Rose (1974) y los valores de su histograma. Los rangos de dichas clases definidas corresponden a: 0.5, 0.55, 0.6, 0.65.
Con este PCA Análisis de Componentes Principales NIR - Red, se generó un diagrama igual que el de dispersión en la sección anterior. Los valores bajos para cada banda permiten la identificación de agua, es decir entre 0.05 y 0.1 para el tojo y el NIR. Por su parte el sector de valores prosimos a 0.4 del NIR y 0.7 del rojo, permitirían la identificación de la vegetación. Teniendo en cuenta que la reflectividad del suelo para las dos bandas tiene valores levemente diferentes, en el caso particular, no se definiría de la mejor manera, ya que se podría decir que éste se halla dentro de la diagonal donde se da la mayor correlación.
## Standard deviations (1, .., p=7):
## [1] 0.0972917835 0.0409560730 0.0380845156 0.0049514868 0.0033425812
## [6] 0.0014080904 0.0004753077
##
## Rotation (n x k) = (7 x 7):
## PC1 PC2 PC3 PC4 PC5
## ultra.blue 0.2890648 -0.30815738 0.17073335 0.340429597 0.47592256
## blue 0.3213378 -0.33255111 0.16907707 0.235397659 0.32689294
## green 0.3483064 -0.30766305 0.14684688 0.001214348 -0.21193257
## red 0.4009954 -0.37434392 0.06351208 -0.268148962 -0.62417122
## NIR 0.3908939 0.63770594 0.64785461 -0.137897586 0.02190789
## SWIR1 0.4873569 0.38930238 -0.56331269 0.506989285 -0.19058673
## SWIR2 0.3751722 0.03786348 -0.42382601 -0.693351579 0.44196191
## PC6 PC7
## ultra.blue -0.143150494 0.655319812
## blue -0.199236282 -0.745384832
## green 0.847078593 0.001749379
## red -0.467885442 0.120934656
## NIR -0.035756750 0.006575598
## SWIR1 -0.009728301 -0.012234139
## SWIR2 0.044505343 0.0116718572.7 Productos preliminares
Este apartado presenta el desarrollo de la clasificación no supervisada por medio del método de agrupamiento k-means. Con este método se buscó delimitar las clases espectrales presentes en la imagen de manera automática a partir de grupos homogéneos. Debido al enfoque de éste método, los resultados son meramente interpretativos. El método se aplicará al caso más sencillo, ya que se emplearon las bandas NIR y roja del proceso de NDVI. La literatura ilustra este procedimiento con gráficos bivariados, donde se hallan igual cantidad de puntos respecto de píxeles localizados de acuerdo con sus niveles digitales por banda.
Los siguientes corresponden a la extracción de valores de NDVI para ser clasificados acorde a la definición de clases.
## num [1:2659851] 0.59 0.598 0.612 0.592 0.541 ...Este bloque simplemente presenta los grupos derivados de la extracción e valores previos. Para el ejercicio se extrajeron 10 clases a partir de 500 iteraciones.
## List of 9
## $ cluster : int [1:2659851] 6 6 4 6 7 1 1 10 10 7 ...
## $ centers : num [1:10, 1] 0.4708 0.0817 0.2053 0.6248 0.3195 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:10] "1" "2" "3" "4" ...
## .. ..$ : NULL
## $ totss : num 31650
## $ withinss : num [1:10] 56.9 72.2 54.7 56 54.1 ...
## $ tot.withinss: num 582
## $ betweenss : num 31068
## $ size : int [1:10] 268622 50144 48336 453620 62980 509486 486926 241326 125255 413156
## $ iter : int 328
## $ ifault : NULL
## - attr(*, "class")= chr "kmeans"## class : RasterLayer
## dimensions : 1247, 2133, 2659851 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 514005, 577995, 634695, 672105 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 1, 10 (min, max)El siguiente ploteo es presenta las 10 clases definidas. Como se citó anteriormente los resultados son meramente interpretativos. Estas categorías en si no se podrían comparar con las clases informaciones que se buscan. De esta imagen es poco el significado temático que acorde a las categorías espectrales se puede inferir.
Para una clasificación supervisada, se debería contar con algunos tipos de cobertura del suelo de los datos espaciales de referencia obtenidos del catálogo de mapas y con los cuales se apoye la identificación de sitios específicos de entrenamiento en el área de estudio. Se debería aplicar el principio de leer coberturas con cierta homogeneidad, no obstante, las condiciones de ciertas de ellas, como sucede con los cuerpos de agua, los cuales en ciertos sectores presentan bastante heterogeneidad, aspecto que se debe a su dinámica. En el caso particular se hubiese podido emplear la misma información de referencia que se usó para la generación de las curvas espectrales. Sin embargo el alcance del presente se dio hasta la clasificación no supervisada.
4. Discusión
En el caso de la clasificación no supervisada, teniendo en cuenta que el (los) algoritmo(s) planteados en la literatura se basan en el agrupamiento de píxeles con rasgos espectrales similares, el potencial de este se halla en la discriminación previa que se hizo por medio de los índices espectrales. Como se plantea en la práctica con los NDVI, NDWI y NDBI, entre otros tantos.
Teniendo en cuenta que K-means genera un grupo etiquetado a partir de los valores de los índices determinados, lo cual implica un ejercicio de ida y vuelta de ráster a etiquetado, y de nuevo a ráster, el cual y, como es obvio, debe tener la misma extensión que el índice tomado de base, para el caso particular, el NDVI. Así, generar un gráfico del ráster derivado del proceso de K-means junto, o sobre una composición RGB o el mismo NDVI, puede ser de gran utilidad cuando no se dispone de datos o información de apoyo para la clasificación. En tal sentido será menester de parte del intérprete, asociar las clases definidas a las coberturas que ofrecen los ploteos restantes, con el propósito de que, a partir de apoyo visual y experticia, pueda redefinir, bien sea por medio de aumento o agrupación de las clases máscaras específicas.
Se presenta el típico problema de traslado de una clase a otra por cuanto en el ejercicio práctico, el analista opera en ciertos sectores de la imagen donde a partir de interpretación visual, considera ciertas clases. Esto está limitado a su vez, por el dominio que se tenga de los elementos que conforman la interpretación visual, por citar algunos de los más relevantes: forma, tamaño, textura, color, tono, localización, entre otros. Pero se pierde el potencial de las bandas susceptibles en las longitudes de onda fuera del rango visible. Precisamente donde es posible identificar fenómenos naturales particulares.
Conclusiones
El desarrollo de las actividades vistas en cada uno de los puntos, permitió una mejor comprensión del comportamiento de las coberturas terrestres en torno al espectro electromagnético.
Lo anterior permite evaluar no solo las ventajas y desventajas de la clasificación unianbada, sino la posibilidad de combinaciones que realmente permitan una mejor discriminación de las cubiertas terrestres; el entender el comportamiento de la firmas espectrales conlleva a una mejor toma de decisiones en cuento a que banda emplear acorde al objeto de estudio
El ejercicio de generar las curvas espectrales para cada cobertura implica en muchos casos, hacer uso de los elementos de interpretación visual para una correcta definición de éstas. Cuando en este procedimiento se involucran datos externos de apoyo, como pueden ser interpretaciones o mapas de cobertura elaborados por instituciones, se afianzan aspectos cómo son el entender que ciertas coberturas pueden llegar a tener un rango de valores más amplio, evitando de esta manera, muestras de pequeños grupos de píxeles con valores próximos entre sí. Este aspecto permite que, en actividades siguientes como son la determinación de las zonas de entrenamiento para procesos de clasificación supervisada, no se limiten al uso de espacios centrales y homogéneos dentro de cada cobertura, lo cual aproxima a que las máscaras generadas se acerquen a una mejor exactitud temática, aun cuando se exija más al algoritmo, aspecto en el que radica el desarrollo de estos mismos. La asertividad en cuanto a la escogencia de las composiciones resulta vital para efectos del apoyo de la percepción remota en proyectos de ingeniería, planeación y ordenamiento del territorio.
La posibilidad de asesorar proyectos agroindustriales como es el caso de la agricultura de precisión, requieren del entendimiento y apropiado manejo no solo de los aspectos teóricos en cuanto a la percepción remota, sino también en cuanto al empleo de los insumos acorde a la serie de variables que se puedan presentar; económicos, logísticos, técnicos, entre otros, que en resumen, pueden brindar resultados que beneficien o afecten una área determinada con el proceso de información geoespacial de estas cualidades.
Si bien el alcance del presente no llegó hasta la clasificación supervisada, es menester para el empleo adecuado de información apoyo la temporalidad y la visualización de esta respecto de los datos a procesar, esto permite tener mayor confianza de los lugares a entrenar sin equivocaciones al ser llamados ciegamente, generando al software inconvenientes con polígonos sobre coberturas bastante diferentes, Ej. Sectores de cambios de curso de cuerpos de agua.
Bibliografía
Colombia, S. d. (16 de Abril de 2020). SIAC: Catálogo de mapas. Obtenido de SIAC Web site: http://www.siac.gov.co/catalogo-de-mapas
Chandra P. Giri. (2016). Remote Sensing of Land Use and Land Cover. In Remote Sensing of Land Use and Land Cover. https://doi.org/10.1201/b11964
CODAZZI, I. G. (20 de Diciembre de 2013). UN-Spider home. Obtenido de UN-Spider Web site: http://www.un-spider.org/sites/default/files/LDCM-L8.R1.pdf
Chuvieco, E. (2010). Teledeteción Ambiental, La observacion de la tierra desde el espacio. Barcelona: Editorial Ariel.
Department of the Interior U.S. Geological Survey. (2016). Landsat 8 Data Users Handbook. Nasa, 8(June), 97. Retrieved from https://landsat.usgs.gov/documents/Landsat8DataUsersHandbook.pdf
ESA. (1 de Enero de 2018). Science tool box explotation plataform. Obtenido de SNAP: https://step.esa.int/main/toolboxes/snap/
GitHub. (21 de Diciembre de 2017). Terradue, rLandsat8. Obtenido de Paquete R para procesar datos del USGS Landsat 8: https://github.com/Terradue/rLandsat8
ITC. (19 de Marzo de 2018). ILWIS: software de detección remota y SIG. Obtenido de ILWIS Web Site: https://www.itc.nl/ilwis/
Robert A. Schowengerdt. (1997). Remote Sensing Models and Methods for Image Processing. In Academic Press.
USA-CERL. (14 de Diciembre de 2019). Grass Gis home. Obtenido de grass web site: https://grass.osgeo.org/
USGS. (18 de Abril de 2020). Data Management and Information Distribution (DMID). Obtenido de USGS Web site: https://lta.cr.usgs.gov/EEHelp/ee_help
Young, N. E., Anderson, R. S., Chignell, S. M., Vorster, A. G., Lawrence, R., & Evangelista, P. H. (2017). A survival guide to Landsat preprocessing. Ecology, 98(4), 920–932. https://doi.org/10.1002/ecy.1730