[GRASS-dev] Re: r.series: threshold/count method?

[Markus Neteler]

On Mon, Sep 7, 2009 at 11:49 PM, Glynn Clements<glynn at gclements.plus.com> wrote:
> Markus Neteler wrote:
>
>> I have an extra issue: For the particular growing degree day calculation
>> I am facing the problem that the values I am comparing with a FP values
>> and that GRASS_EPSILON is useless in this particular case. It works
>> only with
>>
>>  /* EPSILON = GRASS_EPSILON; */
>>   EPSILON = 10.
>>   if (fabs(tval - values[i]) < EPSILON ) {
>>
>> because the pixel values are growing so quickly from one map to the
>> next. Since we need a generic solution, another "precision" or whatever
>> parameter needs to be passed otherwise the thresholding fails completely.
>> EPSILON however depends on the input maps and needs to be user
>> controllable.
>>
>> How to deal with that?
>
> What exactly are you trying to measure?
>
> From your original post, it sounds like you have a monotonically
> increasing input,

Here an example (doy = day of year; gdd = growing degree days):

"doy","gdd"
1,0.0
2,0.0
3,0.0
4,0.0
...
44,3.4
45,3.4
46,4.5
47,6.3
...
99,164.6
100,170.3
101,175.4
102,178.5
...
193,997.8
194,1008.4
195,1018.3
196,1026.9
...
363,1980.1
364,1980.1
365,1980.1

it is a kind of sigmoidal function which differs from year to year and
pixel stack to pixel stack.

> and you want the time (index) at which the value
> first exceeds the threshold, e.g.:
>
>        for (i = 0; i < n; i++) {
>            if (values[i] >= threshold) {
>                *result = (DCELL) i;
>                return;
>            }
>        }
>        Rast_set_d_null_value(result, 1);
>
> If you need an approximate equality comparison, the usual solution is
> to use the ratio of the difference to one of the values, e.g.:
>
>        double EPSILON = 1e-6; /* 1 ppm */
>        if (fabs(values[i] - threshold) < EPSILON * threshold)
>                ...

The problem is the non-linear increase here, I guess: If I would check
for 1000 GDD in above example, I have the input

...
193,997.8
194,1008.4
...

but get
EPSILON = 1e-6
> abs(997 - 1000) < EPSILON * 1000
[1] FALSE
> abs(1008.4 - 1000) < EPSILON * 1000
[1] FALSE

and still won't find it.

Markus