Alternativas de código abierto a ERDAS IMAGINE ATCOR?

Question

He estado con ERDAS IMAGINE ATCOR de extensión para la Corrección Atmosférica de las imágenes ASTER.

Lo opensource que existen opciones para la corrección atmosférica de ASTER Nivel 1B productos?

Estoy interesado principalmente en la conversión registrado resplandor en el sensor a los verdaderos valores de reflectancia.

Answer 1

Una opción de código abierto para atmosférico corrección de ASTER L1B productos, con el fin de convertir al sensor de Luminosidad valores a la parte Superior de la copa Reflectancias, es GRASS GIS' yo.atcorr módulo.

Una aplicación de las 6S algoritmo en GRASS GIS

GRASS GIS cuenta con un módulo específico para la tarea en cuestión se llama i.atcorr (en la HIERBA-SIG versión 7 o en GRASS GIS vesión 6.4.x). El módulo es una implementación de la NASA 6S algoritmo (Segunda Simulación de la Señal de Satélite en el Espectro Solar). Como se indica en la Superficie de la Tierra de Reflectancia de Ciencia de Computación de la Instalación's interfaz web para el código:

El 6s código predice la señal del satélite de 0.25 a 4.0 micrones suponiendo sin nubes de la atmósfera. Los principales efectos atmosféricos (en estado gaseoso la absorción por el vapor de agua,dioxyde de carbono,el oxígeno y el ozono; la dispersión de por moléculas y aerosoles) son tomadas en cuenta. No uniforme las superficies pueden ser considerados, así como bidireccional reflectancias como las condiciones de contorno.

Usando yo.atcorr

Con el fin de utilizar la HIERBA de módulo, uno tiene que preparar un simple archivo de texto que contiene los parámetros necesarios para el algoritmo a ejecutar para cada banda espectral por separado. Específicamente, el archivo de texto debe estar construido de manera que describir en una sola línea (por) cada uno de los siguientes parámetros: Condiciones Geométricas, Meteorológicas Modelo, Aerosoles Modelo, la Visibilidad, el Objetivo Y el Sensor de Altitud, el Sensor de la Banda (espectro de condiciones). Entre los diversos sensor de pre-describe espectral condiciones, el módulo contiene la necesaria relativa de la respuesta espectral de los datos para ASTER adquisiciones.

Tal descripción de los parámetros del archivo (aquí el ejemplo se refiere a los datos de Landsat) aspecto:

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

o incluso con comentarios tales como

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

Todo lo anterior se explicó en el módulo del manual. Elaborar un poco más sobre cómo aumentar el módulo cualitativo en el rendimiento, una visibilidad de estimación puede ser definido. En la ausencia de visibilidad de los mapas, una estimación del espesor Óptico de Aerosoles a 550 nm puede ser derivado de MODIS' MOD08 - Cuadriculado Atmosférica Producto. Este producto incluye la capa Optical_Depth_Land_And_Ocean_Mean, descrito como el Espesor Óptico del Aerosol en el 0,55 micrones para ambos Océanos (el mejor) y la Tierra (corregido): la Media. Algunas notas relevantes de la HIERBA-Wiki inciso Fuentes de AOD estimaciones.

De secuencias de comandos

Repetir el ejercicio de la configuración de un archivo de parámetros para cada banda por separado es, naturalmente, una tarea tediosa. Las secuencias de comandos está recomendado.

en Bash

Una costumbre (y desordenado, que contiene valores codificados) de bash shell script podría ser así:

#!/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

en Python

Una secuencia de comandos de Python "trabajo en progreso" y sujeto a cambios, tratando de automatizar el uso de i.atcorr de imágenes Landsat está disponible aquí estoy.landsat.atcorr. Funciona bien hasta ahora para Landsat OLI, ETM+ y TM sensores.

Fuentes

el Algoritmo

• 6S página web
• El 6S Manual
• Publicaciones

en GRASS GIS:

• yo.atcorr's manuales de ver. 6.4.x, ver. 7

en PASTO-Wiki:

• Corrección atmosférica
• Varias referencias sobre el tema