Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have been using Erdas ATCOR for atmospheric correction of ASTER imagery. What opensource options exist for atmospheric correction of ASTER Level 1B products? I'm primarily interested in converting registered radiance at the sensor to true reflectance values.

share|improve this question

1 Answer 1

up vote 5 down vote accepted

One open source option for atmospherically correcting ASTER L1B products, in order to convert at-sensor Radiance values to Top of Canopy Reflectances, is GRASS GIS' i.atcorr module.

An implementation of the 6S algorithm in GRASS GIS

GRASS GIS features a dedicated module for the task in question called i.atcorr (in GRASS-GIS version 7 or in GRASS GIS vesion 6.4.x). The module is an implementation of NASA's 6S algorithm (Second Simulation of Satellite Signal in the Solar Spectrum). As stated in Land Surface Reflectance Science Computing Facility's web interface to the code:

The 6s code predicts the satellite signal from 0.25 to 4.0 microns assuming cloudless atmosphere. The main atmospheric effects (gaseous absorption by water vapor,carbon dioxyde,oxygen and ozone; scattering by molecules and aerosols) are taken into account. Non-uniform surfaces may be considered, as well as bidirectional reflectances as boundary conditions.

Using i.atcorr

In order to use GRASS' module, one needs to prepare a simple text file that contains the required parameters for the algorithm to run for each spectral band separately. Specifically, the text file should be constructed so as to describe in one line (for) each of the following parameters: Geometrical Conditions, Atmospherical Model, Aerosols Model, Visibility, Target & Sensor Altitude, Sensor's Band (spectral conditions). Among the various sensor's pre-described spectral conditions, the module contains the necessary relative spectral response data for ASTER acquisitions.

Such a parameters description file (here the example concerns Landsat data) might look like:

7
11 22 8.83574 23.83895740 34.59989709
3
1
0
0.111
-0.500
-1000
25

or even with comments such as

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

All of the above are explained in the module's manual. Elaborating a bit more on how to increase the module's qualitative performance, a visibility estimation can be defined. In the absence of visibility maps, an estimation of the Aerosol Optical Depth at 550 nm can be derived from MODIS' MOD08 - Gridded Atmospheric Product. This products includes 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. Some notes in the relevant GRASS-Wiki subsection Sources for AOD estimations.

Scripting

Repeating the exercise of setting up a parameters file for each band separately is naturally a tedious task. Scripting is recommended. A custom (and messy, containing hard-coded values) bash-shell script might look like:

#!/bin/bash
# set -u

# Nikos Alexandris, February 2013
# based on a script provided by Yann Chemin!
# i.atcorr scripting

# Warning: still contains hardcoded stuff!

echo "Usage: $0 [Meant Target Elevation] [AOD]"
echo "Note, the script has to be eXecuted from the directory that contains"
echo "the *MTL.txt metadata and within from the GRASS environment!"

# first, remove access to all mapsets but the current
g.mapsets mapset=`g.mapsets -p sep=newline | grep -v `g.mapset -p` | tr "\n" "," | sed 's:^\(.*\)\,$:\1:'` operation=remove
echo "Access to all mapsets removed..."
echo "Please, grant access to mapsets of interest and rerun the script!"

# # add all mapsets containing LT5 imagery?
# g.mapsets addmapset=`g.mapsets -l --v fs=, | tr -d "\n"` # grass64
# g.mapsets removemapset=`g.mapsets -l --v fs=, | tr -d "\n"` # grass70

# access only to specific mapsets! - currently, by hand!
# g.mapsets mapset=`g.mapsets -l --v sep=comma` operator=remove
Scenes_To_Correct=$(g.mapsets -p)
echo "Correcting scenes: $Scenes_To_Correct"
sleep 1 # wait a bit!

echo 
echo "-------------------------------------------------"
echo "6S Atmospheric Correction algorithm)"
echo "-------------------------------------------------"
# # echo "It will create for *.surf.* images from i.atcorr" 
# # echo "Atmospherically corrected image=surface reflectance"

echo
echo "The following parameters are required"
echo "---------------------------------------------------"
echo

# Location of parameter file(s)
Parameters_Pool="/geo/geodata/imagery/ellas/landsat/icnd_pool"
echo
echo "[The directory containing parameter files: $Parameters_Pool]"
sleep 1 # wait a bit!

# Define fixed parameters

Geometry=7 # Geometrical conditions (L5TM)
echo "Geometrical conditions index for Landsat TM 5: $Geometry"

Sensor_Height=-1000 # for satellite borne sensors set to -1000
echo "The sensor's height:          $Sensor_Height"

