i.topo.corr manual: figure added; test script updated

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@59558 15284696-431f-4ddb-bdfa-cd5b030d7da7

imagery/i.topo.corr/Makefile

@@ -8,3 +8,7 @@ DEPENDENCIES = $(GISDEP) $(RASTERDEP)
 include $(MODULE_TOPDIR)/include/Make/Module.make
 
 default: cmd
+
+fig2pdf:
+	fig2dev -L pdf -p 1 -m 0.99 terrain_correction.fig | sed 's+Times-Bold+Helvetica-Bold+g' | sed 's+Times-Roman+Helvetica+g' | sed 's+Times+Helvetica+g' | sed 's+Helvetica-BoldOblique+Helvetica-BoldItal+g' | sed 's+Helvetica-Oblique+Helvetica-ReguItal+g' | sed 's+NimbusSanL-Bold+Helvetica-Bold+g' | sed 's+Helvetica-Bold-iso+Helvetica-Bold+g' | sed 's+Helvetica-iso+Helvetica+g' | sed 's+Helvetica+Helvetica+g' > terrain_correction.pdf
+

imagery/i.topo.corr/i.topo.corr.html

@@ -3,9 +3,17 @@
 <em>i.topo.corr</em> is used to topographically correct reflectance
 from imagery files, e.g. obtained with <em>i.landsat.toar</em>, using a
 sun illumination terrain model. This illumination model represents the
-cosine of the incident angle, i.e. the  angle between the normal to the
-ground and the sun rays. It can be obtained with <em>r.sun</em>
-(parameter incidout), and then calculating its cosine with float precision.
+cosine of the incident angle <em>i</em>, i.e. the  angle between the normal to the
+ground and the sun rays.
+<p>
+Note: If needed, the sun position can be calculated for a given date with
+<em>r.sunmask</em>.
+
+
+<center>
+<img src="terrain_correction.png" alt="Figure showing terrain and solar angles"><br>
+Figure showing terrain and solar angles
+</center>
 
 <p>
 Using the <b>-i</b> flag and given an elevation basemap (metric),
@@ -14,8 +22,11 @@ Using the <b>-i</b> flag and given an elevation basemap (metric),
   <li> cos_i = cos(s) * cos(z) + sin(s) * sin(z) * cos(a - o) </li>
 </ul>
 where,
-<em>s</em> is the terrain slope angle, <em>z</em> is the solar zenith angle,
-<em>a</em> the solar azimuth angle, <em>o</em> the terrain aspect angle.
+<em>i</em> is the incident angle to be calculated,
+<em>s</em> is the terrain slope angle,
+<em>z</em> is the solar zenith angle,
+<em>a</em> the solar azimuth angle,
+<em>o</em> the terrain aspect angle.
 
 <p>
 For each band file, the corrected reflectance (ref_c) is calculate from
@@ -78,7 +89,10 @@ method:
 <p>
 
 <div class="code"><pre>
+# first pass: create illumination model
 i.topo.corr -i base=SRTM zenith=33.3631 azimuth=59.8897 out=SRTM.illumination
+
+# second pass: apply illumination model
 i.topo.corr base=SRTM.illumination input=toar.5,toar.4,toar.3 out=tcor \ 
   zenith=33.3631 method=c-factor
 </pre></div>
@@ -110,13 +124,15 @@ i.topo.corr base=SRTM.illumination input=toar.5,toar.4,toar.3 out=tcor \
 <a href="i.landsat.toar.html">i.landsat.toar</a>,
 <a href="r.mapcalc.html">r.mapcalc</a>,
 <a href="r.sun.html">r.sun</a>
+<a href="r.sunmask.html">r.sunmask</a>
 </em>
 
-
 <h2>AUTHOR</h2>
 
 E. Jorge Tizado  (ej.tizado unileon es)<br>
 Dept. Biodiversity and Environmental Management, University of Le&oacute;n, Spain
+<p>
+Figure derived from Neteler &amp; Mitasova, 2008.
 
 <p>
 <i>Last changed: $Date$</i>

