Atmospheric correction/es

Vea también la wiki de LANDSAT.

Propósitos de este tutorial:

• Transformar una imagen Landsat ETM+ de los valores de Números Digitales (DN) a valores de radiancia.
• Realizar una corrección atmosférica de las imágenes Landsat usando el módulo i.atcorr en GRASS GIS y convertir la radiancia a reflectancia.

Introducción

Los pasos estándar de pre procesamiento de imágenes satelitales son:

1. Convertir los Números Digitales (DN) a Radiancia por encima de la atmósfera (ToAR).
   • para las imágenes Landsat, aplicar la ecuación de abajo o usar el módulo i.landsat.toar.
   • Los Números Digitales --en realidad valores de pixel-- son el resultado de la cantidad cuantificada de energía obervada y medida en el sensor.
   • El módulo i.landsat.toar soporta todos los sensores Landsat incluyendo MSS, TM y ETM+.
2. Para convertir Radiancia por encima de la atmósfera a Reflectancia por encima de la copa (ToCR) → use el módulo i.atcorr.
   • La conversión a Reflectancias intenta medir y eliminar los efectos atmosféricos, proceso denominado corrección atmosférica.
   • La corrección atmosférica es un tema muy común de discusión en Sensores Remotos. Sin embargo, dependiendo de la calidad de los datos, puede ser un proceso complejo que no siempre da el un buen resultado, en cuando a la exactitud de los siguientes pasos (ej. segmentación o clasificación de las imágenes).

Recuerde revisar si las zonas con agua están representadas por valores >0, dado que la reflectancia siempre es positiva. Si encuentra un valor negativo en la reflectacia tiene un problemilla. Esto significa que las ecuaciones usadas para la transformación, no corresponden a la imagen procesada.

Un resumen de las opciones que tiene para "corregir" las reflectancias espectrales:

 +------------------------------------------------------------------------------+
 |                                                                              |
 | Números digitales                                                            |
 |        |                                                                     |
 |  +-----v-----+                                                               |
 |  |  i.*.toar | ---> Reflectancia                                             |
 |  +-----+-----+     (no corregida)~~~(métodos DOS)--+                         |
 |        |                 +                         |     Reflectancia        |
 |    (bandera -r)     (bandera -r)                   +-->   "Corregida"        |
 |        |                 |                         |                         |
 |        v           +-----v----+                    |                         |
 |     Radiancia -----> i.atcorr +--------------------+                         |
 |                    +----------+                                              |
 |                                                                              |
 +------------------------------------------------------------------------------+

Requisitos

• GRASS 6.4.0 o superior
• Conjunto de datos de muestra de North Carolina (NC) (localización): https://grass.osgeo.org/download/sample-data/

Luego de bajar los datos, debe descomprimirlos y ponerlos en su base de datos SIG. Cuando inicie GRASS GIS, seleccione esta base de datos, abra la lcalización ´nc_spm_08´ como Localización, y ´PERMANENT´ como Directorio de mapas.

Calculando valores de radiancia

El conjunto de datos de muestra de NC, entre otras cosas contiene imágenes de Landsat ETM+ del 24 de mayo de 2002. Cada pixel de esta imagen contiene un número digital (DN) o valor de gris. Para poder hacer los cálculos con las imágenes satelitales, o comparar valores entre diferentes sensores, es necesario covertir estos valores a radiancias o reflectancias. Las fórmulas para hacer esta conversión para las imágenes de Landsat se describen a continuación (o use el módulo i.landsat.toar).

La conversión de Número digitales a radiancias por encima de la atmósfera se hace de la siguiente manera:

Failed to parse (syntax error): {\displaystyle L\lambda = \frac{L_{MAX}λ - L_{MIN}λ}{QCAL_{MAX} - QCAL_{MIN}} * (QCAL - QCAL_{MIN}) + L_{MIN}\lambda}

donde:
• Failed to parse (syntax error): {\displaystyle Lλ} - Radiancia espectral en la apertura del sensor en watts/metro cuadrado * ster * μm), esta es la radiancia aparente registrada por el sensor;
• $QCAL$ - el valor del pixel calibrado en $DN$;
• Failed to parse (syntax error): {\displaystyle LMINλ} - la radiancia espectral que es escalada a $QCALMIN$ en watts/(metro cuadrado * ster * μm);
• Failed to parse (syntax error): {\displaystyle LMAXλ} - la radiancia espectral que es escalada a $QCALMAX$ en watts/(metro cuadrado * ster * μm);
• $QCALMINM$ - el valor mínimo cuantificado de pixel (correspondiente a Failed to parse (syntax error): {\displaystyle LMINλ} ) en $DN$ ;
• $QCALMAX$ - el valor mínimo cuantificado de pixel (correspondiente a Failed to parse (syntax error): {\displaystyle LMAXλ} ) en $DN=255$.

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle LMINλ} y Failed to parse (syntax error): {\displaystyle LMAXλ} son las radiancias relacionadas a los valores máximos y mínimos de los números digitales $DN$, y son reportados en el archivo de metadatos de la imagen, o en la tabla 1. La ganancia superior e inferior también se reporta en el archivo de metadatos de cada imagen. El valor mínimo del $DN$ ($QCALMIN$) es 1 para imágenes Landsat ETM+ (1), y el máximo ($QCALMAX$) es de 255. $QCAL$ es el $DN$ para cada pixel independiente en la imagen Landsat.

