Paso a Paso: La Generación de Mapas de Expansión y Contracción de Cuerpos de Agua

Fase 1: Acceso y preprocesamiento de datos (continuación)

Este paso incluye el análisis de fechas de los eventos climatológicos Niño y Niña, la consulta y descarga de imágenes disponibles, y sus posteriores correcciones geométricas, radiométricas y atmosféricas.

Paso 1.2: Análisis de periodos

La selección de imágenes se realizó teniendo en cuenta dos factores:

1. series históricas de fenómenos Niño y Niña, suministrados por Corpoica, con base en la información oficial reportada por el Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM), y
2. los valores de PDSI (Palmer Drought Severity Index) calculados por Corpoica. Este índice clasifica los períodos de tiempo en seco y húmedo a partir del balance entre el suministro y la demanda de humedad (cf. figura abajo). De esta manera se pueden consultar las imágenes empleadas para el análisis multitemporal de la contracción y expansión de los cuerpos de agua y el análisis de sequías.

Figura : Períodos de análisis de los eventos Niño y Niña. Fuente Corpoica 2013.

Paso 2: Pre-procesamiento - Correcciones geométricas

La georeferenciación es el proceso de asignación de coordenadas a los píxeles en un sistema de referencia establecido. Una transformación geométrica cambia la posición de los píxeles de una imagen sin modificar los valores de los mismos. Este es un proceso importante para los análisis multitemporales puesto que asegura que los valores de reflectancia capturada por el sensor satelital en diferentes fechas, correspondan al mismo punto en la superficie (IGAC, 2004).

Este proceso utiliza diferentes algoritmos de asignación de coordenadas como: vecino más cercano, transformación bilineal y convolución cúbica. En este trabajo se decidió utilizar el primer algoritmo puesto que éste no altera los valores originales del píxel, considerando las siguientes características:

En este caso, utilizaremos la imagen Rapideye ya georeferenciada como base para corregir una imagen SPOT 5 del 2010.

• El primer paso consiste en salvar el nombre del proyecto en un archivo .prj, definiendo el método de modelamiento a seguir (File - New):

El programa se inicia con una ventana que requiere introducir los datos de un proyecto (File –New). En primer lugar se introduce un nombre de referencia para ese proyecto (p.ej. Rapideye_referencia), que puede extenderse con una descripción más cualitativa. Luego se indica qué tipo de corrección se va a abordar, si se trata de fotografía aérea, modelo orbital, polinómico, funciones spline, racionales o simplemente realizar mosaicos. Algunas de estas opciones requieren parámetros adicionales, como el tipo de cámara o de sensor y un MDT para la orthorectificación. En nuestro caso, indicaremos que se trata de un modelo orbital, para corregir imagen a imagen. El paso siguiente es indicar el método de corrección Toutin´s Model el cual emplearía los datos orbitales de nuestra imagen más los puntos de control GCPs extraídos de la imagen de referencia.

• Posteriormente indicamos los datos de la proyección de referencia:

Indicaremos: UTM, Elipsoide E012 (WGS 84), Zona 18, UTM Row S (hemisferio Norte). El tamaño del píxel de salida lo fijaremos a 10x10 metros. La proyección de los puntos de control, será la misma que la de referencia (Set GCP Projection base on Output projection).

• A continuación, dentro de Processing step señalar la opción para seleccionar los datos de entrada (Data Input):

En este caso leeremos los datos a través de la opción Read CD-ROM data. Seleccionaremos el formato de nuestra imagen SPOT 1-5 (DIMAP) seleccionamos el archivo de nuestra imagen original de toda la escena: (imagery), marcamos todos los canales. Esto nos creará una nueva imagen con metadatos en donde estarán todos los datos orbitales y puntos de control del satélite en un solo archivo. La guardamos en la misma carpeta de la imagen con un nombre apropiado: (ej: Imagen Spot_2010), incluyendo una descripción apropiada:

