Provided by: grass-doc_6.4.3-3_all bug

NAME

       i.atcorr  - Performs atmospheric correction using the 6S algorithm.
       6S - Second Simulation of Satellite Signal in the Solar Spectrum.

KEYWORDS

       imagery, atmospheric correction

SYNOPSIS

       i.atcorr
       i.atcorr help
       i.atcorr   [-frabo]   iimg=name   [iscl=min,max]    [ialt=name]    [ivis=name]   icnd=name
       oimg=name  [oscl=min,max]   [--overwrite]  [--verbose]  [--quiet]

   Flags:
       -f
           Output raster is floating point

       -r
           Input map converted to reflectance (default is radiance)

       -a
           Input from ETM+ image taken after July 1, 2000

       -b
           Input from ETM+ image taken before July 1, 2000

       -o
           Try to increase computation speed when altitude and/or visibility map is used

       --overwrite
           Allow output files to overwrite existing files

       --verbose
           Verbose module output

       --quiet
           Quiet module output

   Parameters:
       iimg=name
           Name of input raster map

       iscl=min,max
           Input imagery range [0,255]
           Default: 0,255

       ialt=name
           Input altitude raster map in m (optional)

       ivis=name
           Input visibility raster map in km (optional)

       icnd=name
           Name of input text file

       oimg=name
           Name for output raster map

       oscl=min,max
           Rescale output raster map [0,255]
           Default: 0,255

DESCRIPTION

       i.atcorr performs atmospheric correction on the input raster map using  the  6S  algorithm
       (Second  Simulation  of  Satellite  Signal  in  the  Solar Spectrum). A detailed algorithm
       description is available at  the  Land  Surface  Reflectance  Science  Computing  Facility
       website.

       Important  note:  Current region settings are ignored! The region is adjusted to cover the
       input raster map before the atmospheric correction is performed. The previous settings are
       restored afterwards.

       Because  using  a elevation and/or visibility raster map makes execution time much longer,
       it is advised to use the optimization flag -o.   This  flag  tells  i.atcorr  to  try  and
       speedup calculations.  However, this option will increase memory requirements.

       If flag -r is used, the input raster data are treated as reflectance. Otherwise, the input
       raster data are treated as radiance  values  and  are  converted  to  reflectance  at  the
       i.atcorr runtime. The output data are always reflectance.

       Note that the satellite overpass time has to be specified in Greenwich Mean Time (GMT).

       An example 6S parameters:
       8                            - geometrical conditions=Landsat ETM+
       2 19 13.00 -47.410 -20.234   - month day hh.ddd longitude latitude ("hh.ddd" is in decimal
       hours GMT)
       1                            - atmospheric mode=tropical
       1                            - aerosols model=continental
       15                           - visibility [km] (aerosol model concentration)
       -0.600                       - mean target elevation above sea level [km] (here 600m asl)
       -1000                        - sensor height (here, sensor on board a satellite)
       64                           - 4th band of ETM+ Landsat 7
        If the position is not available in longitude-latitude  (WGS84),  the  m.proj  conversion
       module can be used to reproject from a different projection.

