View source for 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 {{cmd|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:

# 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 {{cmd|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 {{cmd|i.landsat.toar}} soporta todos los sensores Landsat incluyendo MSS, TM y ETM+.
# Para convertir '''Radiancia por encima de la atmósfera''' a '''Reflectancia por encima de la copa''' (ToCR) → use el módulo {{cmd|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 <code>>0</code>, 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:


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

== 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 {{cmd|i.landsat.toar}}).

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

<math>L\lambda = \frac{L_{MAX}λ - L_{MIN}λ}{QCAL_{MAX} - QCAL_{MIN}} * (QCAL - QCAL_{MIN}) + L_{MIN}\lambda</math>


<ul>donde:

* <math>Lλ</math> - Radiancia espectral en la apertura del sensor en watts/metro cuadrado * ster * μm), esta es la radiancia aparente registrada por el sensor;
* <math>QCAL</math> - el valor del pixel calibrado en  <math>DN</math>;
* <math>LMINλ</math> - la radiancia espectral que es escalada a <math>QCALMIN</math> en watts/(metro cuadrado * ster * μm);
* <math>LMAXλ</math> - la radiancia espectral que es escalada a <math>QCALMAX</math> en watts/(metro cuadrado * ster * μm);
* <math>QCALMINM</math> - el valor mínimo cuantificado de pixel (correspondiente a <math>LMINλ</math>) en <math>DN</math> ;
* <math>QCALMAX</math> - el valor mínimo cuantificado de pixel (correspondiente a <math>LMAXλ</math>) en <math>DN=255</math>.
</ul>


<math>LMINλ</math> y <math>LMAXλ</math> son las radiancias relacionadas a los valores máximos y mínimos de los números digitales <math>DN</math>, 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 <math>DN</math> (<math>QCALMIN</math>) es 1 para imágenes Landsat ETM+ (1), y el máximo (<math>QCALMAX</math>) es de 255. <math>QCAL</math> es el <math>DN</math> para cada pixel independiente en la imagen Landsat.

Accesando los metadatos:

 {{cmd|r.info}} <source lang="bash" enclose=none>lsat7_2002_xx</source>

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):

 <source lang="bash" enclose=none>CHAN=1</source>
 {{cmd|r.info}} <source lang="bash" enclose=none>lsat7_2002_${CHAN}0 -h | tr '\n' ' '  | sed 's+ ++g' | tr ':' '\n' | grep "LMIN_BAND${CHAN}\|LMAX_BAND${CHAN}"
</source>

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):

 {{cmd|g.region}} <source lang="bash" enclose=none>rast=lsat7_2002_10 -p</source>
 {{cmd|r.mapcalc}} <source lang="bash" enclose=none>"lsat7_2002_10_rad=((''191.6'' - (''-6.2'')) / (255.0 - 1.0)) * (lsat7_2002_10 - 1.0) + (''-6.2'')"</source>

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

  {{cmd|r.mapcalc}} <source lang="bash" enclose=none>"elev_int = round(elevation)"</source>

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

 {{cmd|r.univar}} <source lang="bash" enclose=none>elev_int</source>