Atmospheric_Mode=3 # here: midlatitude winter
echo "Atmospheric mode:             $Atmospheric_Mode"

Aerosol_Model=1 # here: continental
echo "Aerosol model:                $Aerosol_Model"

# Define non-fixed parameters

# get first parameter $1 or read MTE!
        if [ ! -z $1 ]; # if $1 is not null, set MeanTargetElevation
          then MeanTargetElevation=$1
          else echo "Please provide a Mean Target Elevation (negatively signed value in Km, e.g. -0.5)"
          read MeanTargetElevation
          echo
          echo "Note, this value will be overwritten if a DEM raster has been defined as an input!"
        fi

echo "Mean target elevation parameter:       $MeanTargetElevation Km"

# aerosol depth, summer vs winter
  if [ ! -z $2 ];
    then AOD=$2
    echo "The provided aerosol optical depth (AOD) is: ${AOD} (unitless)"
    else echo "Please provide an AOD estimation"
    read AOD # else, read AOD
  fi

# AOD_Winter=0.111
# AOD_Summer=0.222
# AOD="${AOD_Winter}" # set to winter AOD
#   if [ ${Month} -gt 4 ] && [ ${Month} -lt 10 ]; # compare month of acquisition
#       then AOD="${AOD_Summer}"  # set to summer AOD if...
#   fi
echo "The provided aerosol optical depth: ${AOD} (unitless)"

# if AOD provided, set Visibility to Zero
  if [ ! -z ${AOD} ];
    then Visibility=0
    echo 
    echo "Visibility was set to zero since AOD = ${AOD} is provided!"
    else echo "Please provide a visibility value:"
    read Visibility
    echo
  fi

        # ToDo: elaborate on Visibillity -- lines below from YannC's code
# # vis_list=(10 10 8 9.7 15 8 7 10 10 9.7 12 9.7 7 12 12 12 3 15 12 9.7 6 15)
# # vis_len=${#vis_list[*]}
# # echo $vis_len
# # i=0

# spectral conditions index for L5TM:
# 25,26,27,28,29,30 respectively for bands 1, 2, 3, 4, 5, 7
Satellite_Band_No=25 # 1st L5TM band to undergo atmospheric correction -- counter!

for SCENE in $Scenes_To_Correct
  do

# Here we suppose you have Visibility (Visibility) maps ready

# #  r.mapcalc expression="visibility=${vis_list[$i]}" --overwrite

# #  # Dummy visibility value for atcorr param file
# #  vis=12 
# #  #Increment i
# #  i=$(echo "$i + 1" | bc)

echo "Acquisition metadata"
echo

# scene's basename as in GRASS' db
Basename=$(g.mapset -p) #LT5_Basename=$(basename $BAND _MTL.txt)
echo "The scene's base name is ${Basename}"

# Date and time of the satellite's overpass (month, day, GMT in decimal hours)

# grep from the filename?
# MonthDay=$(echo $L5_Basename | sed 's/\(.*\)_.......\(.*\)/\2/')

# Safer to use metadata content!

# grep acquisition date from the metadata file
Date_Acquired=$(grep -a DATE_ACQUIRED "${Basename}"_MTL.txt | cut -d" " -f7)

Month=$(echo "${Date_Acquired}" | cut -d"-" -f2) # grep Month
#echo "Month of acquisition $Month"

Day=$(echo "${Date_Acquired}" | cut -d"-" -f3) # grep Day
#echo "Day of acquisition $Day"

# GMT in decimal hours
GMT=$(grep -a -e "TIME" "${Basename}"_MTL.txt | cut -d" " -f7 | sed 's:^\(.*\)Z$:\1:' | awk -F: '{print $1 + ( $2 / 60) + ($3 / 3600) }')
# echo "UTC time of the acquisition: $GMT"

# Merge in one variable
MDH="$Month $Day $GMT"
echo "Month, Day and UTC Decimal Hours of the acquisition: ${MDH}"


#Get the scene's central Lat/Long
Long_NonProj=`g.region -c | grep easting | sed 's/center\ easting:\ \(.*\)/\1/'`
Lat_NonProj=`g.region -c | grep northing | sed 's/center\ northing:\ \(.*\)/\1/'`
Long=`g.region -cgl | grep center_long | sed 's/center_long=\(.*\)/\1/'`
Lat=`g.region -cgl | grep center_lat | sed 's/center_lat=\(.*\)/\1/'`

