Bugs in r.sun

[ojni0001]

I am new to this, so some things I write here may not be how it should be. My environment: Win7 x64 GRASS Version (7beta and RC1) and hence 32 bit binary (although I mentioned GRASS7, but it may be present in GRASS 6, and earlier versions)

I have 2 issues, probably needs 2 tickets but I am mentioning here both.

I used a virtual landscape with elevation = constant (i.e. flat)

Issue 1.

I found that if I am using the aspect raster, zero degree is East and 270 is South (as mentioned in the help). But if I am using a single value for aspect, zero or 360 degrees is North and 180 degree is South. I don't know if 90 degree is East or West. So, either help document has to be changed or the algorithm.

Issue 2.

When the slope is more than 60 degrees (probably from 45 degrees) facing North (northwards from east or west and North), the global radiation values are increasing

i.e. the radiation value for slope of 70 degrees is more than when the slope is 60 degrees

the radiation value for slope of 90 degrees is more than when the slope is 80 degrees

Here, is a small table for demonstration (aspect of 0 or 360 is North, as mentioned in issue 1 above) Jan->17 (day=17), Feb->16 (day=47), Mar->16 (day=75) April->15 (day=105), May->15 (day=135)

slope aspect Jan Feb Mar Apr May

10 345 1.61 2.45 3.59 5.24 6.52

20 345 0.93 1.73 2.93 4.70 6.18

30 345 0.36 1.04 2.21 4.03 5.67

40 345 0.13 0.46 1.49 3.26 5.01

50 345 0.12 0.15 0.85 2.42 4.22

60 345 0.11 0.16 0.85 2.38 4.17

70 345 0.10 0.59 1.63 3.38 5.12

80 345 0.54 1.29 2.47 4.30 5.94

90 345 1.20 2.05 3.28 5.11 6.62

10 360 1.59 2.42 3.57 5.23 6.52

20 360 0.87 1.66 2.87 4.67 6.17

30 360 0.27 0.91 2.11 3.99 5.66

40 360 0.13 0.26 1.30 3.21 5.00

50 360 0.12 0.14 0.47 2.36 4.21

60 360 0.11 0.13 0.56 2.48 4.35

70 360 0.10 0.35 1.51 3.49 5.32

80 360 0.38 1.12 2.42 4.41 6.15

90 360 1.08 1.97 3.28 5.22 6.81

For testing I have attached a zipped file (simulation for slope 0 to 90, step 10 degrees and for aspect 0 to 360, step 15 degrees) with the script and sample elevation (flat landscape) file. The table may look a bit different but the pattern will be similar.

Thank you.

Nirmal

[ojni0001]

test elevation, script and output (r.info) files

[wenzeslaus]

Replying to ojni0001:

For testing I have attached a zipped file (simulation for slope 0 to 90, step 10 degrees and for aspect 0 to 360, step 15 degrees) with the script and sample elevation (flat landscape) file. The table may look a bit different but the pattern will be similar.

Can you create a test which would be able to run in a fully automated way? Something like:

• source:grass/trunk/temporal/t.rast.accdetect/testsuite

General information about writing tests is here:

• ​http://grass.osgeo.org/grass71/manuals/libpython/gunittest_testing.html

[ojni0001]

I can try creating a test, but it will take some time. I haven't really used Python before and I created the shell script first time for testing above simulation. However, by looking at the page you provided, I may be able to create one.

---

Nirmal

Replying to wenzeslaus:

Replying to ojni0001:

For testing I have attached a zipped file (simulation for slope 0 to 90, step 10 degrees and for aspect 0 to 360, step 15 degrees) with the script and sample elevation (flat landscape) file. The table may look a bit different but the pattern will be similar.

Can you create a test which would be able to run in a fully automated way? Something like:

• source:grass/trunk/temporal/t.rast.accdetect/testsuite

General information about writing tests is here:

• ​http://grass.osgeo.org/grass71/manuals/libpython/gunittest_testing.html

[ojni0001]

Creating Python script is not possible for me at the moment. I have created a shell script running in Windows7 OS. So, I guess it also runs on Linux machine. The script does not assume that data is present. It will create a raster data with 1 cell and runs the test.

I cannot find a way to attach the code. Please see the script below:

function fcomp() {
    awk -v n1=$1 -v n2=$2 'BEGIN{ if (n1<n2) exit 0; exit 1}'
}
echo "running test for ticket no. 2606" > test_2606.txt
g.mapset PERMANENT
g.proj -c epsg=2449
g.region s=-83700 n=-83600 w=-23200 e=-23100 res=100
r.mapcalc --overwrite "const_elev = 15.0"
days=( 17 47 75 105 135 162 198 228 258 288 318 344)
month=( jan feb mar apr may jun jul aug sep oct nov dec)
aspect=( 0.0 15.0 30.0 45.0 60.0 75.0 90.0 105.0 120.0 135.0 150.0 165.0 180.0 195.0 210.0 225.0 240.0 255.0 270.0 285.0 300.0 315.0 330.0 345.0 360.0)
slope=( 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0)
for k in {0..11}
do 
        for i in {0..24}
        do
                for j in {0..9}
                do              
                        CMD="r.sun --overwrite --quiet elevation=const_elev aspect_value=${aspect[$i]} slope_value=${slope[$j]} linke_value=1.0 glob_rad=at1-g_${month[$k]}_${slope[$j]}_${aspect[$i]} day=${days[$k]} step=0.5" 
                        $CMD
                done
        done    
