[GRASS-dev] [GRASS GIS] #2606: Bugs in r.sun

GDAL · March 2, 2015, 2:56am

#2606: Bugs in r.sun
----------------------+-----------------------------------------------------
Reporter: ojni0001 | Owner: grass-dev@…
Type: defect | Status: new
Priority: normal | Milestone:
Component: Default | Version: unspecified
Keywords: | Platform: MSWindows 7
Cpu: x86-64 |
----------------------+-----------------------------------------------------
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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2606>
GRASS GIS <http://grass.osgeo.org>

GDAL · March 2, 2015, 10:24am

#2606: Bugs in r.sun
----------------------+-----------------------------------------------------
Reporter: ojni0001 | Owner: grass-dev@…
Type: defect | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: unspecified
Keywords: r.sun | Platform: MSWindows 7
Cpu: x86-64 |
----------------------+-----------------------------------------------------
Changes (by martinl):

  * keywords: => r.sun
  * component: Default => Raster
  * milestone: => 7.0.1

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2606#comment:1>
GRASS GIS <http://grass.osgeo.org>

GDAL · March 18, 2015, 5:56pm

#2606: Bugs in r.sun
----------------------+-----------------------------------------------------
Reporter: ojni0001 | Owner: grass-dev@…
Type: defect | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: unspecified
Keywords: r.sun | Platform: MSWindows 7
Cpu: x86-64 |
----------------------+-----------------------------------------------------

Comment(by wenzeslaus):

Replying to [ticket:2606 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

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/2606#comment:2>
GRASS GIS <http://grass.osgeo.org>

GDAL · March 25, 2015, 7:11am

#2606: Bugs in r.sun
----------------------+-----------------------------------------------------
Reporter: ojni0001 | Owner: grass-dev@…
Type: defect | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: unspecified
Keywords: r.sun | Platform: MSWindows 7
Cpu: x86-64 |
----------------------+-----------------------------------------------------

Comment(by 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 [comment:2 wenzeslaus]:
> Replying to [ticket:2606 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

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2606#comment:3>
GRASS GIS <http://grass.osgeo.org>

GDAL · April 7, 2015, 4:30am

#2606: Bugs in r.sun
----------------------+-----------------------------------------------------
Reporter: ojni0001 | Owner: grass-dev@…
Type: defect | Status: new
Priority: normal | Milestone: 7.0.1
Component: Raster | Version: unspecified
Keywords: r.sun | Platform: MSWindows 7
Cpu: x86-64 |
----------------------+-----------------------------------------------------

Comment(by 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:

{{{
#!sh
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"

}}}

--
Ticket URL: <https://trac.osgeo.org/grass/ticket/2606#comment:4>
GRASS GIS <http://grass.osgeo.org>