Incorrect beam radiation for some steep inclined surfaces

[andrewg-cse]

trafficstars

Issue was that direct beam radiation falling on the back of some steep approx north facing (in northern hemisphere) inclined surfaces was included in direct beam results. This was caused by an incorrect value of sunSlopeGeom->longit_l in the rsunlib.c, lumcline2 function. The root cause was: In main.c, in the calculate() function at line 1953, the value of q1 changes sign as the slope increases for north (in northern hemisphere) facing surfaces. This causes tan_lam_l to also change sign and hence the result of atan stored in sunSlopeGeom.longit_l to jump from just less than pi/2 to just more than -pi/2. The required value is pi more than this (which is another correct atan result as the tan plot repeats every pi radians). The fix detects this situation and adds pi to the result. I've done this addition in the rsunlib.c, lumcline2 function, as this is where the time offset is used. Also, there is a possibility of a divide by zero error if q1 is zero. This change also corrects this.

[andrewg-cse]

Well, re the test case using an artificial raster I did look at doing that, and also a C test aka lib/raster3d. The first was a bit far from the actual code but could work for the first fix, but the divide by zero will be messy due to rounding, and the second was going need more re-structuring of the code in r.sun than seemed reasonable. So I ended up writing some unit tests in Java against a cut & pasted copy of the lines I changed which at least gave me confidence, not appropriate for automated tests here though.

[petrasovaa]

I think comparing q1 != 0.0 is not going to work for float numbers, usually you would compare it with a threshold. Please see how GRASS_EPSILON is used.

To help with the review, a simple test case focusing on the first fix would be really helpful.

[marisn]

I think comparing q1 != 0.0 is not going to work for float numbers, usually you would compare it with a threshold. Please see how GRASS_EPSILON is used.

Usually you would be right, but this is not the case as q1 != 0.0 here just guards against division by zero. Thus anything not being 0.0 would be fine to be used for division.

To help with the review, a simple test case focusing on the first fix would be really helpful.

@andrewg-cse a simple "synthetic dem in, check output" test would suffice. No need to write C test to trigger div by zero case explicitly. Of course, if it is too much for you, I don't see a huge issue in merging this PR without tests.

[andrewg-cse]

@marisn yes, re q1 != 0.0, that's exactly right, it just guards against when q1 really is exactly 0.0

• I'll try and find some time to write the test case using test rasters etc

[petrasovaa]

@marisn yes, re q1 != 0.0, that's exactly right, it just guards against when q1 really is exactly 0.0

• I'll try and find some time to write the test case using test rasters etc

That would be great. Artificial rasters would be a good idea. Also adding some comments in the code could help, right now some of the new variables have unclear meanings.

[wenzeslaus]

I have updated to code to the current indent style using the steps from #2544 which resolves the conflict.

[andrewg-cse]

I've added a unit test that covers a range of horizons and uses the North Carolina test location.

[andrewg-cse]

Ok, so the tests I wrote fail when run from the CI system and I don't know why. I'm running them from the grass cli with the current mapset set to one I created under the "nc_basic_spm_grass7" dataset - so will have the same projection Then I go into the raster/r.sun/testsuite and run "python test_rsun.py" ... and all 7 tests pass. I'm using Python 3.9.12

The error from CI just says "AssertionError: The actual maximum (0.0004882812) is greater than the reference one (1e-06) for raster map test_rsun_incorrect_beam_rad_fix_diffs_rad (with minimum 0)" which is saying the biggest difference between the raster that was created and the one I was expecting is too large. So my guess is that the region it's using is somehow bigger and it's not running at the same resolution as when I generated the expected data - as I had this issue before I set the region. But I'm now setting the region on lines 363 - 365 of test_rsun.py to match one of my test rasters. So does anyone have any suggestions? Thanks.

[petrasovaa]

Not sure, when I run it locally, the test fails for me as well. Some notes though:

1.This change seems to be adding compiler warnings:

