r.sun: Incorrect beam radiation for some steep inclined surfaces #2534

andrewg-cse · 2022-08-24T14:54:42Z

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.

…p approx north facing (in northern hemisphere) inclined surfaces is included in direct beam results. This was caused by an incorrect value of sunSlopeGeom->longit_l in the rsunlib.c, lumcline2 function. 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 fix also corrects this.

marisn

Unfortunately I can not fully judge the reasoning, code seems to be fine.
Although it is not mandatory, could you try to add a simple test case? You can take a look at r.in.pdal tests to see how to import small artificial DEM and how to compare generated output to a reference raster.

andrewg-cse · 2022-08-24T16:03:42Z

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 · 2022-08-25T08:01:22Z

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 · 2022-08-25T11:32:00Z

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 · 2022-08-25T11:43:47Z

@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 · 2022-08-26T13:48:41Z

@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.

…tive indent

wenzeslaus · 2022-09-07T08:14:26Z

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

…se/grass into r.sun_facing_away_from_sun_fix

andrewg-cse · 2022-09-12T13:32:57Z

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

andrewg-cse · 2022-09-14T16:25:03Z

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 · 2022-09-19T09:21:38Z

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.

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 · 2023-02-04T19:31:09Z

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

hellik · 2023-03-05T16:57:48Z

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

ninsbl · 2023-03-05T17:46:51Z

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

E.g: #905 (comment)

nilason · 2023-03-06T10:25:15Z

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

hellik · 2023-03-11T16:01:26Z

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 · 2023-03-11T17:53:20Z

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

hellik · 2023-03-11T18:06:50Z

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 · 2023-03-11T18:15:07Z

any thoughts?

hellik · 2023-03-11T19:43:08Z

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 · 2023-03-16T20:45:50Z

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).

…rom_sun_fix

wenzeslaus

The results are convincing. The tests are included. The code looks okay.

Let's merge this and include in 8.3.

