r.sun2 commissioning trials

[hamish]

Hi,

[6.5svn]

I was trying to figure out why r.sun was not using its short module description label. It seemed that somewhere it was reset to null, so the description was used regardless. In 6.3.0 it works. module_info.keywords also seemed to be reset to null.

r.sun --tcltk | head -n 3
begin_dialog {r.sun} {
 label {}
 desc {Computes direct (beam), diffuse and ...

but in r.sun/main.c:

   module->keywords = _("raster");
    module->label = _("Solar irradiance and irradiation model.");
    module->description =
	_("Computes direct (beam), diffuse and ...

same for --help and --interface-description. I noticed that (eg) 'r.info --help' reports the keyword.

When I loaded up which r.sun in Kdbg to follow the action live, I noticed the source listing had the ->keywords and ->label lines missing. wtf?!

aaaaaaahhhhhhh.... now I notice in the Kdbg title bar the full path to the source: ..../raster/r.sun2/main.c not raster/r.sun/main.c!

Apparently r.sun2 does not contain the latest changes to r.sun even though it has replaced it in the Makefile. As keyword support was added to r.sun/main.c on 08/19/06, and that's a while ago now, well, much unhappiness. not sure how far back r.sun2 goes.

attack of the out-of-sync clones & another "svn copy" related tragedy. (granted it was added to svn some months ago, when we were all young about these things)

Hamish

[neteler]

Replying to hamish:

Apparently r.sun2 does not contain the latest changes to r.sun even though it has replaced it in the Makefile. As keyword support was added to r.sun/main.c on 08/19/06, and that's a while ago now, well, much unhappiness. not sure how far back r.sun2 goes.

In r35938, r35939, r35940 (7, 6.4.svn, 6.4.0.svn) I have merged pre-fork updates from old r.sun into r.sun2.

attack of the out-of-sync clones & another "svn copy" related tragedy. (granted it was added to svn some months ago, when we were all young about these things)

No, the first tragedy. The "svn copy" ideal-world-solution came up only after that. Yes, sh*t happens. Note that r.sun2 was forked years ago and merging in changes of r.sun/GRASS took already hours. Oversights are possible.

TODO: understand usage of

if (latin == NULL && lt == NULL && (G_projection() != PROJECTION_LL)) {

which yet differs.

Markus

[neteler]

One more problem:

   r.sun -s elevin="pat_dtm_5m" aspin="pat_dtm_5m.as" aspect=270 slopei\
   n="pat_dtm_5m.sl" slope=0.0 linkein="linke_turbidity04" lin=3.0 alb=\
   0.2 horizon="pat5m_horangle" horizonstep=30 beam_rad="pat_beam_rad10\
   4" insol_time="pat_insol_time104" diff_rad="pat_diff_rad104" refl_ra\
   d="pat_refl_rad104" glob_rad="pat_glob_rad104" day=104 step=0.5 dist\
   =1.0 numpartitions=64

If aspin= and/or slopein= are used, aspect=270 and/or slope=0.0 should be suppressed in the history. How to do that? Currently the history is confusing.

Markus

[hamish]

Replying to neteler:

In r35938, r35939, r35940 (7, 6.4.svn, 6.4.0.svn) I have merged pre-fork updates from old r.sun into r.sun2.

cheers.

Replying to hamish:

attack of the out-of-sync clones & another "svn copy" related tragedy. (granted it was added to svn some months ago, when we were all young about these things)

MN:

No, the first tragedy. The "svn copy" ideal-world-solution came up only after that. Yes, sh*t happens. Note that r.sun2 was forked years ago and merging in changes of r.sun/GRASS took already hours. Oversights are possible.

whoa, I wasn't pointing that at you personally at all. sorry if it came across that way, it was meant as a general comment. I have spent too many hours merging stuff by hand to be anything but humble about it, and am certainly no svn expert at all. FWIW as the past example I was specifically thinking about the pain to merge years of here-and-there bug fixes between the four i.points clones some time ago.

TODO: understand usage of

if (latin == NULL && lt == NULL && (G_projection() != PROJECTION_LL)) {

which yet differs.

I will try and have a look. (just ran indent on r.sun2 to make it easier)

also I wonder about r19184. r.sun2 has the <= part, so I lightly guess the 0.5 part of the patch has been considered & dropped??

Hamish

[hamish]

fyi, my setup for testing r.sun will now be:

#spearfish
g.region -d
r.surf.volcano out=gauss method=gaussian
DAY=180
r.sun '-s'   # shading calc ON

it is much easier with a standardized hill which is simple enough that you can intuitively understand the result.

Hamish

[hamish]

.. make that using day 355 for testing. more noticeable shadows to follow.

• the incidout output map needs to be modified so fully-shaded reports NULL not 0 degrees incidence angle.

• for mode 2 I notice that I have to drop the time step to about 3 minutes to get a smooth output map (step=0.05)

H

[hamish]

I have been doing some experiments with r.sun + a simple gaussain mound as the elevation:

g.region -d   # spearfish dataset region
r.surf.volcano gauss method=gaussian   # module from addons

DAY=355   # longest shadows

I have looked at the effect of using r.horizon and slope/aspect maps. setup:

r.horizon elev=gauss horizonstep=30 dist=0.7 horizon=horangle
r.slope.aspect elevation=gauss aspect=guass.aspect slope=guass.slope

Q1) how small must the horizonstep be to recreate the beam_rad output map effectively as good as e.g. seen in the r.sun wiki page example?

​http://grass.osgeo.org/wiki/r.sun

By using a 30deg angle step are we trading time/maps for accuracy? For some reason the lower east wing of the butterfly pattern is stunted using that. Running tests with horizonstep=1 degree now.

Q2) am I using aspin= correctly? should that be a map created by r.slope.aspect? Using it the processing speeds up by 3x, but the results seem a bit weird -- to the south of the mound is a depression in the light which while possible for reflected light seems that it would be "upstream" & incorrect for beam_rad.

defining a slopein= map (from r.slope.aspect result) doesn't change the processing time appreciably. I haven't tested for its effect on the output map yet.

? Hamish

[hamish]

Replying to hamish:

Q2) am I using aspin= correctly? should that be a map created by r.slope.aspect? Using it the processing speeds up by 3x, but the results seem a bit weird -- to the south of the mound is a depression in the light which while possible for reflected light seems that it would be "upstream" & incorrect for beam_rad.

defining a slopein= map (from r.slope.aspect result) doesn't change the processing time appreciably. I haven't tested for its effect on the output map yet.

correction: using aspin=r.slope.aspect doesn't speed things up, but slope+aspect does. Using slopein= causes differences (missing south west butterfly wing), but slope+aspect makes it really bad.

test script:

time r.sun -s elevin=gauss day=$DAY beam_rad=rad_test.355.beam.try4
real    3m1.750s
user    3m1.327s
sys     0m0.400s
#

time r.sun -s elevin=gauss day=$DAY beam_rad=rad_test.355.beamA.try4  \
  aspin=guass.aspect
real    3m4.432s
user    3m4.092s
sys     0m0.364s
#

time r.sun -s elevin=gauss day=$DAY beam_rad=rad_test.355.beamS.try4 \
  slopein=guass.slope
real    3m7.865s
user    3m6.444s
sys     0m0.612s
#

time r.sun -s elevin=gauss day=$DAY beam_rad=rad_test.355.beamAS.try4 \
  aspin=guass.aspect slopein=guass.slope
real    1m32.279s
user    1m31.686s
sys     0m0.240s
#


r.contour in=gauss out=gauss_200m_contours step=200

plot_stuff()
{
MAP="rad_test.355.beam${MAPEXT}.try4"
d.rast "$MAP"
d.vect gauss_200m_contours color=180:180:180
echo "$1" | d.text color=black
eval `r.univar -g $MAP`
echo "sum: $sum" | d.text color=black at=1,5
d.legend "$MAP"
}

###
d.mon x1
d.resize w=640 h=480
d.split.frame 4
#
d.frame uno
MAPEXT=""
plot_stuff "Std."
#
d.frame dos
MAPEXT="A"
plot_stuff "w/Aspect map"
#
d.frame tres
MAPEXT="S"
plot_stuff "w/Slope map"
#
d.frame cuatro
MAPEXT="AS"
plot_stuff "w/Aspect and Slope maps"
###
d.frame full_screen
#d.out.file out=rsun_aspslp format=png

results in the attached image (rsun_aspslp.png).

Hamish

[hamish]

documenting r.sun aspin= slopein= errors

[hamish]

Q1) how small must the horizonstep be to recreate the beam_rad output map effectively as good as e.g. seen in the r.sun wiki page example?

trials showing differences between no horizon steps; trying to find the ideal "sweet spot" compromise between time and quality.

note that none are perfect. Standard r.sun with no horizon inputs is slightly clipped in lower east wing; the 360x 1 degree horizon map has slightly uneven width between inner east & west bands. note2 color tables are not constant between maps.

see attached screenshot (rsun_horizons.png)

test script:

for DEG in 1 15 30 ; do
   # can take a while to create all maps
   r.horizon elev=gauss horizonstep=$DEG dist=0.7 horizon=horangle$DEG

   time r.sun -s elevin=gauss day=$DAY horizon=horangle$DEG \
     horizonstep=$DEG beam_rad=rad_test.355.beam.Hz${DEG}deg.try5
done


# r.sun with no horizon seeds
real    3m7.571s
user    3m7.444s
sys     0m0.060s

# r.sun with 360x 1 deg seeds
#  sits quite a while at 0% due to loading sheer number of horizon maps!
real    1m20.385s
user    1m18.853s
sys     0m1.272s

# r.sun with 24x 15 deg seeds
real    0m32.185s
user    0m31.774s
sys     0m0.248s

# r.sun with 12x 30 deg seeds
real    0m29.518s
user    0m29.382s
sys     0m0.120s
#

plot_stuff()
{
r.colors $MAP -e color=bcyr
d.rast "$MAP"
d.vect gauss_200m_contours color=180:180:180
echo "$1" | d.text color=black
eval `r.univar -g $MAP`
echo "sum: $sum" | d.text color=black at=1,5
d.legend "$MAP" range=1300,1500
}

###
d.mon x2
d.resize w=1024 h=768
d.split.frame 4
#
d.frame uno
MAP=rad_test.355.beam.try5
plot_stuff "Std."
#
d.frame dos
MAP=rad_test.355.beam.Hz1deg.try5
plot_stuff "w/ 1 degree horizon seeds"
#
d.frame tres
MAP=rad_test.355.beam.Hz15deg.try5
plot_stuff "w/ 15 degree horizon seeds"
#
d.frame cuatro
MAP=rad_test.355.beam.Hz30deg.try5
plot_stuff "w/ 30 degree horizon seeds"
###
d.frame full_screen
#d.out.file out=rsun_horizons format=png

next I will try std and 1deg horizons with a r.sun time step of 3 min instead of 30 minutes. (sun travels ~1deg of sky in 4min)

Hamish

[hamish]

showing effects (/errors) when using different number of horizon slices

[hamish]

first four more step sizes to add to previous 1,15,30 deg.

with 120x 3 deg horizon maps, 'r.sun -s' runs in:
real    0m45.418s
user    0m44.895s
sys     0m0.464s

with 72x 5 deg horizon maps, 'r.sun -s' runs in:
real    0m38.214s
user    0m37.822s
sys     0m0.368s

with 45x 8 deg horizon maps, 'r.sun -s' runs in:
real    0m34.099s
user    0m33.806s
sys     0m0.296s

with 36x 10 deg horizon maps, 'r.sun -s' runs in:
real    0m34.072s
user    0m33.666s
sys     0m0.244s

see attached rsun_horizons3_5_8_10.png

interestingly 36x 10 degree horizon maps is seemingly the most consistent. I wonder if the trick is to match the horizon angle frequency with the r.sun time step= option?

---

next I will try std and 1deg horizons with a r.sun time step of 3 min instead of 30 minutes. (sun travels ~1deg of sky in 4min)

3 minute setup:

time r.sun -s elevin=gauss day=$DAY \
   beam_rad=rad_test.355.beam.3min_step.try6 step=0.05
real    24m34.998s
user    24m32.680s
sys     0m2.488s


time r.sun -s elevin=gauss day=$DAY \
   beam_rad=rad_test.355.beam.3min_step.Hz1deg.try6 \
   step=0.05 horizon=horangle horizonstep=1
real    4m45.455s
user    4m44.054s
sys     0m1.320s

plotting setup:

plot_stuff()
{
r.colors $MAP -e color=bcyr
d.rast "$MAP"
d.vect gauss_200m_contours color=180:180:180
echo "$1" | d.text color=black
eval `r.univar -g $MAP`
echo "sum: $sum" | d.text color=black at=1,5
d.legend "$MAP" range=1300,1500
}

###
d.mon x2
d.resize w=1024 h=768
d.split.frame 4
#
d.frame uno
MAP=rad_test.355.beam
plot_stuff "dt=30min; no horizon seeds"
#
d.frame dos
MAP=rad_test.355.beam.Hz1deg.try5
plot_stuff "dt=30min; 1 degree horizon seeds"
#
d.frame tres
MAP=rad_test.355.beam.3min_step.try6
plot_stuff "dt=3min; no horizon seeds"
#
d.frame cuatro
MAP=rad_test.355.beam.3min_step.Hz1deg.try6
plot_stuff "dt=3min; 1 degree horizon seeds"
###
d.frame full_screen
#d.out.file out=rsun_horizon_dt format=png

see attached rsun_horizon_dt.png.

hmmm... are the lower lobes for 1 degree r.horizon seeds real? maybe that 10 degree step size is actually the worst case. Again this is for the winter solstice so the sun should never get into the NE,NW.

I suspect the standard r.sun @ 3min steps without horizon seeds is in fact the most accurate result.

I suppose the way to test that is to run r.sun in Mode 1 for a loop over the day and then run r.series to sum them and see if the pattern matches. But that's a job for later or someone else, it's the end of the day here. I'll set it up to run overnight with a 0.0025 timestep (9sec) & no horizon seeds.

hoping this is useful, Hamish

[hamish]

r.sun -s fed with r.horizon seeds a 3,5,8, and 10 degree steps

[hamish]

r.sun -s testing 3min vs 30min time steps & 1deg horizon seeds vs none

[hamish]

I built and tried the old r.sun, there slopein= and aspin= are required and it has the same problems as rsun_aspslp.png's lower right pane.

---

I ran the test case with step=0.0025 (9 sec) and step=0.01 (36 sec). They, together with a differences map between them are in the attached pic rsun36_9sec_diff.png. Note the inner eastern wing in the 9sec image has a rather crisp edge- it would appear to be missing some near-sunset iterations??

setup:

time r.sun -s elevin=gauss day=$DAY \
   beam_rad=rad_test.355.beam.9sec_step.try6 step=0.0025
real    479m0.940s
user    476m11.914s
sys     1m12.181s


time r.sun -s elevin=gauss day=$DAY \
   beam_rad=rad_test.355.beam.36sec_step.try7 step=0.01
real    217m34.835s
user    217m12.742s
sys     0m23.561s


r.colors -e col=bcyr map=...

r.mapcalc diff36_9sec = \
   "rad_test.355.beam.36sec_step.try7 - rad_test.355.beam.9sec_step.try6"

r.colors.stddev -z diff36_9sec

plot_stuff()
{
r.colors $MAP -e color=bcyr
d.rast "$MAP"
d.vect gauss_200m_contours color=180:180:180
echo "$1" | d.text color=black
eval `r.univar -g $MAP`
echo "sum: $sum" | d.text color=black at=1,5
d.legend "$MAP" range=1300,1500
}

###
d.mon x2
d.resize w=1024 h=768
d.split.frame 4
#
d.frame uno
MAP=rad_test.355.beam.9sec_step.try6
plot_stuff "dt=0.0025 hr (9 sec), day 355"
#
d.frame dos
MAP=rad_test.355.beam.36sec_step.try7
plot_stuff "dt=0.01 hr (36 sec), day 355"
#
d.frame tres
MAP=diff36_9sec
d.rast $MAP
d.legend $MAP
echo "difference between dt=9,36 sec maps" | d.text color=black
#
d.frame cuatro
MAP=rad_test.355.beam.3min_step.try6
plot_stuff "dt=0.05 hr (3 minutes), day 355"
###
d.frame full_screen
#d.out.file out=rsun36_9sec_diff format=png

Hamish

[hamish]

pic showing differences between runs with 9 and 36 sec step size

[hamish]

r.sun2 seems safe to background without worrying if it will clobber the region settings, so a poor-man's multi-threading script is used to help speed things up on the quad-core.

The next tests are summing a whole bunch of Mode 1 steps with r.series + r.mapcalc and see how it compares to the Mode 2 daily sum map. Trying with a time step of 0.01 from sunrise to sunset on day 355:

r.info from Mode 2 beam_rad map:
Sunrise time min-max (hr.):     7.68
Sunset time min-max (hr.):     16.32

### mode 1 loop ###
SUNRISE=7.67
SUNSET=16.33
STEP=0.01
# | wc -l   867
CORES=4

for TIME in `seq $SUNRISE $STEP $SUNSET` ; do
   echo "time=$TIME"
   CMD="r.sun -s elevin=gauss day=$DAY time=$TIME \
         beam_rad=rad1_test.355_${TIME}_beam --quiet"

   # poor man's multi-threading for a multi-core CPU
   MODULUS=`echo "$TIME $STEP $CORES" | awk '{print $1 % ($2 * $3)}'`
   if [ "$MODULUS" = "$STEP" ] ; then
      # stall to let the background jobs finish
      $CMD
      sleep 2
   else
     $CMD &
   fi

done


for TIME in `seq $SUNRISE $STEP $SUNSET` ; do
   r.colors rad1_test.355_${TIME}_beam col=rules --quiet << EOF
     0 white
     0.000001 blue
     33.333% cyan
     66.667% yellow
     100% red
EOF
done

... in progress ...

Hamish

[hamish]

... continuing on from the last entry in the bug log ...

# combine into hourly sums (867 input maps too many for shell)
for HOUR in `seq 7 16` ; do
   echo " $HOUR o'clock"
   INMAPS=`g.mlist rast patt="rad1_test.${DAY}_$HOUR.??_beam" sep=,`
   r.series in="$INMAPS" out=beam_sum_$HOUR.oclock method=sum
   r.colors beam_sum_$HOUR.oclock color=bcyr
done

# daily
INMAPS=""
for HOUR in `seq 7 16` ; do
   INMAPS="$INMAPS,beam_sum_$HOUR.oclock"
done
INMAPS=`echo "$INMAPS" | sed -e 's/^,//'`

# animate the day
xganim $INMAPS

# daily sum
r.series in="$INMAPS" out=beam_sum_day.355 method=sum
r.mapcalc "beam_irad_day.355 = beam_sum_day.355 * $STEP"

the result is highly similar (but not exactly identical numerically) to the same step size (36sec; dt=0.01) Mode 2 raster as seen in the upper right pane of the rsun36_9sec_diff.jpg attached image. So choosing mode 1 vs mode 2 seems to introduce no errors.

# animate an hour
HOUR=8
INMAPS=`g.mlist rast patt="rad1_test.355_$HOUR.??_beam" sep=,`
xganim $INMAPS

here we see a big jump between runs for time = (8.45,8.46) (8.72,8.73) (8.82,8.83).

Just looking at the biggest lobe-jump (8.45,8.46) I did a few more runs at time=8.455, 8.45875, etc. Results in rsun_8oclock_jump.png (attached) detailing the moment of the sudden disappearance of a lobe on the south-eastern side. These lobes seem to be responsible for the jumpy nature of the Mode 2 output.

I looked at the input map's (gauss) slope and curvature maps (1st & 2nd derivative) from r.slope.aspect; the input map seems perfectly smooth & error free.

Even with this, it should be noted that the output is much cleaner than the old version of r.sun (as long as you don't use the aspin= and slopein= options, then it still hits the same bug as the old version)

??, Hamish

[hamish]

showing transition of sudden jump between 8.45 and 8.46

[hamish]

Hi,

lat=, latin=, and longin= now play nice with pj_proj() in devbr6 and trunk.

• lat= can no longer work due to rsunlib's needs, so I ripped it out. It's kind of lame anyway.
• latin= and longin= need to be given together (except with civiltime=)
• even if you give those seed raster maps it only gives you about a 1.3% speed increase (4 sec over 4.5minutes for my single test). The calc in rsunlib needs Cartesian coords so it only really saves the pj_proj() calc in one place, and that really isn't so expensive anyway. The man page says those lat= opts are there in case you want to run r.sun from a simple XY location. With the new pj_proj() stuff in rsunlib I don't think that's viable anymore, and even if it did work it's not worth the confusion it introduces. Given all that I'm inclined to rip out those options and have civiltime= simply use pj_proj() if it wants to know the local longitude.

see r38770,2.

Hamish

[neteler]

Replying to hamish:

correction: using aspin=r.slope.aspect doesn't speed things up, but slope+aspect does. Using slopein= causes differences (missing south west butterfly wing), but slope+aspect makes it really bad.

Let me suggest to apply the same "This does nothing" approach from lat= also to slopein= and aspin= since most users won't find the online hints easily.

[hamish]

Replying to neteler:

Let me suggest to apply the same "This does nothing" approach from lat= also to slopein= and aspin= since most users won't find the online hints easily.

but they do do something, when used together in the tests in comment:8 they cut processing time in half. (!for that particular case!)

the trouble, and where IMHO we should concentrate our effort, was that the output was buggy -- missing its last iteration(s) of the day or so.

?, Hamish

[dylan]

Has any progress been made on this topic? I recently re-computed an annual global radiance calculation with GRASS7 (previously with GRASS65), using only the elevation input (previously used elev, slope, aspect). The results were astounding. Zooming into a smaller region and checking, it appears as if the reflected radiance raster was all 0s and the diffuse radiance raster was a vertical gradient across the region-- with no apparent influence from the local terrain. When using slope and aspect inputs the output appears to be correct.

I'll work up examples tonight and re-post with figures, code, and the DEM in question.

[hamish]

Replying to dylan:

Has any progress been made on this topic?

a couple of weeks ago I applied the patch which avoided using libproj in the main loop, making things faster and fixing a modern rotational bug, but I don't think that's relevant to what you describe below.

I recently re-computed an annual global radiance calculation with GRASS7 (previously with GRASS65), using only the elevation input (previously used elev, slope, aspect). The results were astounding. Zooming into a smaller region and checking, it appears as if the reflected radiance raster was all 0s and the diffuse radiance raster was a vertical gradient across the region-- with no apparent influence from the local terrain. When using slope and aspect inputs the output appears to be correct.

I'll work up examples tonight and re-post with figures, code, and the DEM in question.

this is surprising news to me. please let us know what you find out. I'd like to sort it out before applying the OpenCL GPU patches (but the autoconf ./configure changes don't have to wait).

thanks, Hamish

[dylan]

beam radiance computed with only elevation input

[dylan]

diffuse radiance computed with only elevation input

[dylan]

beam radiance computed with elevation, slope, and aspect inputs

[dylan]

diffuse radiance computed with elevation, slope, and aspect inputs

[dylan]

Here is the code used to produce the example images. Note that the beam and diffuse radiance maps created by r.sun are clearly incorrect when slope and aspect maps are not specified. The output looks correct when slope and aspect are specified. The beam and diffuse radiance maps are attached as: "beam_rad.015.jpg" / "diff_rad.015.jpg" (only elevation given), and "sa_beam_rad.015.jpg" / "sa_diff_rad.015.jpg" (slope and aspect given). The DEM used can be found here: ​http://basho.dyndns.org/~dylan/temp/elev30.tif

g.region rast=elev30 res=30 -pa
r.slope.aspect elevation=elev30 slope=slope30 aspect=aspect30

r.sun -s --o elevin=elev30  day=15  beam_rad=beam_rad.015 diff_rad=diff_rad.015 refl_rad=refl_rad.015 glob_rad=glob_rad.015

r.sun -s --o elevin=elev30 slopein=slope30 aspin=aspect30 day=15  beam_rad=sa_beam_rad.015 diff_rad=sa_diff_rad.015 glob_rad=sa_glob_rad.015 refl_rad=sa_refl_rad.015

[dylan]

Even stranger... When running r.sun (GRASS7) without shadows, and without slope/aspect input maps, the beam radiance output looks a lot like the input elevation raster...(?)!

Using the DEM referenced upstream:

r.sun --o elevin=elev30  day=15  beam_rad=beam_rad.015

[dylan]

An updated link to the DEM referenced in comment 18 can be found here: ​http://216.86.177.86/~dylan/temp/elev30.tif

[hamish]

Replying to dylan:

The output looks correct when slope and aspect are specified.

One thing I notice about sa_diff_rad.015.jpg is that the diffuse output there seems to incorporate the effects of direct shading, while diff_rad.015.jpg does not.

I would ~ guess that the amount of diffuse light reaching a grid cell would be directly proportional to the visible sky. (stable throughout the day; should be independent of aspect (sa_diff.015.jpg is correlated with aspect); perhaps we should try a new r.surf.volcano surface and see how much diffuse light gets into a deep caldera or well, with the only visible sky a small amount directly above)

Ultimately I think a formal verification of the model using (multiple) real light meters is going to be necessary, but to get the deployment design right we're going to have to become experts in the algorithm. (e.g. if a light meter is shaded by a reflective wall, its beam will be 0 but will it's diffuse be increased vs a black wall? how many bounces does reflected light follow? just one?) ...

here's another pseudo-data test-case, I'll post results once it finishes.

#spearfish (further north than NC so more defined shadows)
g.region -d
r.mapcalc "undulates = (2 + sin( row() * 2 ) + cos( col() * 2 )) * 500"
r.colors undulates color=bcyr
nviz undulates

DAY=355
time r.sun -s elevin=undulates day=$DAY step=0.05 \
  beam_rad=rad_test.und.355.beam \
  diff_rad=rad_test.und.355.diff \
  refl_rad=rad_test.und.355.refl

Hamish

[hamish]

to be continued ...

[rorschach]

Any updates on this for 7.0.0? Thanks.

[wenzeslaus]

Is somebody able to write some tests showing what is wrong and what is correct? It would be important to have these tests no not only to show what is wrong now, but also that the new potential implementation actually fixes the issue and does not break the existing functionality. Read about crating tests here:

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

[neteler]

For a test case, see also #2606

[neteler]

Milestone renamed

[neteler]

Ticket retargeted after milestone closed

[mlennert]

What is the status on this bug report ?

[neteler]

Ticket retargeted after milestone closed

[martinl]

Ticket retargeted after milestone closed

[martinl]

Still relevant?