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-boun...@lists.osgeo.org [mailto:grass-user-boun...@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