[GRASS-dev] r.univar: 'nan' results

[Glynn Clements]

Maciej Sieczka wrote:

> >> (nan = not a number, GRASS uses NULL for the same idea, but NULL is not
> >> the same as nan.)

FWIW, the bit patterns which GRASS uses to represent null values for
FCELL and DCELL (all-ones) are interpreted as NaN values when using
the IEEE-754 floating-point format. However, the converse isn't
necessarily true; IEEE-754 defines any value with an all-ones exponent
and a non-zero mantissa as NaN. IOW, GRASS' null values are IEEE-754
NaN values, but they aren't the only ones.

In particular, the NaN value resulting from e.g. 0.0/0 or domain
errors in math functions won't normally be the GRASS null value.

> > agreed, for me, it seems like a bug (v.to.rast?), in any case a 'nan'
> > value is an *invalid* cell value.
> 
> >> If there are nan in the cell data, the r.univar calculations can't cope.
> >> The question is: how did nan cells get into the input raster map?
> 
> > v.to.rast ?
> 
> >> There is a possibility that there is a bug in r.univar..
> >> Try using r.mapcalc to isolate the nans? What does "r.info -r" say?
> 
> > min=0.000000
> > max=358.264295
> 
> It seems like a bug in v.to.rast. Take the location I attached, and do:
> 
> $ g.region vect=cieki10_2zl res=10 -a
> $ v.to.rast input=cieki10_2zl output=cieki10_2zl_dir use=dir
> 
> There are several nan cells indeed:

Without having tried the data, my first guess is that the NaN values
result from a domain error in the flow-direction calculation,
specifically from the call to acos() in v2angle(), in do_lines.c.

The main problem is that you can't pass GRASS FP null values around
like normal values. When dealing with NaN values, the implementation
isn't required to maintain the specific bit pattern, only that
"NaN-ness" is preserved.

It would be easy enough to add a check to v2angle(), e.g.

	if (costheta < -1 || costheta > 1)
	{
		/* acos() will return NaN */
	}

However, there's no way to report this; we can't just return GRASS'
null value, as the hardware may transform it to a different
bit-pattern en route.

We could try to detect the NaNs before the data is written to file,
but isnan() isn't in C89, and other ways of detecting NaNs are
unreliable. ANSI dictates that "x == x" is false if x is NaN, but some
compilers (particularly for architectures without full hardware FP
support, e.g. Alpha) overlook this (and others may optimise "x == x"
to true with certain optimisation switches).

The only robust and portable solution is to change v2angle(),
deg_angle() and set_dcat() to accept/return the value via a pointer
argument, to ensure that the all-ones bit pattern is preserved, and
change v2angle() to specifically check for out-of-domain values.

-- 
Glynn Clements <glynn at gclements.plus.com>