main.c:1831:9: warning: ‘isBestAM’ may be used uninitialized in this function [-Wmaybe-uninitialized]
main.c:1831:9: warning: ‘shouldBeBestAM’ may be used uninitialized in this function [-Wmaybe-uninitialized]

That's not the problem I guess, but it's worth fixing this. There are similar warnings coming from the existing code, but this PR at least shouldn't add new ones.

2. The test seems little overcomplicated and difficult to interpret. The layers seems to be reprojected. Couldn't this problem be replicated on a simpler surface, and derive the slope and aspect from the surface directly?

[hellik]

seems to be related: https://lists.osgeo.org/pipermail/grass-user/2023-January/083131.html

[hellik]

anyone a hint how to rebase this PR to main? ;-)

[ninsbl]

anyone a hint how to rebase this PR to main? ;-)

E.g: https://github.com/OSGeo/grass/pull/905#issuecomment-1435788300

[nilason]

Rebased, this is now up-to-date with main.

[hellik]

some follow up test with real data of mountain area with high geomorphology variability:

https://lists.osgeo.org/pipermail/grass-user/2023-March/083153.html

[hellik]

see also: https://lists.osgeo.org/pipermail/grass-user/2023-March/083154.html

[hellik]

real data example

g.region -p                                                                     
projection: 99 (MGI / Austria GK West)
zone:       0
datum:      hermannskogel
ellipsoid:  bessel
north:      218281.5
south:      204874.5
west:       49389.5
east:       62610.5
nsres:      1
ewres:      1
rows:       13407
cols:       13221
cells:      177253947

r.horizon elevation="laser_dgm@data" step=30 start=0.0 end=360.0 bufferzone=500 maxdistance=3000 output="horangle" distance=1.0 file="-"

r.sun cmd not patched

r.sun elevation=laser_dgm@data aspect=laser_dgm_aspect@data slope=laser_dgm_slope@data horizon_basename=horangle horizon_step=30 beam_rad=beam_rad_166_not_patched day=166

r.sun cmd patched

r.sun elevation=laser_dgm@data aspect=laser_dgm_aspect@data slope=laser_dgm_slope@data horizon_basename=horangle horizon_step=30 beam_rad=beam_rad_166 day=166

mountain summit with steep slopes in all aspect directions:

aspect

slope

difference: r.sun patched minus r.sun not patched

[hellik]

any thoughts?

[hellik]

left r.sun patched GRASS 8.3.dev, right r.sun GRASS 8.2.1 unpatched

with the patch/PR: not much differences on southern oriented mountain slopes, though steeper northern oriented slopes

[petrasovaa]

I decided to do a test on an artificial surface (saddle shape) and the results look pretty convincing.

Top row: with this PR Bottom row: before this PR Left column: without shadowing effect of terrain (-p) Right column: with shadow effect of terrain

Elevation and slope:

I can write a simpler test based on this experiment, the one in this PR is problematic (complex and didn't run).

[petrasovaa]

@andrewg-cse thank you for this important contribution!

[girishnand]

@petrasovaa @wenzeslaus

Is this patch included in current developmental versions of GRASS 8.3.0 available here - https://wingrass.fsv.cvut.cz/grass83/ ?

Or is there an easy way to apply it to my current 8.2.1 installation on Windows 11?

regards, Girish.

[petrasovaa]

@petrasovaa @wenzeslaus

Is this patch included in current developmental versions of GRASS 8.3.0 available here - https://wingrass.fsv.cvut.cz/grass83/ ?

Yes!

[girishnand]

@petrasovaa @wenzeslaus Is this patch included in current developmental versions of GRASS 8.3.0 available here - https://wingrass.fsv.cvut.cz/grass83/ ?

Yes!

Thanks. Girish.

[girishnand]

@petrasovaa

I downloaded and installed https://wingrass.fsv.cvut.cz/grass83/WinGRASS-8.3.dev-c05cc14387-42-Setup.exe But r.sun from this version still generates output identical to version 8.2.1.

This comment above -

@wenzeslaus wenzeslaus modified the milestones: 8.4.0, 8.3.0 [on Mar 23] (https://github.com/OSGeo/grass/pull/2534#event-8827645503)

suggests that the patch will perhaps be included in version 8.4.0 ?

regards, Girish.

[girishnand]

My elevation map consists of a small square hillock. The central zone is a flat square, with regions towards its north, south, east and west sloping downwards at 15 degrees, respectively towards north, south, east and west. Location Lat-Long are - Latitude : 34.179N & Longitude : -118.556. The location/projection I used is UTM zone 11N, specified as EPSG:32611. The image below shows total global insolation for the month of May using r.sun.daily

    gs.run_command('r.sun.daily', verbose=True, elevation='hotraster', 
                    horizon_basename='horangle', horizon_step='05', 
                    slope='slope', aspect='aspect', lat='hotlatitude', 
                    long='hotlongitude', step='0.25', start_day='121', 
                    end_day='151', glob_rad='MayHor05') 

Slope and Aspect are from r.slope.aspect and lat-long are from r.latlong. Horizon rasters are from r.horizon with a step of 5 degrees.

In the image, r.sun appears to show significantly greater global insolation on the north facing slope as compared to the south facing slope. For a location at latitude 34 deg N, this should not be the case.

Would appreciate if someone can try to replicate my results and point out what I am doing wrong.

regards, Girish.

PS: My input elevation raster is: hot.txt. It is in Arc/Info ASCII Grid format, so you need to use r.in.gdal to read it into GRASS GIS - gs.run_command('r.in.gdal', flags='o', input='hot.txt', output='hotraster')

[nilason]

@petrasovaa

I downloaded and installed https://wingrass.fsv.cvut.cz/grass83/WinGRASS-8.3.dev-c05cc14387-42-Setup.exe But r.sun from this version still generates output identical to version 8.2.1.

This comment above -

@wenzeslaus wenzeslaus modified the milestones: 8.4.0, 8.3.0 [on Mar 23] (#2534 (comment))

suggests that the patch will perhaps be included in version 8.4.0 ?

regards, Girish.

This PR was merged in March into main https://github.com/OSGeo/grass/commit/aac1c5a74ee5512d031132aebadee8aebf7c913f, that is 8.3.dev, the version you tried is https://github.com/OSGeo/grass/commit/c05cc143875f15234a45a231bb1efe6d47ba1e6a (same branch) from a couple of days ago. In short, this change is included in the version you tried.

[girishnand]

@nilason

Thanks.

Let me try out a few other elevation rasters with different slopes to see if I can get different results for v8.3dev (c05cc14) versus v8.2.1.

regards, Girish.

[petrasovaa]

The change from this PR would likely not influence your results, the slopes would have to be much higher.

But I am getting different results:

r.sun.daily elevation=hot aspect=aspect slope=slope start_day=121 end_day=151 glob_rad=global nprocs=10

[girishnand]

@petrasovaa

Appreciate you running my test case.

I get the same results as you do if I just run r.sun.daily. But this seems to change if r.horizon is used.

Can you please try the following sequence of commands -

r.horizon elevation=hotraster step=05 output=horzangle distance=0.5 r.sun.daily elevation=hotraster horizon_basename=horzangle horizon_step=05 slope=slope aspect=aspect start_day=121 end_day=151 glob_rad=global

Or is there something wrong with this?

regards, Girish.

This is the output I get -

[petrasovaa]

I get the same results as you do if I just run r.sun.daily. But this seems to change if r.horizon is used.

Can you please try the following sequence of commands -

r.horizon elevation=hotraster step=05 output=horzangle distance=0.5 r.sun.daily elevation=hotraster horizon_basename=horzangle horizon_step=05 slope=slope aspect=aspect start_day=121 end_day=151 glob_rad=global

Or is there something wrong with this?

When I compare r.sun output with and without using horizons, I get significantly different values for some areas. I don't have a good explanation, so perhaps create a new issue and show how the results of r.sun differ with and without horizons. r.sun.daily is just a wrapper around r.sun.