• A continuación abrimos en Focus los archivos creados: (Imagen Spot_2010_data) y el reporte con los metadatos incluidos:

• Posteriormente seguimos con el siguiente paso, seleccionando en processing step los puntos de control (GCP Collection):

Puede utilizarse como origen de esos puntos, un archivo con coordenadas GPS, introducirlos por teclado o por tableta, extraerlos a partir de una imagen ya georeferenciada o de vectores de referencia.

En esta práctica utilizaremos la imagen Rapideye (rapideye_2009.2) como referencia. Para comenzar el proceso, seleccionamos la imagen no corregida (Rapideye_referencia _data) con la opción Open Image (primer botón del panel, desde la izquierda).  Al abrirla, puede seleccionarse la combinación de bandas que nos resulte de mayor interés, con el tipo de realce que convenga. 

Una vez abierta, se pulsa el botón Collect GCPs Manually, lo que abre otro panel. En ese panel hay que indicar cuál es el origen de las coordenadas de referencia (Ground Control Source). Seleccionar la opción Geocoded image, y pulsar sobre la imagen Rapideye (rapideye_2009.2). Se recomienda utilizar la misma combinación de bandas que hayamos empleado con la imagen sin corregir.

Introducimos el modelo digital del terreno DEM, incluido en el archivo (DEM_30_Geosp). Si todo el proceso ha sido correcto, aparecerán las siguientes ventanas (Fig. 3):

Fig. 3. Ventanas de trabajo en Orthoengine, para la colección de GCP´s.

A la izquierda está la original sin corregir, a la derecha, la de referencia, y en el centro la ventana donde aparecen las coordenadas de los puntos de control que vayamos introduciendo. Se seleccionan de la forma siguiente. Observar puntos en la imagen de referencia que se ven con claridad en la no corregida. Seleccionar uno, con el nivel de zoom y de realce que se necesite, Cuando estemos satisfechos, pulsar Use Point consecutivamente en las dos imágenes. Aparecerán las coordenadas seleccionadas en la ventana central. Luego pulsar Accept en esa ventana para que el punto pase a considerarse válido. Cada punto puede tener una etiqueta propia si se requiere. A partir del 4º punto, el programa ya puede calcular la función de ajuste (en caso de ser de primer grado, si es de segundo, necesita 6 puntos) y facilita el error promedio global y el promedio en X e Y, así como el error de cada punto. Cuando está activada la opción Auto Locate, a partir del 4º punto al seleccionar un punto en la ventana de referencia, aparecerá su posición aproximada en la ventana de la imagen a corregir, o bien seleccionando Extract Elevation para cada punto que hayamos seleccionado.  Esto ayudará a localizar algunos rasgos del mapa en la imagen. Ojo, no aceptar necesariamente ese punto, pues el que estima el programa con el modelo que va construyendo según se le introducen puntos, pero eso no es correcto hasta que están todos los puntos seleccionados, así que habrá que mover el punto sugerido hasta el lugar más preciso que se pueda localizar. El proceso se continuará hasta tener 15-20 puntos y un error promedio inferior a 1 píxel. Si observamos puntos que dan mucho error, pueden corregirse haciendo doble clic sobre ellos, revisando si están las coordenadas bien posicionadas y, si no es el caso, moverlas y volver a pulsar Use Point en cada una de las dos ventanas. También pueden eliminarse si se descubre que la localización no es muy fiable. Hay que señalar también 3 puntos de verificación, de cara a estimar de modo más riguroso el error de la corrección. Basta cambiar el carácter del punto en la opción desplegable GCP, cambiándola a Check Point. Podemos ver el error residual de nuestra corrección, señalando la opción de Residual Report   en GCP TP Collection.

