Difference between revisions of "Atmospheric correction"

From GRASS-Wiki
Jump to: navigation, search
(+ See also: i.topo.corr)
 
(42 intermediate revisions by 4 users not shown)
Line 4: Line 4:
  
 
* Transforming a Landsat ETM+ image from DN values to radiance values
 
* Transforming a Landsat ETM+ image from DN values to radiance values
* Conducting an atmospheric correction of a Landsat image using the {{cmd|i.atcorr}} module in GRASS GIS
+
* Conducting an atmospheric correction of a Landsat image using the {{cmd|i.atcorr}} module in GRASS GIS and convert radiance to reflectance values
  
 
== Introduction ==
 
== Introduction ==
  
The standard steps for processing satellite images are:
+
Standard steps for pre-processing satellite images are:
  
# Convert DN (digital number = pixel values) to ''Radiance at TOA'' (see formula below or use {{cmd|i.landsat.toar|version=65}}). {{cmd|i.landsat.toar|version=65}} supports all Landsat versions from MSS to TM to ETM+.
+
# Converting '''Digital Numbers''' (DN) to '''Top-of-Atmosphere Radiances''' (ToAR)
# Perform atmospheric correction (which converts Radiance at TOA -> ''Reflectance at surface''): done by {{cmd|i.atcorr}}.
+
#* for Landsat imagery, apply the equation below or use the {{cmd|i.landsat.toar}} module
 +
#* DNs -- which are actually the pixel values -- are the result of the quantified amount of energy observed and measured at the sensor
 +
#* the {{cmd|i.landsat.toar}} module supports all Landsat sensors icluding MSS, TM and ETM+
 +
# Converting '''ToA Radiances''' to '''Top-of-Canopy Reflectances''' (ToCR) → use the {{cmd|i.atcorr}} module
 +
#* the conversion to ''Reflectances'' attempts to assess for and remove atmospheric effects, a process known as ''atmospheric correction''
 +
#* the ''atmospheric correction'' step is quite a common topic of discussion in Remote Sensing. However, depending on the quality of the acquisition(s), it can be a complex process that does not always offer a great benefit, in terms of accuracy for subsequent processing steps (e.g. image segmentation and classification).
  
''Note:'' Remember, to check if water areas are > 0, since reflectance is always positive. If negative, you have a "cornichon" (pickle in English ;-0). Means the transform equation inputs are not belong to image processed.
+
'''Remember''' to check if water areas are represented by values <code>>0</code>, since reflectance is always positive. If you encounter negative reflectance values, you have a "cornichon" (pickle in English ;-0). This means that the transformation equations used, do not correspond to the image processed.
 +
 
 +
 
 +
An overview of options to get "corrected" Spectral Reflectances:
 +
 
 +
  <nowiki>+------------------------------------------------------------------------------+
 +
|                                                                              |
 +
| Digital Numbers                                                              |
 +
|        |                                                                    |
 +
|  +-----v-----+                                                              |
 +
|  |  i.*.toar | ---> Reflectance                                              |
 +
|  +-----+-----+    (uncorrected)~~~(DOS methods)--+                          |
 +
|        |                +                        |                          |
 +
|    (-r flag)        (-r flag)                    +--> "Corrected"          |
 +
|        |                |                        |    Reflectance          |
 +
|        v          +-----v----+                  |                          |
 +
|    Radiance ------> i.atcorr +-------------------+                          |
 +
|                    +----------+                                              |
 +
|                                                                              |
 +
+------------------------------------------------------------------------------+</nowiki>
  
 
== Requirements ==
 
== Requirements ==
  
 
* GRASS 6.4.0 or higher
 
* GRASS 6.4.0 or higher
* North Carolina sample dataset (location): http://grass.osgeo.org/download/data.php
+
* North Carolina sample dataset (location): https://grass.osgeo.org/download/sample-data/
  
 
After downloading of the North Carolina dataset, it must be unzipped into your GIS database. When starting GRASS GIS, select this GIS database, open the ‘nc_spm_08’ dataset as LOCATION, and ‘PERMANENT’ as MAPSET.  
 