This is a possible errata case (#2813).

petrasovaa · 2023-03-23T13:31:18Z

@andrewg-cse thank you for this important contribution!

girishnand · 2023-04-24T13:42:10Z

@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 · 2023-04-24T17:43:50Z

@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 · 2023-04-25T08:35:23Z

@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 · 2023-04-25T11:41:22Z

@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.

girishnand · 2023-04-25T11:46:50Z

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 · 2023-04-25T12:09:21Z

@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 aac1c5a, that is 8.3.dev, the version you tried is c05cc14 (same branch) from a couple of days ago.
In short, this change is included in the version you tried.

girishnand · 2023-04-25T13:26:17Z

@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 · 2023-04-25T13:38:50Z

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 · 2023-04-25T14:33:31Z

@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 · 2023-04-25T18:03:42Z

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.

girishnand · 2023-04-25T19:36:59Z

Will do, Thanks.

regards,
Girish.

girishnand · 2023-04-26T14:39:15Z

New issue created here -

#2935

regards,
Girish.

Issue was that direct beam radiation falling on the back of some steep approx north facing (in northern hemisphere) inclined surfaces is included in direct beam results. This was caused by an incorrect value of sunSlopeGeom->longit_l in the rsunlib.c, lumcline2 function. 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 fix also corrects this. Co-authored-by: Anna Petrasova <kratochanna@gmail.com>

marisn approved these changes Aug 24, 2022

View reviewed changes

andrewg-cse and others added 3 commits September 6, 2022 16:50

Additional tests for incorrect beam radiation fix

1d4338a

Indent code

aead872

Merge main with indent into r.sun_facing_away_from_sun_fix with selec…

db76e4f

…tive indent

andrewg-cse added 2 commits September 12, 2022 12:18

Altered to used compressed tifs

294afb2

Merge branch 'r.sun_facing_away_from_sun_fix' of github.com:andrewg-c…

d293e6b

…se/grass into r.sun_facing_away_from_sun_fix

nilason added this to the 8.4.0 milestone Feb 15, 2023

nilason added bug Something isn't working raster Related to raster data processing C Related code is in C labels Feb 15, 2023

nilason added 3 commits March 6, 2023 10:20

apply clang-format

bb547c4

Merge branch 'main' into r.sun_facing_away_from_sun_fix

c31c910

revert bad merge conflict fixes

ed1ac47

fix -Werror=maybe-uninitialized

f739325

petrasovaa added 3 commits March 17, 2023 16:54

rewrite test to make it simple

250552f

Merge remote-tracking branch 'upstream/main' into r.sun_facing_away_f…

5c50191

…rom_sun_fix

black and reformat long string

32dbb7b

wenzeslaus approved these changes Mar 23, 2023

View reviewed changes

petrasovaa merged commit aac1c5a into OSGeo:main Mar 23, 2023

wenzeslaus modified the milestones: 8.4.0, 8.3.0 Mar 23, 2023

wenzeslaus changed the title ~~Incorrect beam radiation for some steep inclined surfaces~~ r.sun: Incorrect beam radiation for some steep inclined surfaces Jun 6, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

r.sun: Incorrect beam radiation for some steep inclined surfaces #2534

r.sun: Incorrect beam radiation for some steep inclined surfaces #2534

andrewg-cse commented Aug 24, 2022

marisn left a comment

andrewg-cse commented Aug 24, 2022

petrasovaa commented Aug 25, 2022

marisn commented Aug 25, 2022

andrewg-cse commented Aug 25, 2022

petrasovaa commented Aug 26, 2022

wenzeslaus commented Sep 7, 2022

andrewg-cse commented Sep 12, 2022

andrewg-cse commented Sep 14, 2022

petrasovaa commented Sep 19, 2022

hellik commented Feb 4, 2023

hellik commented Mar 5, 2023

ninsbl commented Mar 5, 2023

nilason commented Mar 6, 2023

hellik commented Mar 11, 2023

hellik commented Mar 11, 2023

hellik commented Mar 11, 2023

hellik commented Mar 11, 2023

hellik commented Mar 11, 2023

petrasovaa commented Mar 16, 2023

wenzeslaus left a comment

petrasovaa commented Mar 23, 2023

girishnand commented Apr 24, 2023

petrasovaa commented Apr 24, 2023

girishnand commented Apr 25, 2023

girishnand commented Apr 25, 2023

girishnand commented Apr 25, 2023

nilason commented Apr 25, 2023

girishnand commented Apr 25, 2023

petrasovaa commented Apr 25, 2023

girishnand commented Apr 25, 2023

petrasovaa commented Apr 25, 2023

girishnand commented Apr 25, 2023

girishnand commented Apr 26, 2023

r.sun: Incorrect beam radiation for some steep inclined surfaces #2534

r.sun: Incorrect beam radiation for some steep inclined surfaces #2534

Conversation

andrewg-cse commented Aug 24, 2022

marisn left a comment

Choose a reason for hiding this comment

andrewg-cse commented Aug 24, 2022

petrasovaa commented Aug 25, 2022

marisn commented Aug 25, 2022

andrewg-cse commented Aug 25, 2022

petrasovaa commented Aug 26, 2022

wenzeslaus commented Sep 7, 2022

andrewg-cse commented Sep 12, 2022

andrewg-cse commented Sep 14, 2022

petrasovaa commented Sep 19, 2022

hellik commented Feb 4, 2023

hellik commented Mar 5, 2023

ninsbl commented Mar 5, 2023

nilason commented Mar 6, 2023

hellik commented Mar 11, 2023

hellik commented Mar 11, 2023

hellik commented Mar 11, 2023

hellik commented Mar 11, 2023

hellik commented Mar 11, 2023

petrasovaa commented Mar 16, 2023

wenzeslaus left a comment

Choose a reason for hiding this comment

petrasovaa commented Mar 23, 2023

girishnand commented Apr 24, 2023

petrasovaa commented Apr 24, 2023

girishnand commented Apr 25, 2023

girishnand commented Apr 25, 2023

girishnand commented Apr 25, 2023

nilason commented Apr 25, 2023

girishnand commented Apr 25, 2023

petrasovaa commented Apr 25, 2023

girishnand commented Apr 25, 2023

petrasovaa commented Apr 25, 2023

girishnand commented Apr 25, 2023

girishnand commented Apr 26, 2023