#610: r.sun: incidout uses 0 instead of NULL
--------------------+-------------------------------------------------------
Reporter: hamish | Owner: grass-dev@lists.osgeo.org
Type: defect | Status: new
Priority: normal | Milestone: 6.5.0
Component: Raster | Version: 6.4.0 RCs
Keywords: r.sun | Platform: All
Cpu: All |
--------------------+-------------------------------------------------------
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.
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:
#610: r.sun: incidout uses 0 instead of NULL
---------------------+------------------------------------------------------
Reporter: hamish | Owner: grass-dev@lists.osgeo.org
Type: defect | Status: new
Priority: normal | Milestone: 6.5.0
Component: Raster | Version: 6.4.0 RCs
Resolution: | Keywords: r.sun
Platform: All | Cpu: All
---------------------+------------------------------------------------------
Comment (by 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.
#610: r.sun: incidout uses 0 instead of NULL
---------------------+------------------------------------------------------
Reporter: hamish | Owner: grass-dev@lists.osgeo.org
Type: defect | Status: new
Priority: normal | Milestone: 6.5.0
Component: Raster | Version: 6.4.0 RCs
Resolution: | Keywords: r.sun
Platform: All | Cpu: All
---------------------+------------------------------------------------------
Comment (by mmetz):
Replying to [comment:6 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 ?
#610: r.sun: incidout uses 0 instead of NULL
---------------------+------------------------------------------------------
Reporter: hamish | Owner: grass-dev@lists.osgeo.org
Type: defect | Status: new
Priority: normal | Milestone: 6.5.0
Component: Raster | Version: 6.4.0 RCs
Resolution: | Keywords: r.sun
Platform: All | Cpu: All
---------------------+------------------------------------------------------
Comment (by hamish):
Replying to [comment:8 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.