imagery/i.topo.corr/terrain_correction.fig

@@ -0,0 +1,180 @@
+#FIG 3.2  Produced by xfig version 3.2.5c
+Portrait
+Center
+Metric
+A4      
+100.00
+Single
+-2
+1200 2
+5 1 0 1 0 7 0 0 -1 0.000 0 0 0 0 3101.053 2566.737 2990 2472 3120 2422 3240 2522
+5 1 0 1 0 7 0 0 -1 4.000 0 0 0 0 3407.500 2489.500 3410 2372 3510 2432 3510 2547
+5 1 2 1 0 7 0 0 -1 3.000 0 1 0 1 3012.345 2608.550 2415 2522 3025 3212 3615 2642
+	1 0 1.00 60.00 120.00
+5 1 0 1 0 7 0 0 -1 0.000 0 1 0 1 3037.739 2624.146 2692 2692 3105 2970 3390 2617
+	1 0 1.00 60.00 120.00
+5 1 0 1 0 7 0 0 -1 0.000 0 1 1 0 2402.636 2238.516 3105 1077 1985 947 1305 1440
+	1 0 1.00 60.00 120.00
+5 1 0 1 0 7 0 0 -1 0.000 0 0 1 0 2483.438 2359.688 990 2295 1080 1845 1305 1440
+	1 0 1.00 60.00 120.00
+5 1 0 1 0 7 0 0 -1 0.000 0 0 0 1 2581.702 2585.417 1941 1123 2520 990 3105 1077
+	1 0 1.00 60.00 120.00
+5 1 0 1 0 7 0 0 -1 0.000 0 0 0 1 2323.603 2036.782 1467 1536 1647 1311 1923 1129
+	1 0 1.00 60.00 120.00
+6 44 524 510 980
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 361 818 384 825 402 839 408 855 412 867 412 879
+	 418 894 427 916 422 928 415 911 403 898 393 890
+	 376 882 356 875 339 863 332 842
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 337 842 355 855 367 876 368 892 367 906 363 916
+	 363 932 363 956 356 965 354 948 349 930 341 921
+	 328 907 311 894 300 877 301 855
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 304 855 317 874 321 898 316 913 308 925 302 933
+	 296 949 288 971 277 977 282 960 283 941 279 931
+	 271 914 261 897 256 876 265 855
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 266 855 271 879 266 901 256 914 245 923 235 930
+	 225 942 208 960 196 962 208 948 215 931 216 919
+	 215 901 211 880 213 862 229 845
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 223 662 205 649 193 628 193 612 194 598 196 587
+	 197 572 197 547 205 538 206 556 213 574 220 583
+	 232 596 250 609 261 626 259 648
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 198 826 187 845 168 860 151 864 137 865 125 865
+	 109 866 83 870 72 865 91 861 108 852 117 842
+	 128 829 139 810 154 797 177 795
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 179 795 161 809 139 816 121 814 108 810 97 806
+	 80 803 56 798 47 789 66 792 85 789 96 783
+	 112 773 129 761 148 753 171 758
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 173 758 152 766 128 765 111 758 101 748 92 742
+	 77 734 56 721 51 710 68 719 88 722 99 721
+	 118 717 138 710 160 709 178 721
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 178 722 154 723 133 714 121 702 112 690 107 681
+	 96 669 80 649 80 638 93 651 110 662 122 663
+	 141 665 163 665 183 671 196 688
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 195 688 173 681 154 667 148 651 144 638 143 627
+	 137 612 129 590 133 578 141 595 152 609 164 614
+	 180 624 200 631 217 643 224 664
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 258 648 246 630 242 607 248 592 254 579 260 570
+	 266 554 275 533 287 527 280 543 280 562 285 573
+	 291 590 302 608 307 627 297 647
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 297 646 292 625 297 603 307 590 318 580 328 574
+	 339 562 355 544 367 542 356 556 349 572 347 585
+	 349 603 353 622 350 643 334 659
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 331 659 334 637 347 617 360 608 374 603 386 600
+	 400 591 422 579 435 581 419 592 406 605 401 616
+	 396 634 392 653 382 671 361 681
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 359 681 370 662 389 648 405 644 419 643 431 644
+	 448 641 474 638 484 643 465 647 448 655 439 665
+	 429 679 418 696 402 710 379 713
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 376 712 394 697 416 690 433 692 447 696 458 700
+	 474 703 500 708 508 717 489 715 470 717 458 723
+	 443 733 425 746 406 754 384 748
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 381 748 403 740 426 740 441 748 453 756 462 764
+	 476 771 499 784 503 795 485 787 466 783 454 785
+	 436 789 415 796 394 797 374 785
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 376 784 401 783 423 792 435 804 443 815 448 825
+	 458 837 474 855 474 868 461 855 445 845 432 842
+	 413 840 392 839 373 834 358 816
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+3 2 0 1 0 7 0 0 -1 4.000 0 0 0 16
+	 229 849 226 871 213 890 199 899 185 905 174 908
+	 160 915 137 927 124 926 141 915 154 902 158 891
+	 164 874 168 854 177 835 198 827
+	 0.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000
+	 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 -1.000 0.000
+-6
+6 25 2621 1425 3568
+4 0 0 50 0 16 8 0.0000 4 135 1395 30 3218 o: terrain aspect angle\001
+4 0 0 50 0 16 8 0.0000 4 135 1020 30 2888 i:  incident angle\001
+4 0 0 50 0 16 8 0.0000 4 135 990 30 2726 h: horizon angle\001
+4 0 0 50 0 16 8 0.0000 4 135 1260 25 3383 z: solar zenith angle\001
+4 0 0 50 0 16 8 0.0000 4 135 1380 25 3538 a: solar azimuth angle\001
+4 0 0 50 0 16 8 0.0000 4 135 1305 30 3053 s: terrain slope angle\001
+-6
+1 3 0 1 0 7 50 0 -1 0.000 1 0.0000 3150 2500 11 11 3150 2500 3161 2500
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 3150 2520 20 20 3150 2520 3170 2520
+2 1 1 1 0 7 0 0 -1 4.000 0 0 -1 0 0 2
+	 3105 2619 1372 2902
+2 3 0 1 0 7 0 0 -1 0.000 0 0 -1 0 0 4
+	 3960 1940 3960 2480 3105 2615 3960 1940
+2 1 1 1 0 7 0 0 -1 4.000 0 0 -1 0 0 3
+	 2520 2480 3375 2345 3960 2480
+2 1 1 1 0 7 0 0 -1 4.000 0 0 -1 0 0 2
+	 3375 2345 3375 1850
+2 1 0 2 0 7 0 0 -1 0.000 0 0 -1 0 0 2
+	 3105 2617 520 922
+2 1 0 1 0 7 50 0 -1 0.000 0 0 7 1 0 2
+	1 1 1.00 60.00 120.00
+	 3105 2617 4140 2617
+2 1 0 1 0 7 50 0 -1 0.000 0 0 7 1 0 2
+	1 1 1.00 60.00 120.00
+	 3105 2617 3105 627
+2 1 0 1 0 7 50 0 -1 0.000 0 0 7 1 0 2
+	1 1 1.00 60.00 120.00
+	 3105 2617 1625 747
+2 3 0 1 0 0 51 0 6 0.000 0 0 -1 0 0 5
+	 3113 2610 2527 2476 3382 1846 3968 1936 3113 2610
+2 1 2 1 0 7 0 0 -1 3.000 0 0 -1 0 0 3
+	 260 1030 260 2177 3105 2612
+4 0 0 50 0 16 8 0.0000 4 75 75 3420 2515 s\001
+4 0 0 50 0 16 8 0.0000 4 105 315 4205 2677 north\001
+4 0 0 50 0 16 8 0.0000 4 105 420 1395 585 normal\001
+4 0 0 50 0 16 8 0.0000 4 90 390 1395 720 vector\001
+4 0 0 50 0 16 8 0.0000 4 105 465 1395 450 surface\001
+4 0 0 50 0 16 8 0.1571 4 105 960 1438 2876 azimuth normal\001
+4 0 0 50 0 16 8 0.0000 4 75 75 2495 1096 s\001
+4 0 0 50 0 16 8 0.0000 4 105 30 1730 1402 i\001
+4 0 0 50 0 16 8 0.0000 4 105 75 1106 1944 h\001
+4 0 0 50 0 16 8 6.1436 4 105 840 315 2295 solar azimuth\001
+4 0 0 50 0 16 8 0.0000 4 75 75 3000 2940 o\001
+4 0 0 50 0 16 8 0.0000 4 75 75 2205 855 z\001
+4 0 0 50 0 16 8 0.0000 4 75 75 2970 3183 a\001
+4 0 0 50 0 16 8 6.1087 4 135 330 3405 1970 slope\001