Accesando los metadatos:

r.info lsat7_2002_xx

Por debajo de ´coments´, se da el valor míimo y máxmio de la radiancia para cada banda (LMAX y LMIN).

Los usuarios más experimentados puede usar esto para obtener una salida legible (ej. para el canal 1):

CHAN=1
r.info lsat7_2002_${CHAN}0 -h | tr '\n' ' '  | sed 's+ ++g' | tr ':' '\n' | grep "LMIN_BAND${CHAN}\|LMAX_BAND${CHAN}"

La conversión a radiancia (aquí se hace solamente para la banda 1, para las otras bandas, los números en itálicas tienen que ser reemplazados por valores adecuados):

g.region rast=lsat7_2002_10 -p
r.mapcalc "lsat7_2002_10_rad=((''191.6'' - (''-6.2'')) / (255.0 - 1.0)) * (lsat7_2002_10 - 1.0) + (''-6.2'')"

Para cálculos más rápidos, usar un MDE de enteros:

 r.mapcalc "elev_int = round(elevation)"

Encontrar la elevación promedio para iniciar los cálculos (necesario para el archivo 'icnd' que se muestra abajo):

r.univar elev_int

Estimando el tiempo de sobrevuelo

Este valor es necesario para el archivo de control. Para más detalles, vea la wikipedia en Scientific decimal time. En caso de que este parámetro no esté reportado en el archivo de metadatos, hay dos maneras de obtenerlo.

Desde la posición del sol

El tiempo de sobrevuelo del satélite puede ser estimado de manera precisa por la posición del sol reportada por los metadatos, usando el módulo r.sunmask. El archivo de metadatos para el siguiente ejemplo contiene: SUN_AZIMUTH = 120.8810347, SUN_ELEVATION = 64.7730999.

1. de manera iterativa cambiar la hora y minuto para obtener un buen resultado, es necesario corregir la zona horaria:

r.sunmask  -s elev=elevation out=dummy year=2002 month=5 day=24 hour=10 min=42 sec=7 timezone=-5

El comando reporta:

Solar position: sun azimuth: 121.342461, sun angle above horz.(refraction corrected): 65.396652

El tiempo local de sobrevuelo resultante es 10:42:07 que corresponde a 15:42 en GMT, que corresponde a 15.70 en horas GMT decimales (minutos decimales: $42*100/60$).

A partir de la herramienta web de la NASA

El NASA LaRC Satellite Overpass Predictor puede calcular el tiempo de sobrevuelo del satélite a partir de la fecha de adquisición y las coordenadas centrales de la escena. El archivo de metadatos para el siguiente ejemplo contiene: ACQUISITION_DATE = 2002-05-24, SCENE_CENTER_LAT = +36.0512847, SCENE_CENTER_LON = -79.3280820.

1. Seleccione la zona del mundo en y de clic en Go >>.
2. Ingrese las coordenadas Lat/Long (grados decimales, sin signo más pero sí con el signo menos)
3. Seleccione el satélite adecuado
4. Ingrese la fecha de adquisición
5. Opcional: seleccione "Day" o "Night" para restingir el cálculo a sobrevuelos de día o de noche (locales) del satélite

Ejemplo: calculo para la escena Landsata de NC: vea aquí here.

--- Solamente para LANDSAT ---

El USGS provee la página web Landsat Bulk Metadata Service de donde es posible extraer los metadatos (incluyendo el tiempo de sobrevuelo) para todas las misiones Landsat.

Corrección atmosférica

Los valores de radiancia pueden ser convertidos a reflectancias al realizar la corrección atmosférica aplicando el algoritmo $6S$. El algoritmo transformará los valores de Radiancia por encima de la atmósfera a valores de Reflectancia por encima de la copa (Top-of-canopy), usando la información del contenido de aerosoles y la composición atmosférica del momento en que la imagen fue tomada. Lo que sigue describe el método para usar este algoritmo en GRASS GIS. De nuevo, solamente se hacen los cálculos para la banda 1, los número(s) a cambiar para otras bandas espectrales, se muestran en red. El archivo de control icnd_lsat1.txt consiste en los siguientes parámetros, y es escrito en un editor de texto:

8                          # indica que es una imagen ETM+
05 24 15.67 -78.691 35.749 # imagen tomada el 24 de mayo, a las 15:42 GMT en decimales; el centro de la imagen está en  78.691°O y 35.749°N 
2                          # the midlatitude summer atmospheric model
1                          # el modelo de aerosoles continentales
50                         # la visibilidad para el modelo de aerosoles [km]
-0.110                     # el terreno está 110 metros por encima del nivel del mar [km] * -1
-1000                      # imagen tomada de un sensor de satélite (1000)
61                         # banda espectral, aquí la 1: azul

De manera alternativa, en caso de que no se tenga un mapa de visibilidad, se puede usar una estimación de Profundidad óptica de aerosoles (Aerosol Optical Depth estimation) a 550nm. El archivo de control debe ser modificado de acuerdo a esto, ingresando 0 para la visibilidad, y en la siguiente línea, ingresar la profundidad óptica de aerosoles. El archivo de control debe ser estructurado de la siguiente manera:

8                          # indica que es una imagen ETM+
05 24 15.67 -78.691 35.749 # imagen tomada el 24 de mayo, a las 15:42 GMT en decimales; el centro de la imagen está en  78.691°O y 35.749°N 
2                          # the midlatitude summer atmospheric model
1                          # el modelo de aerosoles continentales
0                          # establece visibilidad = 0
0.112                      # ej. profundidad óptica de aerosoles
-0.110                     # el terreno está 110 metros por encima del nivel del mar [km] * -1
-1000                      # imagen tomada de un sensor de satélite (1000)
61                         # banda espectral, aquí la 1: azul

El archivo se usa después en el módulo i.atcorr:

i.atcorr -a -o iimg=lsat7_2002_10_rad ialt=elev_int icnd=icdn_lsat1.txt oimg=lsat7_2002_10_atcorr

o, usando GRASS 7.x:

i.atcorr -a input=lsat7_2002_10_rad elevation=elev_int parameters=icdn_lsat1.txt output=lsat7_2002_10_atcorr

Donde:

• -a se refiere a una imagen Landsat tomada luego de julio 2000
• -o activa la aceleración caché (solamente GRASS 6.x)
• iimg/input es la imagen a corregir
• ialt/elevation es el mapa de altitud que sobreescribe el valor inicial de 110 metros
• icnd/parameters es la ruta al archivo icnd.txt
• oimg/output es el nombre de la imagen de salida

Se puede encontrar información detallada en el manual de i.atcorr.

Nota, i.atcorr no se quejará si el modelo de elevación definido (parámetro elevation=) no cubre toda la extensión de la escena Landsat a que se sujeta el proceso -- dará valores NaN (Not a Number). [Maybe this can be requested as an enhancement?]

Fuentes para estimaciones de profundidad óptica de aerosoles

• AERONET - provee observaciones distribuidas globalmente de profundidad óptica de aerosoles (AOD), productos de inversión (inversion products), y agua precipitable en diversos regímenes de aerosoles.
• Producto MODIS MOD08 MOD08 - Gridded Atmospheric Product.
  • Los productos MOD08 diarios pueden ser descargados directamente a partir del servidor FTP de la NASA: ftp://ladsweb.nascom.nasa.gov/allData/51/MOD08_D3/.
  • Están disponibles las especificaciones del Archivo y capa, por ejemplo, en MODIS HDF File Specification MOD08_D3: MODIS Level 3 Daily Atmosphere Gridded Product
  • La capa Optical_Depth_Land_And_Ocean_Mean, descrita como Aerosol Optical Thickness at 0.55 microns for both Ocean (best) and Land (corrected): Mean, probablemente pueda ser usada como entrada para i.atcorr. Nota, los valores de la capa deben ser escalados usando el factor de escalamiento Optical_Depth_Land_And_Ocean_Mean:scale_factor = 0.001d, como se reporta en las especificaciones.
• How to estimate Aerosol Optical Thickness

Referencias

Algoritmo 6S

• 6S Web site
• Land Surface Reflectance Science Computing Facility website - 6S

Radiancia/reflectacia espectral

• How is radiance calculated?, question on USGS' Landsat Missions - Frequently Asked Questions
• (1) http://landsathandbook.gsfc.nasa.gov/handbook/handbook_htmls/chapter11/chapter11.html
• (2) http://landsat.usgs.gov/documents/L5TM_postcal_v11.pdf

Whitepapers

• Various Whitepapers about conversion and atmospheric correction of High Resolution Satellite Imagery
• Song, C., Woodcock, C.E., Seto, K.C., Lenney, M.P., Macomber, S.A., 2001. Classification and Change Detection Using Landsat TM Data: When and How to Correct Atmospheric Effects? Remote Sensing of Environment 75, 230-244. PDF

Herramientas on-line relacionadas

• NASA LaRC Satellite Overpass Predictor
• Landsat Bulk Metadata Service
• NASA's Atmospheric Correction Parameter Calculator

Cómo añadir sensores a i.atcorr

Cómo añadir sensores a i.atcorr:

• ver http://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.atcorr/README