Al terminar el proceso, cerrar la ventana. En el menú principal de Orthoengine, convendrá salvar los puntos para un futuro uso, imprimirlos o modificarlos. Esto se obtiene con el menú Options Export GCP. Pueden exportarse a un segmento del archivo de imagen, de otra imagen o a un fichero ASCII. Salvarlo a un archivo de texto y como un segmento de la imagen a corregir, para el caso de que se quieran luego modificar los puntos o añadir otros. En caso de que queramos recuperar puntos de un archivo de texto, puede utilizarse la opción del menú GCP denominada Import GCP from File. Esta opción permite leer coordenadas de puntos de control en distintos formatos. Hay que indicarle al programa en qué orden están separadas las variables de cada punto: típicamente será identificador del punto, coordenada Píxel, Línea, X, Y (esto equivale al formato IPLXY del módulo Import GCP).

Una vez que el error sea adecuado, imprimir las ecuaciones de ajuste y los puntos finalmente seleccionados, tanto de control como de verificación. Ésto se obtiene en la pestaña Project Report, del menú Reports - Processing Steps. Ahí se indica la posibilidad de salvar a un archivo ASCII las ecuaciones (Geometric Model), los puntos (Point Information) y el área georeferenciada. Finalmente, ya estamos en condiciones de pasar a la última fase, el remuestreo de la imagen original a la nueva posición. Para ello seleccionamos Geometric Correction desde la ventana Processing Steps. Luego seleccionar el botón Schedule Ortho Generation, que abre una nueva ventana, donde se señalan los parámetros del remuestreo que vamos a realizar a partir de las ecuaciones previamente obtenidas:

En esta nueva ventana, se seleccionan las imágenes a corregir, que pasan pulsando la flecha de la lista Available Images a Images to Process. Seguidamente, se introducen los canales que quieren corregirse (las bandas que tendrá la nueva imagen, en principio todas las originales), el nombre de la imagen de salida (indicar un nombre descriptivo adecuado), y las coordenadas de la zona que se va a corregir. Introducir las que coincidan con la ventana de la imagen de referencia. El tamaño del píxel ya se indicó al iniciar el proceso, así como los parámetros de la proyección de salida. Finalmente se señala el método de interpolación (Nearest, Bilineal o Cubic, señalar el cúbico), y se pulsa el botón Correct Images. Este último proceso puede programarse para otra hora, si el volumen de las imágenes a corregir es muy grande, lo cual no es nuestro caso.

Ahora puede comprobarse la calidad de la corrección, abriendo la nueva imagen en Focus y superponiéndola sobre la original.

Paso 3.1: Correcciones atmosféricas/radiométricas

La radiación electromagnética que llega al sensor es almacenada como niveles digitales ND relacionados con las propiedades radiométricas de las coberturas, sin embargo estas no representan magnitudes físicas. Por ello es necesario convertir los ND a valores de radiancia Lλ, y así posteriormente obtener la reflectancia de superficie.

Sin embargo Lλ no depende únicamente de las propiedades de la cubierta, por lo que se requiere convertir de nuevo a una magnitud física como es la reflectancia (ᴩλ) de la cubierta en la banda λ captada por el sensor.

El uso de ᴩλ reduce las variaciones de la radiancia, permitiendo la comparación entre escenas de diferentes satélites adquiridas en fechas distintas. Para este trabajo emplearemos el modelo basado en el método simple de Chavez (1996) y en la ecuación propuesta por Pons & Sole (1994) para convertir los niveles digitales ND a niveles de radiancia y así obtener los niveles de reflectancia:

 (Chandler et al., 2009, p. 900)

Donde:

• ᴩλ = Reflectancia planetaria TOA o reflectancia en el techo de la atmósfera (Top of Atmosphere Reflectance) [sin unidades]
• θs= angulo zenital de sol. Los valores de θs se encuentran en los metadatos (archivos .MTL o .WO) de cada imagen tomada por cada sensor.
• π = constante matemática ~3.14150 [sin unidades]
• Lλ = radiancia espectral a la abertura del sensor [W/(m2 sr μm)]
• d=distancia de la tierra al sensor [unidades astronómicos]. EL valor de d, depende del día y año cuando se realizó la captura de la imagen.
• ESUNλ = medio irradiación solar exoatmosférica (mean exoatmospheric solar irradiance) [W/(m2 μm)]. ESUNλ se puede encontrar en resúmenes de trabajos previos (Chandler et al., 2009).