6S CODE PARAMETER CHOICES

   A. Geometrical conditions
       | Code | Description | Details
       | 1 | meteosat observation | enter month,day,decimal hour (universal time-hh.ddd)

       n. of column,n. of line. (full scale 5000*2500)
       | 2 | goes east observation | enter month,day,decimal hour (universal time-hh.ddd)

       n. of column,n. of line. (full scale 17000*12000)c
       | 3 | goes west observation | enter month,day,decimal hour (universal time-hh.ddd)

       n. of column,n. of line. (full scale 17000*12000)
       | 4 | avhrr (PM noaa) | enter month,day,decimal hour (universal time-hh.ddd)

       n. of column(1-2048),xlonan,hna

       give long.(xlonan) and overpass hour (hna) at

       the ascendant node at equator
       | 5 | avhrr (AM noaa) | enter month,day,decimal hour (universal time-hh.ddd)

       n. of column(1-2048),xlonan,hna

       give long.(xlonan) and overpass hour (hna) at

       the ascendant node at equator
       | 6 | hrv (spot) | enter month,day,hh.ddd,long.,lat. *
       | 7 | tm (landsat) | enter month,day,hh.ddd,long.,lat. *
       | 8 | etm+ (landsat7) | enter month,day,hh.ddd,long.,lat. *
       | 9 | liss (IRS 1C) | enter month,day,hh.ddd,long.,lat. *
       | 10 | aster | enter month,day,hh.ddd,long.,lat. *
       | 11 | avnir | enter month,day,hh.ddd,long.,lat. *
       | 12 | ikonos | enter month,day,hh.ddd,long.,lat. *
       | 13 | RapidEye | enter month,day,hh.ddd,long.,lat. *
       | 14 | VGT1 (SPOT4) | enter month,day,hh.ddd,long.,lat. *
       |  15 | VGT2 (SPOT5) | enter month,day,hh.ddd,long.,lat. * * NOTE: for HRV, TM, ETM+, LISS
       and ASTER experiments, longitude and latitude are the coordinates  of  the  scene  center.
       Latitude  must  be > 0 for northern hemisphere and < 0 for southern. Longitude must be > 0
       for eastern hemisphere and < 0 for western.

   B. Atmospheric model
       | Code | Meaning
       | 0 | no gaseous absorption
       | 1 | tropical
       | 2 | midlatitude summer
       | 3 | midlatitude winter
       | 4 | subarctic summer
       | 5 | subarctic winter
       | 6 | us standard 62
       | 7 | Define your own atmospheric model as a set of the following 5  parameters  per  each
       measurement:
       altitude [km]
       pressure [mb]
       temperature [k]
       h2o density [g/m3]
       o3 density [g/m3]
       For  example: there is one radiosonde measurement for each altitude of 0-25km at a step of
       1km, one measurment for each altitude of  25-50km  at  a  step  of  5km,  and  two  single
       measurements  for altitudes 70km and 100km. This makes 34 measurments. In that case, there
       are 34*5 values to input.
       | 8 | Define your own atmospheric model providing values of  the  water  vapor  and  ozone
       content:
       uw [g/cm2]
       uo3 [cm-atm]
        The profile is taken from us62.

   C. Aerosols model
       | Code | Meaning | Details
       | 0 | no aerosols |
       | 1 | continental model |
       | 2 | maritime model |
       | 3 | urban model |
       | 4 | shettle model for background desert aerosol |
       | 5 | biomass burning |
       | 6 | stratospheric model |
       | 7 | define your own model | Enter the volumic percentage of each component:
       c(1) = volumic % of dust-like
       c(2) = volumic % of water-soluble
       c(3) = volumic % of oceanic
       c(4) = volumic % of soot
       All values between 0 and 1.
       |  8  | define your own model | Size distribution function: Multimodal Log Normal (up to 4
       modes).
       | 9 | define your own model | Size distribution function: Modified gamma.
       | 10 | define your own model | Size distribution function: Junge Power-Law.
       | 11 | define your own model | Sun-photometer measurements, 50 values max, entered as:
       r and d V / d (logr)
       where r is the radius [micron], V is the volume, d V / d (logr) [cm3/cm2/micron].
       Followed by:
       nr and ni for each wavelength
       where nr and ni are respectively the real and imaginary part of the refractive index.

   D. Aerosol concentration model (visibility)
       If you have an estimate of the meteorological parameter visibility v, enter  directly  the
       value  of v [km] (the aerosol optical depth (AOD) will be computed from a standard aerosol
       profile).

       If you have an estimate of aerosol optical depth, enter 0 for  the  visibility  and  in  a
       following  line  enter  the  aerosol  optical depth at 550nm (iaer means 'i' for input and
       'aer' for aerosol), for example:

       0                            - visibility
       0.112                        - aerosol optical depth 550 nm

       NOTE: if iaer is 0, enter -1 for visibility.

   E. Target altitude (xps), sensor platform (xpp)
       Target altitude (xps, in negative [km]): xps >= 0 means the target is at the sea level.
       otherwise xps expresses the altitude of the target (e.g., mean elevation) in  [km],  given
       as negative value

       Sensor platform (xpp, in negative [km] or -1000):
       xpp = -1000 means that the sensor is on board a satellite.
       xpp = 0 means that the sensor is at the ground level.
       -100  <  xpp  <  0  defines the altitude of the sensor expressed in [km]; this altitude is
       given relative to the target altitude as negative value.

       For aircraft simulations only (xpp is neither equal to 0  nor  equal  to  -1000):  puw,po3
       (water vapor content,ozone content between the aircraft and the surface)
       taerp (the aerosol optical thickness at 550nm between the aircraft and the surface)

       If these data are not available, enter negative values for all of them.  puw,po3 will then
       be interpolated from the us62 standard profile according  to  the  values  at  the  ground
       level. taerp will be computed according to a 2km exponential profile for aerosol.

   F. Sensor band
       There  are two possibilities: either define your own spectral conditions (codes -2, -1, 0,
       or 1) or choose a code indicating the band of one of the pre-defined satellites.

       Define your own spectral conditions:
       | Code | Meaning
       | -2 | Enter wlinf, wlsup.
       The filter function will be equal to 1 over the whole band (as iwave=0) but step  by  step
       output will be printed.
       | -1 | Enter wl (monochr. cond, gaseous absorption is included).
       | 0 | Enter wlinf, wlsup.
       The filter function will be equal to 1over the whole band.
       |  1  |  Enter  wlinf,  wlsup  and  user's  filter  function  s(lambda)  by step of 0.0025
       micrometer.

       Pre-defined satellite bands:
            | Code    | Meaning
            | 2  | meteosat vis band (0.350-1.110)
            | 3  | goes east band vis (0.490-0.900)
            | 4  | goes west band vis (0.490-0.900)
            | 5  | avhrr (noaa6) band 1 (0.550-0.750)
            | 6  | avhrr (noaa6) band 2 (0.690-1.120)
            | 7  | avhrr (noaa7) band 1 (0.500-0.800)
            | 8  | avhrr (noaa7) band 2 (0.640-1.170)
            | 9  | avhrr (noaa8) band 1 (0.540-1.010)
            | 10 | avhrr (noaa8) band 2 (0.680-1.120)
            | 11 | avhrr (noaa9) band 1 (0.530-0.810)
            | 12 | avhrr (noaa9) band 1 (0.680-1.170)
            | 13 | avhrr (noaa10) band 1 (0.530-0.780)
            | 14 | avhrr (noaa10) band 2 (0.600-1.190)
            | 15 | avhrr (noaa11) band 1 (0.540-0.820)
            | 16 | avhrr (noaa11) band 2 (0.600-1.120)
            | 17 | hrv1 (spot1) band 1 (0.470-0.650)
            | 18 | hrv1 (spot1) band 2 (0.600-0.720)
            | 19 | hrv1 (spot1) band 3 (0.730-0.930)
            | 20 | hrv1 (spot1) band pan (0.470-0.790)
            | 21 | hrv2 (spot1) band 1 (0.470-0.650)
            | 22 | hrv2 (spot1) band 2 (0.590-0.730)
            | 23 | hrv2 (spot1) band 3 (0.740-0.940)
            | 24 | hrv2 (spot1) band pan (0.470-0.790)
            | 25 | tm (landsat5) band 1 (0.430-0.560)
            | 26 | tm (landsat5) band 2 (0.500-0.650)
            | 27 | tm (landsat5) band 3 (0.580-0.740)
            | 28 | tm (landsat5) band 4 (0.730-0.950)
            | 29 | tm (landsat5) band 5 (1.5025-1.890)
            | 30 | tm (landsat5) band 7 (1.950-2.410)
            | 31 | mss (landsat5) band 1 (0.475-0.640)
            | 32 | mss (landsat5) band 2 (0.580-0.750)
            | 33 | mss (landsat5) band 3 (0.655-0.855)
            | 34 | mss (landsat5) band 4 (0.785-1.100)
            | 35 | MAS (ER2) band 1 (0.5025-0.5875)
            | 36 | MAS (ER2) band 2 (0.6075-0.7000)
            | 37 | MAS (ER2) band 3 (0.8300-0.9125)
            | 38 | MAS (ER2) band 4 (0.9000-0.9975)
            | 39 | MAS (ER2) band 5 (1.8200-1.9575)
            | 40 | MAS (ER2) band 6 (2.0950-2.1925)
            | 41 | MAS (ER2) band 7 (3.5800-3.8700)
            | 42 | MODIS band 1 (0.6100-0.6850)
            | 43 | MODIS band 2 (0.8200-0.9025)
            | 44 | MODIS band 3 (0.4500-0.4825)
            | 45 | MODIS band 4 (0.5400-0.5700)
            | 46 | MODIS band 5 (1.2150-1.2700)
            | 47 | MODIS band 6 (1.6000-1.6650)
            | 48 | MODIS band 7 (2.0575-2.1825)
            | 49 | avhrr (noaa12) band 1 (0.500-1.000)
            | 50 | avhrr (noaa12) band 2 (0.650-1.120)
            | 51 | avhrr (noaa14) band 1 (0.500-1.110)
            | 52 | avhrr (noaa14) band 2 (0.680-1.100)
            | 53 | POLDER band 1 (0.4125-0.4775)
            | 54 | POLDER band 2 (non polar) (0.4100-0.5225)
            | 55 | POLDER band 3 (non polar) (0.5325-0.5950)
            | 56 | POLDER band 4 P1 (0.6300-0.7025)
            | 57 | POLDER band 5 (non polar) (0.7450-0.7800)
            | 58 | POLDER band 6 (non polar) (0.7000-0.8300)
            | 59 | POLDER band 7 P1 (0.8100-0.9200)
            | 60 | POLDER band 8 (non polar) (0.8650-0.9400)
            | 61 | etm+ (landsat7) band 1 (0.435-0.520)
            | 62 | etm+ (landsat7) band 2 (0.506-0.621)
            | 63 | etm+ (landsat7) band 3 (0.622-0.702)
            | 64 | etm+ (landsat7) band 4 (0.751-0.911)
            | 65 | etm+ (landsat7) band 5 (1.512-1.792)
            | 66 | etm+ (landsat7) band 7 (2.020-2.380)
            | 67 | etm+ (landsat7) band 8 (0.504-0.909)
            | 68 | liss (IRC 1C) band 2 (0.502-0.620)
            | 69 | liss (IRC 1C) band 3 (0.612-0.700)
            | 70 | liss (IRC 1C) band 4 (0.752-0.880)
            | 71 | liss (IRC 1C) band 5 (1.452-1.760)
            | 72 | aster  band 1 (0.480-0.645)
            | 73 | aster band 2 (0.588-0.733)
            | 74 | aster band 3N (0.723-0.913)
            | 75 | aster band 4 (1.530-1.750)
            | 76 | aster band 5 (2.103-2.285)
            | 77 | aster band 6 (2.105-2.298)
            | 78 | aster band 7 (2.200-2.393)
            | 79 | aster band 8 (2.248-2.475)
            | 80 | aster band 9 (2.295-2.538)
            | 81 | avnir band 1 (0.390-0.550)
            | 82 | avnir band 2 (0.485-0.695)
            | 83 | avnir band 3 (0.545-0.745)
            | 84 | avnir band 4 (0.700-0.925)
            | 85 | ikonos Green band (0.350-1.035)
            | 86 | ikonos Red band (0.350-1.035)
            | 87 | ikonos NIR band (0.350-1.035)
            | 88 | RapidEye Blue band (0.438-0.513)
            | 89 | RapidEye Green band (0.463-0.594)
            | 90 | RapidEye Red band (0.624-0.690)
            | 91 | RapidEye RedEdge band (0.500-0.737)
            | 92 | RapidEye NIR band (0.520-0.862)
            | 93 | VGT1 (SPOT4) band 0 (0.400-0.500)
            | 94 | VGT1 (SPOT4) band 2 (0.580-0.782)
            | 95 | VGT1 (SPOT4) band 3 (0.700-1.030)
            | 96 | VGT1 (SPOT4) MIR band (1.450-1.800)
            | 97 | VGT2 (SPOT5) band 0 (0.400-0.550)
            | 98 | VGT2 (SPOT5) band 2 (0.580-0.780)
            | 99 | VGT2 (SPOT5) band 3 (0.700-1.000)
            | 100     | VGT2 (SPOT5) MIR band (1.450-1.800)

