r.mapcalc: Allow rounding of floating numbers

[pvanbosgeo]

The round() function in r.mapcalc always returns an integer, regardless of its argument types. Integers are always 32-bit, so the result is limited to the range +/- 2147483647 (2^31-1).

Suggestion: extend the function to allow to round numbers outside the integer range, and to round with a specific number of decimal places, like e.g., the function round() in R.

Various options have been discussed on the mailing list, see this ​email thread

[mmetz]

Replying to pvanbosgeo:

The round() function in r.mapcalc always returns an integer, regardless of its argument types. Integers are always 32-bit, so the result is limited to the range +/- 2147483647 (2^31-1).

Suggestion: extend the function to allow to round numbers outside the integer range, and to round with a specific number of decimal places, like e.g., the function round() in R.

Try trunk r56313.

The output type of round() is now the same like the input type. Rounding to a given number of decimal places is supported with round(x, y) with y = number of decimal places. The new function round(x, y) supports a negative number of decimal places: for example, round(119, -1) results in 120, and round(119, -2) results in 100.

[pvanbosgeo]

Works, great work!!

[glynn]

Replying to mmetz:

The output type of round() is now the same like the input type.

This changes long-standing behaviour in a way which could break scripts. E.g. if k is an integer,

round(x) / k

would always evaluate to an integer, whereas now it may result in a fraction.

Rounding to a given number of decimal places is supported with round(x, y)

From the mailing list discussion ...

Rounding to a given number of decimals is unnecessarily limiting. The algorithm for generalised rounding is essentially:

	roundTo(x, k) = round(x / k) * k.

Rounding to N decimal places is just a case of using k=1/10^N. If you allow k to be specified directly, then you can round to arbitrary steps (e.g. k=5 would round to the nearest multiple of 5, etc).

However: there's a slight problem with doing it that way: 0.1 isn't exactly representable in binary, so e.g. x/0.1 isn't equal to x*10; it would be more accurate to use:

	roundTo(x, k) = round(x * k) / k

where k is the reciprocal of the step, so k=10^N to round to N decimal places (or k=2 to round to 1/2).

The downside is that the interface is less useful if you want to round to something other than a fixed number of decimal places. E.g. if you wanted to round to the nearest multiple of 45 degrees, you'd need to use k=1.0/45, which isn't exactly representable.

Unless someone has a better idea, I plan to change the round() function so that the second argument is the step value, and add an roundi() (round-inverse) function where the second argument is the reciprocal of the step value (to avoid the rounding error when using a step of 10^-N).

[mmetz]

Replying to glynn:

Replying to mmetz:

The output type of round() is now the same like the input type.

This changes long-standing behaviour in a way which could break scripts.

OK. Restoring the original behaviour for round(x) is easy. But it would be nice to have a round(x, y) function that preserves the data type of x in order to have a possibility to avoid integer overflow.

Rounding to a given number of decimal places is supported with round(x, y)

From the mailing list discussion ...

Rounding to a given number of decimals is unnecessarily limiting. The algorithm for generalised rounding is essentially:

	roundTo(x, k) = round(x / k) * k.

Rounding to N decimal places is just a case of using k=1/10^N. If you allow k to be specified directly, then you can round to arbitrary steps (e.g. k=5 would round to the nearest multiple of 5, etc).

I was just looking at the function round() in R which rounds to decimal places. Generalised rounding makes more sense.

However: there's a slight problem with doing it that way: 0.1 isn't exactly representable in binary, so e.g. x/0.1 isn't equal to x*10; it would be more accurate to use:

	roundTo(x, k) = round(x * k) / k

where k is the reciprocal of the step, so k=10^N to round to N decimal places (or k=2 to round to 1/2).

Unless someone has a better idea, I plan to change the round() function so that the second argument is the step value, and add an roundi() (round-inverse) function where the second argument is the reciprocal of the step value (to avoid the rounding error when using a step of 10^-N).

Sounds good to me.

[glynn]

Replying to glynn:

Unless someone has a better idea, I plan to change the round() function so that the second argument is the step value,

Done in r56365. An optional third argument is the start value, so e.g. round(x,1.0,0.5) will round to the nearest something-point-five value.

and add an roundi() (round-inverse) function where the second argument is the reciprocal of the step value (to avoid the rounding error when using a step of 10^-N).

I haven't bothered with this. If the step value can't be represented exactly, then in the general case, neither can the rounded value. If it's desired, it would be better to clone xround.c and modify the i_round() function (swap the multiplication and division) than to try to get yet another case into that file.

Rounding to a given number of decimal places (as opposed to rounding to a multiple of 10^-N) is something which really needs to be done during conversion to a string. Attempting to round a floating-point value to a given number of decimal places will inevitably add rounding errors which may be visible if the value is subsequently converted to a string using sufficient precision.

[pvanbosgeo]

Hi Glynn, this looks and works great. Initially I wasn't really clear about the implementation, but the round(x,y) seems to do exactly what I hoped for (including rounding numbers outside the integer range), great!

The explanations for round(x,y) and round(x,y,z) in the r.mapcalc help file are perhaps not immediately clear (but that might be me). Maybe it is possible to get some examples in the help file? In any case, great work.

[neteler]

Replying to pvanbosgeo:

The explanations for round(x,y) and round(x,y,z) in the r.mapcalc help file are perhaps not immediately clear (but that might be me). Maybe it is possible to get some examples in the help file?

An example for each round(x,y) and round(x,y,z) would be appreciated (if here, I am happy to add them the manual).

E.g. convert degree Celsius map/floating point to 10*degC (like BIOCLIM) as integer?

[martinl]

Milestone renamed

[martinl]

What is status of this issue?

[neteler]

Ticket retargeted after milestone closed

[Nikos Alexandris]

I also think that two simple examples in the manual would be useful in understanding the use of round(x,y) and round(x,y,z).

I tested locally adding

split_window_expression = '(round({swe}, 2, 0.5))'.format(swe=split_window_expression)

in ​https://gitlab.com/NikosAlexandris/i.landsat8.swlst/blob/master/i.landsat8.swlst.py after a bit of trial and error.

[martinl]

What is the state of the ticket?

[martinl]

All enhancement tickets should be assigned to 7.6 milestone.

[martinl]

Ticket retargeted after milestone closed

[martinl]

Ticket retargeted after milestone closed