done

fail1=0
fail2=0
fail3=0
fail4=0
count1=0
count2=0
count3=0
for k in {0..11}
do 
        #for i in {0..24}
        #lets check at aspect 0 and 180 degrees, the radiation value should be equal if 270 is south
        for i in 0 12
        do
                ((i2=i+12))
                for j in {0..9}
                do
                        x1=`r.info -r at1-g_${month[$k]}_${slope[$j]}_${aspect[$i]} | grep "max" | cut -d =  -f 2 `
                        x2=`r.info -r at1-g_${month[$k]}_${slope[$j]}_${aspect[$i2]} | grep "max" | cut -d =  -f 2 `
                        diff_val=`awk "BEGIN {print $x1-$x2; exit}"`
                        echo "  month: ${month[$k]}     aspect: ${aspect[$i]}; ${aspect[$i2]}   slope=${slope[$j]}" >> test_2606.txt
                        echo "          x1=$x1  x2=$x2  difference: $diff_val" >> test_2606.txt
                        ((count1=count1+1))
                        if [ $diff_val = 0 ];then
                                echo "          test passed" >> test_2606.txt
                        else
                                echo "          test failed; bug remaining, east/west is not 0/180, south is not 270 degrees" >> test_2606.txt
                                ((fail1=fail1+1))
                        fi
                done
        done
        wait
done
wait
for k in {0..11}
do 
        #for i in {0..24}
        #lets check at aspect 90 and 270 degrees to check whether 90/270 is east/west or west/east
        for i in 6
        do
                ((i2=i+12))
                for j in {0..9}
                do
                        x1=`r.info -r at1-g_${month[$k]}_${slope[$j]}_${aspect[$i]} | grep "max" | cut -d =  -f 2 `
                        x2=`r.info -r at1-g_${month[$k]}_${slope[$j]}_${aspect[$i2]} | grep "max" | cut -d =  -f 2 `
                        diff_val=`awk "BEGIN {print $x1-$x2; exit}"`
                        ((count2=count2+1))
                        echo "  month: ${month[$k]}     aspect: ${aspect[$i]}; ${aspect[$i2]}   slope=${slope[$j]}" >> test_2606.txt
                        echo "          x1=$x1  x2=$x2  difference: $diff_val" >> test_2606.txt
                        if [ $diff_val = 0 ];then
                                echo "          values are equal at 90/270 bug remaining, east/west is 90/270 degrees" >> test_2606.txt
                                ((fail2=fail2+1))
                        fi
                done
        done
        wait
done
wait
for k in {0..11}
do 
        for i in {0..24}
        #for i in 6
        do
                #for j in {0..9}
                for j in 7
                do
                        ((j2=j+1))
                        ((j3=j+2))
                        x1=`r.info -r at1-g_${month[$k]}_${slope[$j]}_${aspect[$i]} | grep "max" | cut -d =  -f 2 `
                        x2=`r.info -r at1-g_${month[$k]}_${slope[$j2]}_${aspect[$i]} | grep "max" | cut -d =  -f 2 `
                        x3=`r.info -r at1-g_${month[$k]}_${slope[$j3]}_${aspect[$i]} | grep "max" | cut -d =  -f 2 `
                        diff_val1=`awk "BEGIN {print $x1-$x2; exit}"`
                        diff_val2=`awk "BEGIN {print $x2-$x3; exit}"`
                        ((count3=count3+1))
                        echo "  month: ${month[$k]}     aspect: ${aspect[$i]}   slope=${slope[$j]}      ${slope[$j2]}   ${slope[$j3]}" >> test_2606.txt
                        echo "          x1=$x1  x2=$x2  x3=$x3  differences: $diff_val1 $diff_val2" >> test_2606.txt
                        #if [ $diff_val1 > 0 ];then
                        if fcomp 0 $diff_val1; then
                                echo "          $diff_val1 > 0; OK" >> test_2606.txt
                        else
                                echo "          $diff_val1 < 0;bug remaining, ${slope[$j]} has less radiance than ${slope[$j2]}" >> test_2606.txt
                                ((fail3=fail3+1))
                        fi
                        #if [ $diff_val2 > 0 ];then
                        if fcomp 0 $diff_val2; then
                                echo "          $diff_val2 > 0; OK" >> test_2606.txt
                        else
                                echo "          $diff_val2 < 0;bug remaining, ${slope[$j2]} has less radiance than ${slope[$j3]}" >> test_2606.txt
                                ((fail4=fail4+1))
                        fi
                done
        done
        wait
done
wait
echo "e-w test fail:$fail1 of $count1;  south test fail:$fail2 of $count2" >> test_2606.txt
echo "e-w test fail:$fail1 of $count1;  south test fail:$fail2 of $count2"
echo "values for slope 70 and 80 fail:$fail3 of $count3;        values for slope 80 and 90 fail:$fail4 of $count3" >> test_2606.txt
echo "values for slope 70 and 80 fail:$fail3 of $count3;        values for slope 80 and 90 fail:$fail4 of $count3"

[neteler]

Ticket retargeted after 7.0.1 milestone closed

[neteler]

Ticket retargeted after milestone closed

[neteler]

Ticket retargeted after milestone closed

[neteler]

Ticket retargeted after 7.0.3 milestone closed

[martinl]

Still relevant?

[neteler]

see also #498