Resending as attachment...
-----Original Message-----
From: grass-user-bounces@lists.osgeo.org
[mailto:grass-user-bounces@lists.osgeo.org] On Behalf Of Chemin, Yann
(IWMI)
Sent: Monday, May 23, 2011 8:37 AM
To: grass-user@lists.osgeo.org
Cc: Elena Mezzini
Subject: RE: [GRASS-user] Re: i.atcorr returns all NULL values
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@lists.osgeo.org
[mailto:grass-user-bounces@lists.osgeo.org] On Behalf Of Elena Mezzini
Sent: Saturday, May 14, 2011 12:14 AM
To: grass-user@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
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user
(attachments)
L5TM_import_TOAR_ATCORR.sh (4.99 KB)