imagery/i.topo.corr/test_i.topo.corr_synthetic_DEM_NC.sh

@@ -3,7 +3,7 @@
 # Script to test i.topo.corr with a synthetic map
 #
 # Use North Carolina location to test:
-#   grass64 ~/grassdata/nc_spm_08/user1
+#   grass71 ~/grassdata/nc_spm_08_grass7/user1
 
 if test "$GISBASE" = ""; then
  echo "You must be in GRASS to run this program."
@@ -25,22 +25,21 @@ 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 --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.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.mapcalc "myplane_pyr = double(min(myplane90,myplane270,myplane0,myplane180)/10. + 8600.)"
 
 # nviz
 # nviz myplane_pyr
 
 # get sun position
-eval `r.sunmask -s -g output=dummy elev=myplane_pyr year=$YEAR month=8 day=$DAY hour=$HOUR minute=$MIN second=$SEC timezone=$TIMEZ`
+eval `r.sunmask -s -g elev=myplane_pyr year=$YEAR month=8 day=$DAY hour=$HOUR minute=$MIN timezone=$TIMEZ`
 
 solarzenith=`echo $sunangleabovehorizon | awk '{printf "%f", 90. - $1}'`
 echo "Sun position ($DATETIME): solarzenith: $solarzenith, sunazimuth: $sunazimuth"