EXAMPLES

   Atmospheric correction of a LANDSAT-7 channel
       The example is based on the North Carolina sample dataset (GMT -5 hours).   First  we  set
       the computational region to the satellite map, e.g. channel 4:
       g.region rast=lsat7_2002_40 -p
         It  is  important  to verify the available metadata for the sun position which has to be
       defined for the atmospheric correction. An option is to check the satellite overpass  time
       with  sun  position  as  reported in metadata. For the North Carolina sample dataset, they
       have also been stored for each channel and can be retrieved like this:
       r.info lsat7_2002_40
        In this case, we have: SUN_AZIMUTH = 120.8810347, SUN_ELEVATION = 64.7730999.

       If the sun position metadata are unavailable, we can also calculate them from the overpass
       time as follows (r.sunmask uses SOLPOS):
       r.sunmask  -s  elev=elevation  out=dummy  year=2002  month=5  day=24  hour=10 min=42 sec=7
       timezone=-5
       # .. reports: sun  azimuth:  121.342461,  sun  angle  above  horz.(refraction  corrected):
       65.396652
        If the overpass time is unknown, use the Satellite Overpass Predictor.

       Convert  DN  (digital number = pixel values) to Radiance at top-of-atmosphere (TOA), using
       the formula
          L&lambda;  =  ((LMAX&lambda;  -  LMIN&lambda;)/(QCALMAX-QCALMIN))  *  (QCAL-QCALMIN)  +
       LMIN&lambda;
        Where:

                      L&lambda;  =  Spectral  Radiance  at  the  sensor's aperture in Watt/(meter
                     squared * ster * &micro;m), the apparent radiance as seen by  the  satellite
                     sensor;

                      QCAL = the quantized calibrated pixel value in DN;

                      LMIN&lambda;  =  the  spectral  radiance  that  is  scaled  to  QCALMIN  in
                     watts/(meter squared * ster * &micro;m);

                      LMAX&lambda;  =  the  spectral  radiance  that  is  scaled  to  QCALMAX  in
                     watts/(meter squared * ster * &micro;m);

                      QCALMIN  =  the  minimum quantized calibrated pixel value (corresponding to
                     LMIN&lambda;) in DN;

                      QCALMAX = the maximum quantized calibrated pixel  value  (corresponding  to
                     LMAX&lambda;) in DN=255.
       LMIN&lambda;  and  LMAX&lambda;  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 (see Landsat handbook), and  the  maximal  DN
       value  (QCALMAX)  is  255.  QCAL  is  the DN value for every separate pixel in the Landsat
       image.

       We extract the coefficients and apply them in order to obtain the radiance map:
       CHAN=4
       r.info lsat7_2002_${CHAN}0 -h | tr '\n' '  '  |  sed  's+  ++g'  |  tr  ':'  '\n'  |  grep
       "LMIN_BAND${CHAN}\|LMAX_BAND${CHAN}"
       LMAX_BAND4=241.100,p016r035_7x20020524.met
       LMIN_BAND4=-5.100,p016r035_7x20020524.met
       QCALMAX_BAND4=255.0,p016r035_7x20020524.met
       QCALMIN_BAND4=1.0,p016r035_7x20020524.met
         Conversion  to  radiance  (this calculation is done for band 4, for the other bands, the
       numbers in italics need to be replaced with their related values):
       r.mapcalc "lsat7_2002_40_rad=((241.1 - (-5.1)) / (255.0 - 1.0)) * (lsat7_2002_40 - 1.0)  +
       (-5.1)"

       #  find  mean  elevation  (target above sea level, used as initialization value in control
       file)
       r.univar elevation
        Create a control file 'icnd.txt' for channel 4 (NIR), based on metadata. For the overpass
       time, we need to define decimal hours:
       10:42:07  NC  local  time  = 10.70 decimal hours (decimal minutes: 42 * 100 / 60) which is
       15.70 GMT:
       8                            - geometrical conditions=Landsat ETM+
       5 24 15.70 -78.691 35.749    - month day hh.ddd longitude latitude  ("hh.ddd"  is  in  GMT
       decimal hours)
       2                            - atmospheric mode=midlatitude summer
       1                            - aerosols model=continental
       50                           - visibility [km] (aerosol model concentration)
       -0.110                       - mean target elevation above sea level [km]
       -1000                        - sensor on board a satellite
       64                           - 4th band of ETM+ Landsat 7
         Finally, run the atmospheric correction (-r for reflectance input map; -a for date >July
       2000; -o to use cache acceleration):
       i.atcorr    -r    -a    -o    lsat7_2002_40_rad     ialt=elevation     icnd=icnd_lsat4.txt
       oimg=lsat7_2002_40_atcorr
         Note  that  the  altitude  value  from 'icnd_lsat4.txt' file is read at the beginning to
       compute the initial transform. It is necessary to give a value which  could  be  the  mean
       value  of  the  elevation  model. For the atmospheric correction then the raster elevation
       values are used from the map.

       Note that the process is computationally intensive.
       Note also, that i.atcorr reports solar elevation angle above  horizon  rather  than  solar
       zenith angle.

