[GRASS-SVN] r47482 - grass/branches/releasebranch_6_4/imagery/i.topo.corr

[svn_grass at osgeo.org]

Author: neteler
Date: 2011-08-07 09:48:13 -0700 (Sun, 07 Aug 2011)
New Revision: 47482

Modified:
   grass/branches/releasebranch_6_4/imagery/i.topo.corr/test_i.topo.corr_synthetic_DEM_NC.sh
Log:
added cosine test; cosmetics

Modified: grass/branches/releasebranch_6_4/imagery/i.topo.corr/test_i.topo.corr_synthetic_DEM_NC.sh
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.topo.corr/test_i.topo.corr_synthetic_DEM_NC.sh	2011-08-07 14:45:13 UTC (rev 47481)
+++ grass/branches/releasebranch_6_4/imagery/i.topo.corr/test_i.topo.corr_synthetic_DEM_NC.sh	2011-08-07 16:48:13 UTC (rev 47482)
@@ -12,7 +12,7 @@
 
 export GRASS_OVERWRITE=1
 
-# we use the NC location
+# we use the NC location, time zone UTM-5
 g.region n=308500 s=215000 w=630000 e=704800 nsres=250 ewres=250 -pa
 
 # note: no daylight saving time in summer 1954!
@@ -29,12 +29,12 @@
 TIMEZ=`echo $DATETIME | cut -d' ' -f5 | awk '{printf "%d", $1/100}'`
 unset TMPTIME
 
-# create synthetic DEM (kind of pyramid)
+# create synthetic DEM (kind of roof)
 r.plane --o myplane0 dip=45 az=0 east=637500 north=221750 elev=1000 type=float
 r.plane --o myplane90 dip=45 az=90 east=684800 north=221750 elev=1000 type=float
 r.plane --o myplane180 dip=45 az=180 east=684800 north=260250 elev=1000 type=float
 r.plane --o myplane270 dip=45 az=270 east=684800 north=221750 elev=1000 type=float
-r.mapcalc "myplane_pyr = double(min(myplane90,myplane270,myplane0,myplane180)/10 + 8600)"
+r.mapcalc "myplane_pyr = double(min(myplane90,myplane270,myplane0,myplane180)/10. + 8600.)"
 
 # nviz
 # nviz myplane_pyr
@@ -43,27 +43,35 @@
 eval `r.sunmask -s -g output=dummy elev=myplane_pyr year=$YEAR month=8 day=$DAY hour=$HOUR minute=$MIN second=$SEC timezone=$TIMEZ`
 
 solarzenith=`echo $sunangleabovehorizon | awk '{printf "%f", 90. - $1}'`
+echo "Sun position ($DATETIME): solarzenith: $solarzenith, sunazimuth: $sunazimuth"
 
 # shade relief
 r.shaded.relief map=myplane_pyr shadedmap=myplane_pyr_shaded altitude=$sunangleabovehorizon azimuth=$sunazimuth
-d.mon stop=x0
+d.mon stop=x0 > /dev/null
 d.mon x0
 d.rast myplane_pyr_shaded
 
 # pre-run: illumination map
-i.topo.corr -i input=myplane_pyr output=myplane_pyr_illumination basemap=myplane_pyr zenith=$solarzenith azimuth=$sunazimuth method=c-factor
+i.topo.corr -i input=myplane_pyr output=myplane_pyr_illumination \
+	    basemap=myplane_pyr zenith=$solarzenith azimuth=$sunazimuth 
+r.colors myplane_pyr_illumination color=gyr
 
 # show original
 d.erase -f
-r.colors myplane_pyr col=grey
+r.colors myplane_pyr color=grey
 d.rast.leg myplane_pyr
+echo "Original" | d.text color=black
 
-#making the 'band' reflectance file
+# nviz with illumination draped over
+# nviz myplane_pyr color=myplane_pyr_illumination
+
+# making the 'band' reflectance file from the shade map
 r.mapcalc "myplane_pyr_band = double((myplane_pyr_shaded - 60.)/18.)"
+r.colors myplane_pyr_band color=gyr
 d.mon stop=x0 > /dev/null
 d.mon x0
 d.rast.leg myplane_pyr_band
-echo "Original" | d.text color=black
+echo "Band reflectance" | d.text color=black
 
 ## test it:
 # percent
@@ -90,4 +98,12 @@
 d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
 echo "METHOD=c-factor" | d.text color=black
 
+# cosine
+METHOD=cosine
+i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+d.mon stop=x4 > /dev/null
+d.mon x4
+d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
+echo "METHOD=cosine" | d.text color=black
+
 exit 0