r.sun: incidout uses 0 instead of NULL

[hamish]

Hi,

r.sun(2)'s incidout= output map uses 0 instead of NULL for shadowed areas which have no direct line of sight to the sun. This harms any sort of numerical analysis on the output because it biases towards zero, and makes the creation of MASKs harder.

Hamish

[neteler]

... but it uses G_set_f_null_value(): ​http://trac.osgeo.org/grass/browser/grass/trunk/raster/r.sun2/main.c#L1221

[mmetz]

... but what about ​http://trac.osgeo.org/grass/browser/grass/trunk/raster/r.sun2/rsunlib.c#L439

shouldn't that be

return UNDEFZ;

? A wild guess.

[hamish]

test case:

G65> g.region -d  # spearfish: 44:30N  477x634
G65> r.surf.volcano out=gauss method=gaussian kurtosis=1

just after sunrise on a frosty winter's morn

G65> r.sun -s elevin="gauss" incidout="day355rad.inci.08hr" \
  day=355 time=8.00

note western boundary cells are all NULL (!??) and no other:

G65> r.univar day355rad.inci.08hr
 100%
total null and non-null cells: 302418
total null cells: 477

Of the non-null cells:
----------------------
n: 301941
minimum: 0
maximum: 2.85292
range: 2.85292
mean: 2.35481
mean of absolute values: 2.35481
standard deviation: 1.04271
variance: 1.08725
variation coefficient: 44.2802 %
sum: 711013.3009517193

real data range is 2.7-2.8deg above horizon note large number of 0.0 cells:

G65> r.stats day355rad.inci.08hr nsteps=10 -c
 100%
0-0.011188 49481
2.752233-2.763421 695
2.763421-2.774609 33611
2.774609-2.785797 29075
2.785797-2.796985 27769
2.796985-2.808173 34948
2.808173-2.819361 42280
2.819361-2.830549 43446
2.830549-2.841737 40635
2.841737-2.852925 1
* 477

work-around:

r.null day355rad.inci.08hr setnull=0.0

then you get nice d.histogram output, etc.

mmetz wrote:

r.sun2/rsunlib.c#L439 shouldn't that be return UNDEFZ; ? A wild guess.

no luck; same result.

Hamish

[hamish]

correction, the above r.stats output is the result of:

G65> r.stats day355rad.inci.08hr nsteps=255 -c

range bins containing no cells are not printed in the output (even though -n,-N flags are not given)

Hamish

[mmetz]

Next try. IIUR, incidout is calculated in rsunlib.c by lumcline(). For a shadowed area, sunVarGeom->isShadow is apparently set to 1 and incidence is left as 0. For an area where there is still night, incidence angle is negative and set to 0 (r.sun2/rsunlib.c#L439). That's why my first wild guess didn't work, it was probably setting night to NULL, not shadow. If any of the input maps is NULL, incidence (as all other output maps) is set to NULL (says the manual). To keep a maximum of information in the incidence output map, you could e.g. set incidence of shadows to -1, make night -2, keep special case incidence == 0 with no shadow, and propagate any input NULL cells as NULL cells. Then optionally use r.null to set shadows and/or night to NULL. Or discard all this info and make input NULL cells, night, shadow and special case incidence == 0 with no shadow all NULL (if (s <= 0) return UNDEFZ;)

Not tested, all interpreted from the code and the manual.

Markus M

[neteler]

I contacted Jaro:

On Fri, May 22, 2009 at 11:09, Jaro Hofierka wrote:

... I believe that it is in rsunlib.c in lumcline2() function. If s < 0 (shadow) it returns 0, what in fact should be null.

Also please note that the s variable in lumcline() is initiated to 0.

    double s = 0;

(this covers also other causes of shadows)

So what is returned from the lumcline() as 0 should be written as null in the output.

[mmetz]

Replying to neteler:

I contacted Jaro:

On Fri, May 22, 2009 at 11:09, Jaro Hofierka wrote:

... I believe that it is in rsunlib.c in lumcline2() function. If s < 0 (shadow) it returns 0, what in fact should be null.

Also please note that the s variable in lumcline() is initiated to 0.

    double s = 0;

(this covers also other causes of shadows)

So what is returned from the lumcline() as 0 should be written as null in the output.

As simple as

if (s <= 0)
    return UNDEFZ;

in trunk/raster/r.sun2/rsunlib.c#L438 - #L439 ?

Markus M

[hamish]

Replying to mmetz:

As simple as

if (s <= 0)
    return UNDEFZ;

in trunk/raster/r.sun2/rsunlib.c#L438 - #L439 ?

yes, as simple as that!

now, any ideas about why the western most column is all nulls? consequence of not supplying aspect/slope maps??

Hamish

[hamish]

Replying to hamish:

now, any ideas about why the western most column is all nulls? consequence of not supplying aspect/slope maps??

weird, if I zoom in to be entirely within elevation.dem, then run with slope=slope and aspect=aspect, now the southern most row and eastern most column are not NULL. I guess that's just because those two edges have only NULL between them and the sun rising in the south-east that day, so no way of knowing that they're on a downslope or not.

Hamish

[hamish]

fixed in r38184 etc. with patch from Thomas Huld.