@@ -68,7 +67,7 @@ echo "Original" | d.text color=black
 # 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.mapcalc "myplane_pyr_band = double((myplane_pyr_shaded - 60.)/180.)"
 echo "Band map statistics: reflectance values:"
 r.univar -g myplane_pyr_band
 r.colors myplane_pyr_band color=gyr
@@ -78,10 +77,13 @@ r.colors myplane_pyr_band color=gyr
 #d.legend myplane_pyr_band
 #echo "Band reflectance" | d.text color=black
 
-## test it:
+## step 2
+## perform terrain flattening of image:
 # percent
 METHOD=percent
-i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} \
+  basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+r.colors myplane_pyr_topocorr_${METHOD}.myplane_pyr_band color=grey
 #d.mon wx2
 #sleep 5 # this is rather annoying
 #d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
@@ -89,7 +91,9 @@ i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap
 
 # minnaert
 METHOD=minnaert
-i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} \
+  basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+r.colors myplane_pyr_topocorr_${METHOD}.myplane_pyr_band color=grey
 #d.mon wx3
 #sleep 5 # this is rather annoying
 #d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
@@ -97,7 +101,9 @@ i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap
 
 # 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
+i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} \
+  basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+r.colors myplane_pyr_topocorr_${METHOD}.myplane_pyr_band color=grey
 #d.mon wx4
 #sleep 5 # this is rather annoying
 #d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band
@@ -105,7 +111,9 @@ i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap
 
 # cosine
 METHOD=cosine
-i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+i.topo.corr input=myplane_pyr_band output=myplane_pyr_topocorr_${METHOD} \
+  basemap=myplane_pyr_illumination zenith=$solarzenith method=$METHOD
+r.colors myplane_pyr_topocorr_${METHOD}.myplane_pyr_band color=grey
 #d.mon wx5
 #sleep 5 # this is rather annoying
 #d.rast.leg myplane_pyr_topocorr_${METHOD}.myplane_pyr_band