[GRASS-user] Re: i.atcorr returns all NULL values

[Chemin, Yann (IWMI)]

p192 r030 image of July 2003 of Italy (L5TM) successfully corrected for
band 1 and 2. 
The following script works well in GRASS GIS Trunk SVN.
Please try.

#!/bin/bash

echo "Part 0 (Import in GRASS GIS)"
echo "RUN from the MTL.txt directory and within the GRASS environment"
echo "-----------------------------------"
echo "It will create for *.[1-7] images from r.in.gdal"
echo "image=DN"

# DEM file name
dem=dem
r.mapcalc expression=dem=5.0
# r.in.gdal input=$dem output=$dem


for file in L5*.TIF
do
       out=$(echo $file | sed
's/\(.*\)_\(.*\)_B\(.*\)0.TIF/\1\_\2\.\3/g')
       echo $out
       r.in.gdal --overwrite input=$file output=$out
done

echo "Part 1 (After DN)"
echo "-----------------------------------"
echo "It will create for *.toar.* images from i.landsat.toar"
echo "image=top of atmosphere reflectance"

for L5_MTL_file in L5*_MTL.txt
do
       L5_prefix=$(echo $L5_MTL_file | sed 's/\(.*\)_MTL.txt/\1/')
       i.landsat.toar -t input_prefix=$L5_prefix\.
output_prefix=$L5_prefix\.toar.  metfile=$L5_MTL_file sensor=tm5
done


echo "Part 2 (After TOAR)"
echo "-----------------------------------"
echo "It will create for *.surf.* images from i.atcorr"
echo "Atmospherically corrected image=surface reflectance"
#-----------------------------------------------------
# For i.atcorr scripting
#-----------------------------------------------------
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

# Location of parameter file
root=~/
# Basic script for i.atcorr for L 5 TM
#Geometrical conditions (L5TM)
geom=7
#Sensor height (satellite is -1000)
sens_height=-1000
#Atmospheric mode
atm_mode=6 #us standard 62 (for lack of more precise model)
#Aerosol model
aerosol_mode=1 #continental
#satellite band number (L5TM [25,26,27,28,29,30])
satbandno=25 #Band 1 of L5TM is first to undergo atmospheric correction

for file in $(g.mlist type=rast pattern=*.toar.1)
do
 #Here we suppose you have altitude (DEM) and Visibility (VIS) 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)
 #Altitude dummy value (in Km should be negative in this param file)
 #(overwritten by DEM raster input)
 alt=-1.200
 # L5 basename as stored in GRASS GIS and used by i.landsat.toar
 L5basename=$(echo $file | sed 's/\(.*\)\.\(.*\)\.\(.*\)/\1/')
 echo $L5basename

 #---------------------------
 # Please change as you need
 #---------------------------

 #datetime of satellite overpass (month, day, GMT decimal hour)
 echo "Input: GMT (i.e. 6.30)"
 read gmt
 #mdh="6 03 6.30"
 monthday=$(echo $L5basename | sed 's/\(.*\)_.......\(.*\)/\2/')
 day=$(echo $L5basename | sed 's/\(.*\)_.........\(.*\)/\2/')
 month=$(echo "($monthday - $day) / 100" | bc)
 mdh="$month $day $gmt"
 echo $mdh

 # Central Lat/Long
 north=$(g.region -p | grep north | sed 's/north:\ \(.*\)/\1/' | bc)
 south=$(g.region -p | grep south | sed 's/south:\ \(.*\)/\1/' | bc)
 east=$(g.region -p | grep east | sed 's/east:\ \(.*\)/\1/' | bc)
 west=$(g.region -p | grep west | sed 's/west:\ \(.*\)/\1/' | bc)
 Lat_nonproj=$(echo "(($north - $south)/2.0) + $south" | bc )
 Long_nonproj=$(echo "(($east - $west)/2.0) + $west" | bc )
 echo "$Long_nonproj $Lat_nonproj" > tempfile.txt
 Long=$(m.proj -o -d input=tempfile.txt | sed
's/\(.*\)|\(.*\)|\(.*\)/\1/')
 Lat=$(m.proj -o -d input=tempfile.txt | sed
's/\(.*\)|\(.*\)|\(.*\)/\2/')
 echo $Long_nonproj $Lat_nonproj
 echo $Long $Lat

 for bandno in 1 2 3 4 5 7
 do # Generate the parameterization file
 echo "$geom                            - geometrical
conditions=Landsat 5 TM" > $root/param_L5.txt
 echo "$mdh $Long $Lat   - month day hh.ddd longitude latitude
(\"hh.ddd\" is in decimal hours GMT)" >>  $root/param_L5.txt
 echo "$atm_mode                            - atmospheric
mode=tropical" >> $root/param_L5.txt
 echo "$aerosol_mode                            - aerosols
model=continental" >> $root/param_L5.txt
 echo "$vis                           - visibility [km] (aerosol
model concentration), not used as there is raster input" >>
$root/param_L5.txt
 echo "$alt                       - mean target elevation above sea
level [km] (here 600m asl), not used as there is raster input" >>
$root/param_L5.txt
 echo "$sens_height                        - sensor height (here,
sensor on board a satellite)" >> $root/param_L5.txt
 echo "$satbandno                           - 'i'th band of TM
Landsat 5" >> $root/param_L5.txt
 # Process band-wise atmospheric correction with 6s
 cat $root/param_L5.txt
 echo "i.atcorr -r input=$L5basename.toar.$bandno elevation=$dem
visibility=visibility parameters=$root/param_L5.txt
output=$L5basename.surf.$bandno range=0,1 rescale=0,1 --overwrite"
 i.atcorr -r input=$L5basename.toar.$bandno elevation=$dem
visibility=visibility parameters=$root/param_L5.txt
output=$L5basename.surf.$bandno range=0,1 rescale=0,1 --overwrite
 satbandno=$((satbandno+1))
 if [ $satbandno -gt 30 ]
  then satbandno=25
 fi
 done
done


-----Original Message-----
From: grass-user-bounces at lists.osgeo.org
[mailto:grass-user-bounces at lists.osgeo.org] On Behalf Of Elena Mezzini
Sent: Saturday, May 14, 2011 12:14 AM
To: grass-user at lists.osgeo.org
Subject: [GRASS-user] Re: i.atcorr returns all NULL values

I have the same problem, both with band 1 and 2.
I did i.landsat.toar and then i.atcorr but all the bands worked except
for band 1 and 2, these have values min=-nan and max=-nan.
There is someone who can help us?

Thank you!
Elena