After downloading of the North Carolina dataset, it must be unzipped into your GIS database. When starting GRASS GIS, select this GIS database, open the ‘nc_spm_08’ dataset as LOCATION, and ‘PERMANENT’ as MAPSET.  
Line 24: Line 48:
 
== Calculating radiance values ==
 
== Calculating radiance values ==
  
The NC data set mapset contains, amongst others, a Landsat ETM+ image of 24th of May 2002. Every pixel of this image contains a DN or grey value. In order to be able to make calculations with satellite imagery, or compare values amongst different sensors, these values have to be converted to radiances or reflectances. The formulas to do this conversion are described here for Landsat images (or use [[GRASS_AddOns#i.landsat.toar|i.landsat.toar]]).
+
The NC data set mapset contains, amongst others, a Landsat ETM+ image of 24th of May 2002. Every pixel of this image contains a DN or grey value. In order to be able to make calculations with satellite imagery, or compare values amongst different sensors, these values have to be converted to radiances or reflectances. The formulas to do this conversion are described here for Landsat images (or use {{cmd|i.landsat.toar}}).
 +
 
 +
 
 
The conversion from DN values to top-of-atmosphere radiances is done as follows:
 
The conversion from DN values to top-of-atmosphere radiances is done as follows:
  
<!-- why doesn't TeX markup work??
+
<math>L\lambda = \frac{L_{MAX}λ - L_{MIN}λ}{QCAL_{MAX} - QCAL_{MIN}} * (QCAL - QCAL_{MIN}) + L_{MIN}\lambda</math>
<math>L\lambda = \frac{L_{MAX}λ - L_{MIN}λ}{QCAL_{MAX}-QCAL_{MIN}} * (QCAL-QCAL_{MIN}) + L_{MIN}\lambda</math>
+
 
  
<math>\sqrt{2}</math>
+
<ul>where:
-->
+
 
Lλ = ((LMAXλ - LMINλ)/(QCALMAX-QCALMIN)) * (QCAL-QCALMIN) + LMINλ
+
* <math></math> - spectral Radiance at the sensor's aperture in watts/(meter squared * ster * μm), the apparent radiance as seen by the satellite sensor ;
 +
* <math>QCAL</math> - the quantized calibrated pixel value in <math>DN</math> ;
 +
* <math>LMINλ</math> - the spectral radiance that is scaled to <math>QCALMIN</math> in watts/(meter squared * ster * μm) ;
 +
* <math>LMAXλ</math> - the spectral radiance that is scaled to <math>QCALMAX</math> in watts/(meter squared * ster * μm) ;
 +
* <math>QCALMINM</math> - the minimum quantized calibrated pixel value (corresponding to <math>LMINλ</math>) in <math>DN</math> ;
 +
* <math>QCALMAX</math> - the maximum quantized calibrated pixel value (corresponding to <math>LMAXλ</math>) in <math>DN=255</math>.
 +
</ul>
  
Where:
 
  
* Lλ - spectral Radiance at the sensor's aperture in watts/(meter squared * ster * μm), the apparent radiance as seen by the satellite sensor;
+
<math>LMINλ</math> and <math>LMAXλ</math> are the radiances related to the minimal and maximal <math>DN</math> value, and are reported in the metadata file for each image, or in the table 1. High gain or low gain is also reported in the metadata file of each Landsat image. The minimal <math>DN</math> value (<math>QCALMIN</math>) is 1 for Landsat ETM+ images (1), and the maximal <math>DN</math> value (<math>QCALMAX</math>) is 255. <math>QCAL</math> is the <math>DN</math> value for every separate pixel in the Landsat image.
* QCAL - the quantized calibrated pixel value in DN;
 
* LMINλ - the spectral radiance that is scaled to QCALMIN in watts/(meter squared * ster * μm);
 
* LMAXλ - the spectral radiance that is scaled to QCALMAX in watts/(meter squared * ster * μm);
 
* QCALMIN - the minimum quantized calibrated pixel value (corresponding to LMINλ) in DN;
 
* QCALMAX - the maximum quantized calibrated pixel value (corresponding to LMAXλ) in DN=255.
 
  
LMINλ and LMAXλ are the radiances related to the minimal and maximal DN value, and are reported in the metadata file for each image, or in the table 1. High gain or low gain is also reported in the metadata file of each Landsat image. The minimal DN value (QCALMIN) is 1 for Landsat ETM+ images (1), and the maximal DN value (QCALMAX) is 255. QCAL is the DN value for every separate pixel in the Landsat image.
 
  
 
Accessing the metadata:
 
Accessing the metadata:
  
  {{cmd|r.info}} lsat7_2002_xx
+
  {{cmd|r.info}} <source lang="bash" enclose=none>lsat7_2002_xx</source>
  
 
Under ‘comments’, the maximal and minimal radiance (LMAX and LMIN) for each band are given.  
 
Under ‘comments’, the maximal and minimal radiance (LMAX and LMIN) for each band are given.  
Line 53: Line 78:
 
Power users may use this for nicely readable output (example for channel 1):
 
Power users may use this for nicely readable output (example for channel 1):
  
  CHAN=1
+
  <source lang="bash" enclose=none>CHAN=1</source>
  {{cmd|r.info}} lsat7_2002_${CHAN}0 -h | tr '\n' ' '  | sed 's+ ++g' | tr ':' '\n' | grep "LMIN_BAND${CHAN}\|LMAX_BAND${CHAN}"
+
  {{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>
  
 
Conversion to radiance (this calculation is done for band 1, for the other bands, the numbers in italics need to be replaced with their correct values):
 
Conversion to radiance (this calculation is done for band 1, for the other bands, the numbers in italics need to be replaced with their correct values):
  
  {{cmd|g.region}} rast=lsat7_2002_10 -p
+
  {{cmd|g.region}} <source lang="bash" enclose=none>rast=lsat7_2002_10 -p</source>
  {{cmd|r.mapcalc}} "lsat7_2002_10_rad=((''191.6'' - (''-6.2'')) / (255.0 - 1.0)) * (lsat7_2002_10 - 1.0) + (''-6.2'')"
+
  {{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>
  
 
For faster computations, use an integer DEM:
 
For faster computations, use an integer DEM:
  
   {{cmd|r.mapcalc}} "elev_int = round(elevation)"
+
   {{cmd|r.mapcalc}} <source lang="bash" enclose=none>"elev_int = round(elevation)"</source>
  
 
Find mean elevation to initialize computations (needed in 'icnd' file below)
 
Find mean elevation to initialize computations (needed in 'icnd' file below)
  
  {{cmd|r.univar}} elev_int
+
  {{cmd|r.univar}} <source lang="bash" enclose=none>elev_int</source>
  
== Estimating the overpass time from the sun position ==
+
== Estimating the overpass time ==
 +
This value is needed for the control file. For details, refer to Wikipedia's entry on [http://en.wikipedia.org/wiki/Decimal_time#Scientific_decimal_time Scientific decimal time].
 +
In case this parameter is not reported in the metadata file, we have two ways to obtain it.
  
The satellite overpass time can be estimated rather precisely from the sun position reported in metadata using {{cmd|r.sunmask}}:
+
=== From sun position ===
The [ftp://ftp.glcf.umiacs.umd.edu/glcf/Landsat/WRS2/p016/r035/p016r035_7x20020524.ETM-EarthSat-Orthorectified/p016r035_7x20020524.met metadata] of this example report: SUN_AZIMUTH = 120.8810347, SUN_ELEVATION = 64.7730999
 
  
# iteratively change hour and minutes to obtain a good result, timezone needs to be correct of course:
+
The satellite overpass time can be estimated rather precisely from the sun position reported in metadata using {{cmd|r.sunmask}}.
{{cmd|r.sunmask}} -s elev=elevation out=dummy year=2002 month=5 day=24 hour=10 min=42 sec=7 timezone=-5
+
The [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/2002/p016r035_7x20020524.met.gz metadata file] for the following example contains: <source lang="bash" enclose="none">SUN_AZIMUTH = 120.8810347, SUN_ELEVATION = 64.7730999</source>.
  
Reports: sun azimuth: 121.342461, sun angle above horz.(refraction corrected): 65.396652
+
# iteratively change hour and minutes to obtain a good result, timezone needs to be correct of course
 +
{{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>
  
The resulting overpass local time 10:42:07 corresponds to 15:42 in GMT which corresponds to 15.67 in decimal GMT hours (decimal minutes: 42 * 100 / 60). This value is needed for the control file.
+
The above command reports:
 +
 
 +
Solar position: sun azimuth: 121.342461, sun angle above horz.(refraction corrected): 65.396652
 +
 
 +
The resulting overpass local time <code>10:42:07</code> corresponds to <code>15:42</code> in GMT which corresponds to <code>15.70</code> in decimal GMT hours (decimal minutes: <math>42 * 100 / 60</math>).
 +
 
 +
=== From NASA web tool ===
 +
 
 +
The [http://cloudsgate2.larc.nasa.gov/cgi-bin/predict/predict.cgi NASA LaRC Satellite Overpass Predictor] can compute satellite overpass time from d''ate of acquisition'' and ''scene center coordinates''. The [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/2002/p016r035_7x20020524.met.gz metadata file] for the following example contains: <source lang="bash" enclose="none">ACQUISITION_DATE = 2002-05-24, SCENE_CENTER_LAT = +36.0512847, SCENE_CENTER_LON = -79.3280820</source>.
 +
 
 +
# Select the World zone in the bottom-right side of the page and click <code>Go >></code>
 +
# Insert Lat/Long coordinates (decimal degrees, without ''plus'' sign but with ''minus'' sign)
 +
# Select the proper satellite
 +
# Insert the date of acquisition
 +
# ''Optional:'' select "Day" or "Night" to restrict computation to (local) day/night overpass of the satellite
 +
 
 +
''Example:'' Calculation for NC Landsat scenes: see [http://www.grassbook.org/wp-content/uploads/ncexternal/landsat/landsat_overpass_time_list_NC.txt here]
 +
 
 +
--- For LANDSAT only ---
 +
 
 +
USGS provide the [http://landsat.usgs.gov/consumer.php Landsat Bulk Metadata Service] web page where is possible to extract metadata (included overpass time) for all the Landsat missions.
  
 
== Atmospheric correction ==
 
== Atmospheric correction ==
  
This radiance image can be used for the atmospheric correction with the '''6S algorithm'''. The algorithm will transform the top-of-the-atmosphere radiance values to top-of-canopy reflectance values using predefined information on the aerosol content and atmospheric composition of the time the image was taken. What follows describes the method to use this algorithm in GRASS GIS. Again, only the calculations for band 1 are shown, and the numbers to be changed for the other bands are indicated in red. The 'icnd_lsat1.txt' control file consists of the following parameters, and is written with a text editor:
+
Radiance values can be converted to reflectance values in undergoing an atmospheric correction by applying the <math>6S</math> algorithm available in {{cmd|i.atcorr}}. The algorithm will transform the Top-of-Atmosphere radiance values to Top-of-Canopy reflectance values using predefined information on the aerosol content and atmospheric composition of the time the image was taken. What follows describes the method to use this algorithm in GRASS GIS. Again, only calculations for band 1 are demonstrated here, and the number(s) to be changed for indicating other spectral bands, are <span style="color:#FF0000">red</span>. The <code>icnd_lsat1.txt</code> control file consists of the following parameters, and is written with a text editor:
  
 
  8                          # indicates that it is an ETM+ image
 
  8                          # indicates that it is an ETM+ image
Line 92: Line 140:
 
  -0.110                    # the terrain lies 110 meters above sea level [km] * -1
 
  -0.110                    # the terrain lies 110 meters above sea level [km] * -1
 
  -1000                      # image taken of a satellite sensor (1000)
 
  -1000                      # image taken of a satellite sensor (1000)
  61                        # spectral band, here 1: blue
+
  <span style="color:#FF0000">61</span>                         # spectral band, here 1: blue
  
''TODO: understand if the [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.''
+
Alternatively, in case visibility maps are not available, an Aerosol Optical Depth estimation at 550nm can be used. The control file should be altered accordingly by entering 0 for the visibility and, in a following line, entering the aerosol optical depth. The control file should be structured as follows:
 +
 
 +
8                          # indicates that it is an ETM+ image
 +
05 24 15.67 -78.691 35.749 # image taken on the 24th of May, at 15:42 GMT in decimals; the center of the image lies at 78.691°W and 35.749°N
 +
2                          # the midlatitude summer atmospheric model
 +
1                          # the continental aerosol model
 +
0                          # set visibility = 0
 +
0.112                      # e.g. aerosol optical depth
 +
-0.110                    # the terrain lies 110 meters above sea level [km] * -1
 +
-1000                      # image taken of a satellite sensor (1000)
 +
<span style="color:#FF0000">61</span>                        # spectral band, here 1: blue
  
  
 
This file is then used in the {{cmd|i.atcorr}} module:
 
This file is then used in the {{cmd|i.atcorr}} module:
  
  {{cmd|i.atcorr}} -a -o iimg=lsat7_2002_10_rad ialt=elev_int icnd=icdn_lsat1.txt oimg=lsat7_2002_10_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>
 +
 
 +
or, using 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>
  
 
Where:
 
Where:
* <tt>-a</tt> refers to a Landsat image taken after July 2000
+
* <tt>'''-a'''</tt>   refers to a Landsat image taken after July 2000
* <tt>-o</tt> activates the cache acceleration
+
* <tt>'''-o'''</tt>   activates the cache acceleration (GRASS 6.x only)
* <tt>iimg</tt> is the image to be corrected
+
* <tt>'''iimg/input'''</tt> is the image to be corrected
* <tt>ialt</tt> is the altitude map which overrides the initialization value of 110 meters
+
* <tt>'''ialt/elevation'''</tt> is the altitude map which overrides the initialization value of 110 meters
* <tt>icnd</tt> is the path to the icnd.txt file
+
* <tt>'''icnd/parameters'''</tt> is the path to the icnd.txt file
* <tt>oimg</tt> is the name of the output image
+
* <tt>'''oimg/output'''</tt> is the name of the output image
 +
 
 +
Detailed information can be found in {{cmd|i.atcorr}}'s manual.
 +
 
 +
'''Note,''' {{cmd|i.atcorr}} wont complain if the defined elevation model (parameter <code>elevation=</code>) does not cover the extent of the Landsat scene that is the subject of the process -- it will result in NaN's. ''[Maybe this can be requested as an enhancement?]''
 +
 
 +
== Sources for aerosol optical depth (AOD) estimations ==
 +
 
 +
* [http://aeronet.gsfc.nasa.gov AERONET] - provides globally distributed observations of spectral aerosol optical depth (AOD), inversion products, and precipitable water in diverse aerosol regimes.
 +
* MODIS MOD08 product [http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php?MOD_NUMBER=08 MOD08 - Gridded Atmospheric Product].
 +
** Daily MOD08 products can be downloaded directly from NASA's FTP server: [ftp://ladsweb.nascom.nasa.gov/allData/51/MOD08_D3/ '''ftp://ladsweb.nascom.nasa.gov/allData/51/MOD08_D3/'''].
 +
** File and layer specifications are available, for example, at [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]
 +
** The layer '''Optical_Depth_Land_And_Ocean_Mean''', described as ''Aerosol Optical Thickness at 0.55 microns for both Ocean (best) and Land (corrected): Mean'', can be probably used as an input for {{cmd|i.atcorr}}. '''Note,''' the layer's values have to be scaled by using the scaling factor ''Optical_Depth_Land_And_Ocean_Mean:scale_factor ='' '''0.001d''' as reported in the specifications.
 +
* [http://www.cesbio.ups-tlse.fr/multitemp/?p=1710 How to estimate Aerosol Optical Thickness]
  
More information on the use of {{cmd|i.atcorr}} for other images can be found in the {{cmd|i.atcorr}} manual.
+
=== TοDο ===
 +
* Cross-check with [http://atmcorr.gsfc.nasa.gov  Landsat Atmospheric Correction Parameter Calculator]
 +
* understand if the [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
  
 
== References ==
 
== References ==
 +
 +
=== 6S algorithm ===
  
 
* [http://6s.ltdri.org 6S Web site]
 
* [http://6s.ltdri.org 6S Web site]
 
* [http://modis-sr.ltdri.org/ Land Surface Reflectance Science Computing Facility website - 6S]
 
* [http://modis-sr.ltdri.org/ Land Surface Reflectance Science Computing Facility website - 6S]
 +
 +
=== Spectral radiance/reflectance ===
 +
 +
* [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
 
* (1)  http://landsathandbook.gsfc.nasa.gov/handbook/handbook_htmls/chapter11/chapter11.html
 
* (2)  http://landsat.usgs.gov/documents/L5TM_postcal_v11.pdf
 
* (2)  http://landsat.usgs.gov/documents/L5TM_postcal_v11.pdf
* [http://www-air.larc.nasa.gov/tools/predict.htm NASA LaRC Satellite Overpass Predictor]
 
  
Adding new sensors:
+
=== 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]
 +
 
 +
=== Related on-line tools ===
 +
 
 +
* [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]
 +
 
 +
=== How to add new sensors to i.atcorr ===
 +
 
 +
How to add new sensors to {{cmd|i.atcorr}}:
 +
<ul>
 
* see http://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.atcorr/README
 
* see http://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.atcorr/README
 +
</ul>
  
 
[[Category: Documentation]]
 
[[Category: Documentation]]
Line 126: Line 225:
 
[[Category: Image processing]]
 
[[Category: Image processing]]
  
 +
== See also ==
  
==Other==
+
* Terrain correction: {{cmd|i.topo.corr}}
[http://thetvtopc.com/Reverse_Cell_Phone_Lookup_Number reverse phone lookup]
 

Latest revision as of 05:06, 16 February 2018

See also the LANDSAT wiki page.

Aim of this tutorial:

  • Transforming a Landsat ETM+ image from DN values to radiance values
  • Conducting an atmospheric correction of a Landsat image using the i.atcorr module in GRASS GIS and convert radiance to reflectance values

Introduction

Standard steps for pre-processing satellite images are:

  1. Converting Digital Numbers (DN) to Top-of-Atmosphere Radiances (ToAR)
    • for Landsat imagery, apply the equation below or use the i.landsat.toar module
    • DNs -- which are actually the pixel values -- are the result of the quantified amount of energy observed and measured at the sensor
    • the i.landsat.toar module supports all Landsat sensors icluding MSS, TM and ETM+
  2. Converting ToA Radiances to Top-of-Canopy Reflectances (ToCR) → use the i.atcorr module
    • the conversion to Reflectances attempts to assess for and remove atmospheric effects, a process known as atmospheric correction
    • the atmospheric correction step is quite a common topic of discussion in Remote Sensing. However, depending on the quality of the acquisition(s), it can be a complex process that does not always offer a great benefit, in terms of accuracy for subsequent processing steps (e.g. image segmentation and classification).

Remember to check if water areas are represented by values >0, since reflectance is always positive. If you encounter negative reflectance values, you have a "cornichon" (pickle in English ;-0). This means that the transformation equations used, do not correspond to the image processed.


An overview of options to get "corrected" Spectral Reflectances:

 +------------------------------------------------------------------------------+
 |                                                                              |
 | Digital Numbers                                                              |
 |        |                                                                     |
 |  +-----v-----+                                                               |
 |  |  i.*.toar | ---> Reflectance                                              |
 |  +-----+-----+     (uncorrected)~~~(DOS methods)--+                          |
 |        |                 +                        |                          |
 |    (-r flag)         (-r flag)                    +--> "Corrected"           |
 |        |                 |                        |     Reflectance          |
 |        v           +-----v----+                   |                          |
 |     Radiance ------> i.atcorr +-------------------+                          |
 |                    +----------+                                              |
 |                                                                              |
 +------------------------------------------------------------------------------+

Requirements

After downloading of the North Carolina dataset, it must be unzipped into your GIS database. When starting GRASS GIS, select this GIS database, open the ‘nc_spm_08’ dataset as LOCATION, and ‘PERMANENT’ as MAPSET.

Calculating radiance values

The NC data set mapset contains, amongst others, a Landsat ETM+ image of 24th of May 2002. Every pixel of this image contains a DN or grey value. In order to be able to make calculations with satellite imagery, or compare values amongst different sensors, these values have to be converted to radiances or reflectances. The formulas to do this conversion are described here for Landsat images (or use i.landsat.toar).


The conversion from DN values to top-of-atmosphere radiances is done as follows:

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


    where:
    • Lλ - spectral Radiance at the sensor's aperture in watts/(meter squared * ster * μm), the apparent radiance as seen by the satellite sensor ;
    • QCAL - the quantized calibrated pixel value in DN ;
    • LMINλ - the spectral radiance that is scaled to QCALMIN in watts/(meter squared * ster * μm) ;
    • LMAXλ - the spectral radiance that is scaled to QCALMAX in watts/(meter squared * ster * μm) ;
    • QCALMINM - the minimum quantized calibrated pixel value (corresponding to LMINλ) in DN ;
    • QCALMAX - the maximum quantized calibrated pixel value (corresponding to LMAXλ) in DN=255.


LMINλ and LMAXλ are the radiances related to the minimal and maximal DN value, and are reported in the metadata file for each image, or in the table 1. High gain or low gain is also reported in the metadata file of each Landsat image. The minimal DN value (QCALMIN) is 1 for Landsat ETM+ images (1), and the maximal DN value (QCALMAX) is 255. QCAL is the DN value for every separate pixel in the Landsat image.


Accessing the metadata:

r.info lsat7_2002_xx

Under ‘comments’, the maximal and minimal radiance (LMAX and LMIN) for each band are given.

Power users may use this for nicely readable output (example for channel 1):

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

Conversion to radiance (this calculation is done for band 1, for the other bands, the numbers in italics need to be replaced with their correct values):

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'')"

For faster computations, use an integer DEM:

 r.mapcalc "elev_int = round(elevation)"

Find mean elevation to initialize computations (needed in 'icnd' file below)

r.univar elev_int

Estimating the overpass time

This value is needed for the control file. For details, refer to Wikipedia's entry on Scientific decimal time. In case this parameter is not reported in the metadata file, we have two ways to obtain it.

From sun position

The satellite overpass time can be estimated rather precisely from the sun position reported in metadata using r.sunmask. The metadata file for the following example contains: SUN_AZIMUTH = 120.8810347, SUN_ELEVATION = 64.7730999.

# iteratively change hour and minutes to obtain a good result, timezone needs to be correct of course
r.sunmask  -s elev=elevation out=dummy year=2002 month=5 day=24 hour=10 min=42 sec=7 timezone=-5

The above command reports:

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

The resulting overpass local time 10:42:07 corresponds to 15:42 in GMT which corresponds to 15.70 in decimal GMT hours (decimal minutes: 42 * 100 / 60).

From NASA web tool

The NASA LaRC Satellite Overpass Predictor can compute satellite overpass time from date of acquisition and scene center coordinates. The metadata file for the following example contains: ACQUISITION_DATE = 2002-05-24, SCENE_CENTER_LAT = +36.0512847, SCENE_CENTER_LON = -79.3280820.

  1. Select the World zone in the bottom-right side of the page and click Go >>
  2. Insert Lat/Long coordinates (decimal degrees, without plus sign but with minus sign)
  3. Select the proper satellite
  4. Insert the date of acquisition
  5. Optional: select "Day" or "Night" to restrict computation to (local) day/night overpass of the satellite

Example: Calculation for NC Landsat scenes: see here

--- For LANDSAT only ---

USGS provide the Landsat Bulk Metadata Service web page where is possible to extract metadata (included overpass time) for all the Landsat missions.

Atmospheric correction

Radiance values can be converted to reflectance values in undergoing an atmospheric correction by applying the 6S algorithm available in i.atcorr. The algorithm will transform the Top-of-Atmosphere radiance values to Top-of-Canopy reflectance values using predefined information on the aerosol content and atmospheric composition of the time the image was taken. What follows describes the method to use this algorithm in GRASS GIS. Again, only calculations for band 1 are demonstrated here, and the number(s) to be changed for indicating other spectral bands, are red. The icnd_lsat1.txt control file consists of the following parameters, and is written with a text editor:

8                          # indicates that it is an ETM+ image
05 24 15.67 -78.691 35.749 # image taken on the 24th of May, at 15:42 GMT in decimals; the center of the image lies at 78.691°W and 35.749°N 
2                          # the midlatitude summer atmospheric model
1                          # the continental aerosol model
50                         # the visibility for the aerosol model [km]
-0.110                     # the terrain lies 110 meters above sea level [km] * -1
-1000                      # image taken of a satellite sensor (1000)
61                         # spectral band, here 1: blue

Alternatively, in case visibility maps are not available, an Aerosol Optical Depth estimation at 550nm can be used. The control file should be altered accordingly by entering 0 for the visibility and, in a following line, entering the aerosol optical depth. The control file should be structured as follows:

8                          # indicates that it is an ETM+ image
05 24 15.67 -78.691 35.749 # image taken on the 24th of May, at 15:42 GMT in decimals; the center of the image lies at 78.691°W and 35.749°N 
2                          # the midlatitude summer atmospheric model
1                          # the continental aerosol model
0                          # set visibility = 0
0.112                      # e.g. aerosol optical depth
-0.110                     # the terrain lies 110 meters above sea level [km] * -1
-1000                      # image taken of a satellite sensor (1000)
61                         # spectral band, here 1: blue


This file is then used in the i.atcorr module:

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

or, using GRASS 7.x:

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

Where:

  • -a refers to a Landsat image taken after July 2000
  • -o activates the cache acceleration (GRASS 6.x only)
  • iimg/input is the image to be corrected
  • ialt/elevation is the altitude map which overrides the initialization value of 110 meters
  • icnd/parameters is the path to the icnd.txt file
  • oimg/output is the name of the output image

Detailed information can be found in i.atcorr's manual.

Note, i.atcorr wont complain if the defined elevation model (parameter elevation=) does not cover the extent of the Landsat scene that is the subject of the process -- it will result in NaN's. [Maybe this can be requested as an enhancement?]

Sources for aerosol optical depth (AOD) estimations

TοDο

References

6S algorithm

Spectral radiance/reflectance

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

Related on-line tools

How to add new sensors to i.atcorr

How to add new sensors to i.atcorr:

See also