REMAINING DOCUMENTATION ISSUES

       1.  The  influence and importance of the visibility value or map should be explained, also
       how to obtain an estimate for either visibility or aerosol optical depth at 550nm.

SEE ALSO

       GRASS Wiki page about Atmospheric correction

        r.info, r.mapcalc, r.univar

REFERENCES

                      Vermote, E.F., Tanre, D., Deuze, J.L., Herman,  M.,  and  Morcrette,  J.J.,
                     1997,  Second  simulation of the satellite signal in the solar spectrum, 6S:
                     An overview., IEEE Trans. Geosc. and Remote Sens. 35(3):675-686.

                      6S Manual: PDF1, PDF2, and PDF3

                     RapidEye sensors have been provided by RapidEye AG, Germany

AUTHORS

       Original version of the program for GRASS 5:
       Christo Zietsman, 13422863(at)sun.ac.za

       Code clean-up and port to GRASS 6.3, 15.12.2006:
       Yann Chemin, ychemin(at)gmail.com

       Documentation clean-up + IRS LISS sensor addition 5/2009:
       Markus Neteler, FEM, Italy

       ASTER sensor addition 7/2009:
       Michael Perdue, Canada

       AVNIR, IKONOS sensors addition 7/2010:
       Daniel Victoria, Anne Ghisla

       RapidEye sensors addition 11/2010:
       Peter L&ouml;we, Anne Ghisla

       VGT1 and VGT2 sensors addition from 6SV-1.1 sources, addition 07/2011:
       Alfredo Alessandrini, Anne Ghisla

       Last changed: $Date: 2013-02-15 14:04:18 -0800 (Fri, 15 Feb 2013) $

       Full index

       © 2003-2013 GRASS Development Team