I have been using ERDAS IMAGINE's ATCOR extension 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.
Answer
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.
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.
Repeating the exercise of setting up a parameters file for each band separately is naturally a tedious task. Scripting is recommended.
in Bash
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
in Python
A Python script "work-in-progress" & subject to changes, attempting to automate the use of i.atcorr for Landsat imagery is available here i.landsat.atcorr. It works well so far for Landsat OLI, ETM+ and TM sensors.
the Algorithm
in GRASS GIS:
in GRASS-Wiki:
No comments:
Post a Comment