|
@@ -3,7 +3,7 @@
|
|
|
# Script to test i.topo.corr with a synthetic map
|
|
|
#
|
|
|
# Use North Carolina location to test:
|
|
|
-# grass71 ~/grassdata/nc_spm_08_grass7/user1
|
|
|
+# grass64 ~/grassdata/nc_spm_08/user1
|
|
|
|
|
|
if test "$GISBASE" = ""; then
|
|
|
echo "You must be in GRASS to run this program."
|
|
@@ -25,34 +25,36 @@ YEAR=`echo $DATETIME | cut -d' ' -f3 | awk '{printf "%d", $1}'`
|
|
|
TMPTIME=`echo $DATETIME | cut -d' ' -f4 | awk '{printf "%d", $1}'`
|
|
|
HOUR=`echo $TMPTIME | cut -d':' -f1 | awk '{printf "%d", $1}'`
|
|
|
MIN=`echo $TMPTIME | cut -d':' -f2 | awk '{printf "%d", $1}'`
|
|
|
+SEC=`echo $TMPTIME | cut -d':' -f3 | awk '{printf "%d", $1}'`
|
|
|
TIMEZ=`echo $DATETIME | cut -d' ' -f5 | awk '{printf "%d", $1/100}'`
|
|
|
unset TMPTIME
|
|
|
|
|
|
# create synthetic DEM (kind of roof)
|
|
|
-r.plane myplane0 dip=45 az=0 east=637500 north=221750 elev=1000 type=FCELL
|
|
|
-r.plane myplane90 dip=45 az=90 east=684800 north=221750 elev=1000 type=FCELL
|
|
|
-r.plane myplane180 dip=45 az=180 east=684800 north=260250 elev=1000 type=FCELL
|
|
|
-r.plane myplane270 dip=45 az=270 east=684800 north=221750 elev=1000 type=FCELL
|
|
|
+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.)"
|
|
|
|
|
|
# nviz
|
|
|
# nviz myplane_pyr
|
|
|
|
|
|
# get sun position
|
|
|
-eval `r.sunmask -s -g elev=myplane_pyr year=$YEAR month=8 day=$DAY hour=$HOUR minute=$MIN timezone=$TIMEZ`
|
|
|
+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 input=myplane_pyr output=myplane_pyr_shaded altitude=$sunangleabovehorizon azimuth=$sunazimuth
|
|
|
-d.mon stop=wx0 2> /dev/null
|
|
|
-d.mon wx0
|
|
|
-sleep 5 # this is rather annoying
|
|
|
-d.rast myplane_pyr_shaded
|
|
|
+# show raw map as shaded map
|
|
|
+#d.mon wx0
|
|
|
+#sleep 5 # this is rather annoying
|
|
|
+#d.rast myplane_pyr_shaded
|
|
|
+#echo "Original as shaded map" | d.text color=black
|
|
|
|
|
|
# pre-run: illumination map
|
|
|
-i.topo.corr -i input=myplane_pyr output=myplane_pyr_illumination \
|
|
|
+i.topo.corr -i output=myplane_pyr_illumination \
|
|
|
basemap=myplane_pyr zenith=$solarzenith azimuth=$sunazimuth
|
|
|
r.colors myplane_pyr_illumination color=gyr
|
|
|
|
|
@@ -67,47 +69,44 @@ echo "Original" | d.text color=black
|
|
|
|
|
|
# making the 'band' reflectance file from the shade map
|
|
|
r.mapcalc "myplane_pyr_band = double((myplane_pyr_shaded - 60.)/18.)"
|
|
|
+echo "Band map statistics: reflectance values:"
|
|
|
+r.univar -g myplane_pyr_band
|
|
|
r.colors myplane_pyr_band color=gyr
|
|
|
-d.mon stop=wx1 2> /dev/null
|
|
|
-d.mon wx1
|
|
|
-sleep 5 # this is rather annoying
|
|
|
-d.rast myplane_pyr_band
|
|
|
-d.legend myplane_pyr_band
|
|
|
-echo "Band reflectance" | d.text color=black
|
|
|
+#d.mon wx1
|
|
|
+#sleep 5 # this is rather annoying
|
|
|
+#d.rast myplane_pyr_band
|
|
|
+#d.legend myplane_pyr_band
|
|
|
+#echo "Band reflectance" | d.text color=black
|
|
|
|
|
|
## test it:
|
|
|
# percent
|
|
|
METHOD=percent
|
|
|
i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
|
|
|
-d.mon stop=wx2 2> /dev/null
|
|
|
-d.mon wx2
|
|
|
-sleep 5 # this is rather annoying
|
|
|
-d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
|
|
|
-echo "METHOD=percent" | d.text color=black
|
|
|
+#d.mon wx2
|
|
|
+#sleep 5 # this is rather annoying
|
|
|
+#d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
|
|
|
+#echo "METHOD=percent" | d.text color=black
|
|
|
|
|
|
# minnaert
|
|
|
METHOD=minnaert
|
|
|
i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
|
|
|
-d.mon stop=wx3 2> /dev/null
|
|
|
-d.mon wx3
|
|
|
-sleep 5 # this is rather annoying
|
|
|
-d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
|
|
|
-echo "METHOD=minnaert" | d.text color=black
|
|
|
+#d.mon wx3
|
|
|
+#sleep 5 # this is rather annoying
|
|
|
+#d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
|
|
|
+#echo "METHOD=minnaert" | d.text color=black
|
|
|
|
|
|
# c-factor
|
|
|
METHOD=c-factor
|
|
|
i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
|
|
|
-d.mon stop=wx4 2> /dev/null
|
|
|
-d.mon wx4
|
|
|
-sleep 5 # this is rather annoying
|
|
|
-d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
|
|
|
-echo "METHOD=c-factor" | d.text color=black
|
|
|
+#d.mon wx4
|
|
|
+#sleep 5 # this is rather annoying
|
|
|
+#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=wx5 2> /dev/null
|
|
|
-d.mon wx5
|
|
|
-sleep 5 # this is rather annoying
|
|
|
-d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
|
|
|
-echo "METHOD=cosine" | d.text color=black
|
|
|
+#d.mon wx5
|
|
|
+#sleep 5 # this is rather annoying
|
|
|
+#d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
|
|
|
+#echo "METHOD=cosine" | d.text color=black
|