Arriba: Angulo zenital del sol en los metadatos de RapidEye.

Este procedimiento llamado también corrección atmosférica se desarrolló bajo el software Geomatica PCI 2013, debido a su sencillez para modelar y su capacidad de gestión. 

Paso 3.2: Generación del modelo de corrección

Las imágenes Spot 5 fueron sometidas a la evaluación de varios índices, realizando posteriormente un análisis de sensibilidad, descartando así qué índice se ajusta mejor para la extracción de superficies de agua. Para estos procesos se generó un modelo que fue llamado CAIN (Correcciones Atmosféricas e Índices). Para este proceso utilizaremos el módulo llamado: Modeler, el cual se utiliza para generar funciones adaptadas a las necesidades del usuario, mediante un lenguaje de programación en macros; llamados Pipelines o tuberías. PCI modeler está diseñado para ser una poderosa herramienta de análisis y además es fácil de usar, ya que posee una gran variedad de algoritmos que se puede interconectar mediante tubos o pipeline (Street, Hill, & Lb, 2003).

Dicho módulo se encuentra desplegado el panel superior de iconos Geomática 2013, tal y como se aprecia en la siguiente figura:

Figura: Menú de inicio de PCI y ubicación del módulo de programación Modeler en rojo.

Para poner en marcha este modelo, primero hay que describirlo y comprenderlo:

• Iniciaremos creando un nuevo modelo y salvándolo con un nombre apropiado desde el comando file:

• Posteriormente cargaremos cada una de las librerías de algoritmos necesarios para la implementación del modelo CAIN. Dichas librerías se encuentran alojadas en el menú Module Librarian, ubicado en la parte superior del menú principal:

Las principales librerías a utilizar son las siguientes: COMMENT, IMPORT, ARICONST, CRHN32R, MODEL, VIEWER Y EXPORT:

• Posteriormente procederemos a construir el modelo de corrección atmosférica comentado en la ecuación inicial, conectando mediante tuberías cada uno de los algoritmos importados desde Modeler (Dichas tuberías definen cada uno de los canales Raster o bandas a utilizar).
• Ahora empezaremos a construir el modelo paso a paso, en la ventana de Modeler  cargando una a una las librerías  a utilizar:

• Click derecho >> Common Modules

• Common Modules >> Import

Este módulo importa el conjunto de datos e imágenes necesarios que serán procesados posteriormente,tal y como se muestra en el siguiente diagrama (figura abajo), en donde aparece toda la estructura del modelo CAIN a través de la red de algoritmos conectados:

A continuación describiremos cada uno de los pasos que componen la estructura del modelo CAIN:

• En este primer paso, se cargará el archivo que contiene todas las bandas de la imagen SPOT.

Posteriormente separamos cada una de las bandas de la imagen:

• Click derecho >> Common Modules

• Common Modules >> Split

El módulo SPLIT divide los canales o bandas en vías de procesamiento separado.

El siguiente paso consiste en buscar la ruta de la imagen y daremos aceptar. Realmente para estas correcciones trabajaremos solo con las Bandas 2 (rojo) y 3 (infrarojo) de Spot 5 (para RapidEye son bandas 3 y 5), ya que son las necesarias para crear el índice espectral.

A continuación se une mediante una tubería el módulo IMPORT y el modulo SPLIT. Al unirse el SPLIT toma tonalidad verde. (Ésto quiere decir que el canal ya ha sido asignado correctamente).

Seguidamente agregamos el siguiente algoritmo de procesamiento:

• Menú>> View >> Module Librarian

• Find >>Escribir ARICONST>> Arrastra al espacio de trabajo

Definimos el valor de la constante y el operador a procesar con la capa raster:

• Find >>Escribir ARICONST>> Arrastra al espacio de trabajo:

Posteriormente implementar del mismo módulo ARICONST 5 veces.

En el primer ARICONST,  comenzaremos  a buscar los parámetros incluidos en los metadatos de la imagen.

Este caso pondremos el valor de ganancia para la banda a procesar, el cual aparece en los metadatos como GAIN, o valor de ganancia, este número hace referencia al coeficiente multiplicativo.

En el segundo módulo ARICONST  incluiremos el valor del Coeficiente Aditivo, también llamado BIAS, al igual que GAIN se obtiene en el metadato.

En el tercer módulo incluiremos el valor de (π), el cual quedará siempre constante:

Parta el siguiente módulo, observaremos también dentro de los metadatos el valor de SUNDIST como el factor que considera la distancia del Sol. Este valor será elevado al cuadrado antes de incorporarlo al modelo y se refiere a la distancia entre el Sol y la Tierra para el día de la toma.

Para el caso de nuestra imagen SPOT  será 1.031…

Posteriormente en el módulo marcado como WMO, incluiremos el valor que hace referencia a la energía exoatmosférica para cada una de las bandas del sensor: Y, que en nuestro caso para la banda 2 será de 1859.8

Posteriormente en el módulo marcado como COS (SZ)= SEN(AngElev), adicionaremos el valor que hace referencia al seno del ángulo de elevación, esta información se encuentra en el metadato

Hasta aquí nuestro espacio de trabajo de PCI Modeler. Tendría que aparecer de la siguiente manera:

A continuación incluiremos el módulo que permite salvar y exportar los archivos resultantes:

• Click derecho >> Commom Modules >> Export

Este módulo almacena el conjunto de datos después del procesamiento.

Daremos una ruta y nombre apropiado al nuevo archivo generado:

Posteriormente corremos el modelo generado:

De esta manera completaríamos todo el modelo de procesamiento CAIN para la correción atmosférica y obtención de imágenes en valores de reflectancia:

Paso 4: Procesamiento -Máscaras de Nubes

Con el fin de identificar y eliminar aquellas áreas con presencia de nubes dentro de las imágenes, procederemos a la creación de una máscara que nos permita excluir la información no deseada.

Este procedimiento lo desarrollaremos desde PCI Geomatica 2013. Como cualquier programa profesional de teledetección, FOCUS permite generar al usuario sus propias rutinas de procesamiento. Para ello dispone de dos opciones. Por un lado, un lenguaje propio de programación, similar en muchos aspectos al Basic, que se denomina EASI Modeling, accesible desde el menú Tools de la barra del programa. En este lenguaje se pueden escribir macros utilizando una sintaxis sencilla. En primer lugar es necesario indicar sobre qué archivo se va a aplicar el modelo (Input File).

• Menu >> Tools >> Easy Modelling

Éste es el conjunto estándar de operaciones aritméticas disponibles en las expresiones de la modelación:

Cada banda de ese archivo se identifica en el programa por un símbolo (%) seguido del número de canal. Para realizar una operación con los ND de los píxeles de una banda basta incluir en esa fórmula el identificador %X, siendo X la banda de interés. Por ejemplo, para dividir los ND de los píxeles almacenados en el canal 1 y el 2 de un archivo, salvando el resultado en el canal 5, se introduciría: %5=%1/%2. Conviene tener claro que el resultado debe enviarse a un canal vacío; de lo contrario, se estarían sobrescribiendo bandas con información de interés. Otro tema importante es que el resultado de la fórmula tiene que enviarse a un canal con la resolución apropiada. Por ejemplo, si dividimos dos bandas y enviamos el resultado a un canal de 8 bits, probablemente no veremos casi nada como resultado, ya que un canal de 8 bits no admite decimales y, por tanto, el resultado se redondea al entero más cercano perdiéndose la mayor parte de los valores resultantes. En consecuencia, antes de abordar cualquier proceso, habrá que crear bandas vacías sobre el archivo que queramos utilizar. Estas bandas serán de 8, 16 ó 32 bits dependiendo de lo que se quiera almacenar.