== Estimando el tiempo de sobrevuelo ==
Este valor es necesario para el archivo de control. Para más detalles, vea la wikipedia en [http://en.wikipedia.org/wiki/Decimal_time#Scientific_decimal_time 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 {{cmd|r.sunmask}}.
El  [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/2002/p016r035_7x20020524.met.gz archivo de metadatos] para el siguiente ejemplo contiene: <source lang="bash" enclose="none">SUN_AZIMUTH = 120.8810347, SUN_ELEVATION = 64.7730999</source>.

# de manera iterativa cambiar la hora y minuto para obtener un buen resultado, es necesario corregir la zona horaria:
 {{cmd|r.sunmask}}  <source lang="bash" enclose=none>-s elev=elevation out=dummy year=2002 month=5 day=24 hour=10 min=42 sec=7 timezone=-5</source>

El comando reporta:

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

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

=== A partir de la herramienta web de la NASA ===

El [http://cloudsgate2.larc.nasa.gov/cgi-bin/predict/predict.cgi 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 [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/2002/p016r035_7x20020524.met.gz archivo de metadatos] para el siguiente ejemplo contiene: <source lang="bash" enclose="none">ACQUISITION_DATE = 2002-05-24, SCENE_CENTER_LAT = +36.0512847, SCENE_CENTER_LON = -79.3280820</source>.


# Seleccione la zona del mundo en y de clic en <code>Go >></code>.
# Ingrese las coordenadas Lat/Long (grados decimales, sin signo ''más'' pero sí con el signo ''menos'')
# Seleccione el satélite adecuado
# Ingrese la fecha de adquisición
# ''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í [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/landsat_overpass_time_list_NC.txt here].

--- Solamente para LANDSAT ---

El USGS provee la página web [http://landsat.usgs.gov/consumer.php 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 <math>6S</math>. 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 <span style="color:#FF0000">red</span>. El archivo de control <code>icnd_lsat1.txt</code> 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)
 <span style="color:#FF0000">61</span>                         # 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)
 <span style="color:#FF0000">61</span>                         # banda espectral, aquí la 1: azul


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

 {{cmd|i.atcorr}} <source lang="bash" enclose=none>-a -o iimg=lsat7_2002_10_rad ialt=elev_int icnd=icdn_lsat1.txt oimg=lsat7_2002_10_atcorr</source>

o, usando GRASS 7.x:

 {{cmd|i.atcorr}} <source lang="bash" enclose=none>-a input=lsat7_2002_10_rad elevation=elev_int parameters=icdn_lsat1.txt output=lsat7_2002_10_atcorr</source>

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

Se puede encontrar información detallada en el manual de {{cmd|i.atcorr}}.

'''Nota,''' {{cmd|i.atcorr}} no se quejará si el modelo de elevación definido (parámetro <code>elevation=</code>) 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  ==

* [http://aeronet.gsfc.nasa.gov 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 [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=08 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/ '''ftp://ladsweb.nascom.nasa.gov/allData/51/MOD08_D3/'''].
** Están disponibles las especificaciones del Archivo y capa, por ejemplo, en [http://www.atmos.washington.edu/~robwood/modis/filespecs/MOD08_D3.CDL.fs 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 {{cmd|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.
* [http://www.cesbio.ups-tlse.fr/multitemp/?p=1710 How to estimate Aerosol Optical Thickness]

<!--
=== TοDο ===
* entender si el  [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=04 MODIS Aerosol Product (MOD04)] could be useful to estimate the visibility for the aerosol model.
* provide an example on using [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=08 MODIS Gridded Atmospheric Product (MOD08)].
* explain if it is possible and how to use ''Rayleigh Optical Depths estimated at 0km altitude for six different atmosphere models'' from the paper [http://www.opticsinfobase.org/ao/abstract.cfm?uri=ao-34-15-2765  Rayleigh-scattering calculations for the terrestrial atmosphere] by Anthony Buchholz, table 4, page 2270 -->

== Referencias ==

=== Algoritmo 6S ===

* [http://6s.ltdri.org 6S Web site]
* [http://modis-sr.ltdri.org/ Land Surface Reflectance Science Computing Facility website - 6S]

=== Radiancia/reflectacia espectral ===

* [http://landsat.usgs.gov/how_is_radiance_calculated.php How is radiance calculated?], question on USGS' [http://landsat.usgs.gov/tools_faq.php 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 [http://apollomapping.com/about-us/whitepapers 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. [http://www.unc.edu/courses/2008spring/geog/577/001/www/Song01_RSE.pdf PDF]

=== Herramientas on-line relacionadas ===

* [http://cloudsgate2.larc.nasa.gov/cgi-bin/predict/predict.cgi NASA LaRC Satellite Overpass Predictor]
* [http://landsat.usgs.gov/consumer.php Landsat Bulk Metadata Service]
* NASA's [http://atmcorr.gsfc.nasa.gov/ Atmospheric Correction Parameter Calculator]

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

Cómo añadir sensores a {{cmd|i.atcorr}}:
<ul>
* ver http://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.atcorr/README
</ul>

[[Category: Documentation]]
[[Category: Tutorial]]
[[Category: Image processing]]
[[Category: Languages/es]]