echo "The scene's center geographic coordinates (in decimal degrees): ${Long_NonProj} $Lat_NonProj ($Long $Lat)"

# loop over Landsat bands in question
for Band_No in 1 2 3 4 5 7
  do # Generate the parameterization file (icnd_landsat5)
  echo
  echo "Processing band $Band_No (spectral conditions index: $Satellite_Band_No)"

  #   echo "$Geometry - geometrical conditions=Landsat 5 TM" > $Parameters_Pool/${Basename}_icnd_${Band_No}.txt
  echo $Geometry > $Parameters_Pool/${Basename}_icnd_${Band_No}.txt

  #   echo "$MDH $Long $Lat - month day hh.ddd longitude latitude" >>  $Parameters_Pool/${Basename}_icnd_${Band_No}.txt
  echo $MDH $Long $Lat >>  $Parameters_Pool/${Basename}_icnd_${Band_No}.txt

  #   echo "$Atmospheric_Mode - atmospheric mode=tropical" >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt
  echo $Atmospheric_Mode >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt

  #   echo "$Aerosol_Model - aerosols model=continental" >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt
  echo $Aerosol_Model >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt

  #   echo "$Visibility - visibility [km] (aerosol model concentration), not used as there is raster input" >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt
  echo $Visibility >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt

  #   echo "$AOD - aerosol optical depth at 550nm" >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt
  echo $AOD >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt

  # echo "$MeanTargetElevation - mean target elevation above sea level [km] (here 600m asl), not used as there is raster input" >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt
  echo $MeanTargetElevation >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt

  #   echo "$Sensor_Height - sensor height (here, sensor on board a satellite)" >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt
  echo $Sensor_Height >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt

  #   echo "$Satellite_Band_No - band ${Band_No} of TM Landsat 5" >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt
  echo $Satellite_Band_No >> $Parameters_Pool/${Basename}_icnd_${Band_No}.txt

  # Process band-wise atmospheric correction with 6s
  echo "Using the following parameters:"
  cat $Parameters_Pool/${Basename}_icnd_${Band_No}.txt | tr "\n" " "

# echo "Executing:"
# echo "i.atcorr -r \ "
# echo "input=Basename.ToAR.$Band_No \ "
# echo "elevation=$MeanTargetElevation visibility=$Visibility \ "
# echo "parameters=$Parameters_Pool/${Basename}_icnd_${Band_No}.txt \ "
# echo "output=B.Trimmed.ToCR.$Band_No \ "
# echo "range=0,1 rescale=0,1 \ "
# echo "--overwrite"

# Attention!
  ## input is radiance or reflectance?
  ## flag "-r" ?
  ## does the elevation cover the area of the images?
  ## input and output is converted in to float values anyway!

# i.atcorr -r --overwrite \
# input=B.Trimmed.ToAR.$Band_No \
# range=0,1 \
# elevation=aster_gdem2_ellas_smooothed_threshold_08_iterations_10 \
# parameters=$Parameters_Pool/${Basename}_icnd_${Band_No}.txt \
# output=Test_$Basename_$Band_No \
# rescale=0,1

# without elevation here!
i.atcorr -r --overwrite --v \
input=B.Trimmed.DOS1.ToCR.$Band_No \
range=0,1 \
parameters=${Parameters_Pool}/${Basename}_icnd_${Band_No}.txt \
output=B.Trimmed.ToCR.${Band_No} \
rescale=0,1

r.info -r B.Trimmed.ToCR.${Band_No} 
sleep 1

Satellite_Band_No=$(($Satellite_Band_No+1))
  if [ $Satellite_Band_No -eq 30 ]
   then break
   elif [ $Satellite_Band_No -gt 30 ]
   then Satellite_Band_No=25
  fi

 done
done

Sources

the Algorithm

in GRASS GIS:

in GRASS-Wiki:

share|improve this answer
1  
A Python script would be more attractive as it would be (practically) cross-platform. –  Nikos Alexandris Oct 12 '13 at 10:26
1  
Many thanks, this is a good alternative to Erdas. –  Aaron Oct 16 '13 at 22:00
1  
@Aaron Also, the process to add new sensors is not that hard though it takes some time and requires, of course, to be checked by more than one person so as to ensure quality. On the tracks are the high-resolution sensors QuickBird-2 and WorldView-2 -- hopefully they will be added officially soon (there is a ticket already for those). Hopefully the module will receive further improvements, e.g. in what concerns retrieving various metadata (at least for easily accessible/public domain data such as Landsat and MODIS imagery). –  Nikos Alexandris Oct 16 '13 at 23:04

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.