Paso 4.1: Procesamiento - Crear umbral para la discriminación de nubes

• Adicionar un nuevo canal  Raster de 32bits

Apoyándonos en la respuesta espectral de las coberturas nubosas en las diferentes longitudes de onda de las bandas de la imagen SPOT, esta característica permite crear un umbral para la exclusión de nubes y sombras de acuerdo a la geometría de iluminación. De esta manera desarrollamos un algoritmo para discriminación del áreas cubiertas por nubes bajo el módulo de programación “easy modeling” contenido en el programa PCI Geomatica 2013.

Dicha descripción se puede observar a continuación:

%6=0

if %3>6400 then

%6=255

endif

%6=0

if %1>3500 then

%6=255

else

if %2>3500 then

%6=255

endif

endif

Su visualización dentro del módulo “easy modeling”:

Nota: Cabe destacar que estos valores fueron empleados para esta escena SPOT en particular, y que sus resultados pueden variar o ser sensibles dependiendo de algunos factores como: geometría de iluminación, tipo de sensor, hora de la tomas entre otros.

• El resultado final visto como seudocolor:

Una comparación más detallada se puede apreciar en la siguiente figura :

Paso 4.2: Procesamiento - Creación de polígonos

Una vez creada la máscara en PCI Geomatica 2013, exportamos el canal a un formato raster .tiff con la finalidad de realizar un procedimiento de reclasificación en ArcMap.

• Exportar el canal creado como nuevo archivo .tiff desde el comando TRANSLATE de FOCUS:

A continuación abrimos el software ArcMap 10.2 e importamos la imagen de la máscara. Posteriormente accedemos a la caja de heramienta en donde se encuentran alojados los módulos de procesamiento y anális espacial “Toolbox”:

• Buscar en Toolbox, la herramienta Reclassify
• ArcToolboox >> Spatial Analyst Tools >> Reclass >> Reclassify

Input Raster: la imagen tiff de máscara de nubes. Creamos dos clases o categorías (Nubes y no nubes). Lo demás por defecto. Y damos OK.

Conversion de Raster a Polígono

Una vez reclasificada la imagen, procedemos a convertir nuestra capa Raster en una capa Vectorial de formato shape (.shp).

Para ello buscamos nuevamente en nuestro módulo de herramientas ToolBox el algoritmo Raster to Polygon:

• ArcToolboox >> Conversion Tools >> From Raster>> Raster to Polygon

• Input Raster: la imagen tiff reclasificada. Lo demás por defecto.OK

Paso 4.3: Visualización

Con el fin de mejorar el realce visual, procedemos a realizar un visualizado en categorías de un único valor de GRIDCODE.

• Click derecho en la capa vector recién creada >> Propierties…

• Layer Properties >> Symbology  

• Value Field: gridcode. Click en Add All Value. Aceptar

El resultado final nos muestra una salida gráfica con la visualización total de nubes para la escena seleccionada:

Nota: Para la extracción de polígonos, se tendrá en cuenta que para escalas de trabajo a 1:25.000, la unidad mínima cartografiable es de 1 Ha. De este modo, una vez extraigamos los polígonos debemos de asegurarnos que los polígonos son del tamaño adecuado. Para ello procederemos con las siguientes herramientas:

Reclasify

Raster to polygon

Eliminate, Majority

Paso 4.4: Eliminación de áreas mínimas

Una vez que tengamos el vectorial, crearemos en su tabla de atributos una columna para calcular el área en hectáreas, de este modo seleccionaremos los polígonos menos a 1 Ha, para eliminarlos con la siguientes herramientas:

• ArcTools Box >> Data Management Tools >> Generalization >> Eliminate
• ArcTools Box >> Spatial Analyst Tools >> Generalization >> Majority Filter