I am working with a raster with very large integers. That are stored as 64 bit floats (imported from a GTiff). I wish to extract part of these values with an integer division. The first approach was:
r.mapcalc ‘short_ints = int(long_ints)/1000000’
However, this results in something that is not an integer division. In fact I have no idea what the result is. Eventually, I succeeded with the following formulation:
r.mapcalc ‘psu_ids = int(ssu_ids/1000000)’
I have been going through the manual back and forth and can’t understand why these two expressions produce different results. Any insights?
čt 11. 4. 2024 v 8:26 odesílatel Luís Moreira de Sousa via grass-user
<grass-user@lists.osgeo.org> napsal:
I am working with a raster with very large integers. That are stored as 64 bit floats (imported from a GTiff). I wish to extract part of these values with an integer division. The first approach was:
r.mapcalc 'short_ints = int(long_ints)/1000000'
However, this results in something that is not an integer division. In fact I have no idea what the result is. Eventually, I succeeded with the following formulation:
r.mapcalc 'psu_ids = int(ssu_ids/1000000)'
I have been going through the manual back and forth and can't understand why these two expressions produce different results. Any insights?
This is probably due to the 32-bit size of `int()`. If your values are
higher than `2 ** 31 -1` (i.e. 214783467; `31` because half of the
numbers go to the negatives; `-1` because there is also NULL), the
counter starts to rotate. It means that `2 ** 31` will be NULL and `2
** 31 + 1` will be `- (2 ** 31 - 1)`. You can encounter this behaviour
also in other software, such as in `numpy`.
Indeed, there is no mention of that in the docs. I will update it
today to include this information. It is true that it's easy to break
your teeth on this.
On Thursday, April 11th, 2024 at 9:49 AM, Ondřej Pešek <pesej.ondrek@gmail.com> wrote:
Hi Luís,
čt 11. 4. 2024 v 8:26 odesílatel Luís Moreira de Sousa via grass-user
grass-user@lists.osgeo.org napsal:
> I am working with a raster with very large integers. That are stored as 64 bit floats (imported from a GTiff). I wish to extract part of these values with an integer division. The first approach was:
>
> r.mapcalc 'short_ints = int(long_ints)/1000000'
>
> However, this results in something that is not an integer division. In fact I have no idea what the result is. Eventually, I succeeded with the following formulation:
>
> r.mapcalc 'psu_ids = int(ssu_ids/1000000)'
>
> I have been going through the manual back and forth and can't understand why these two expressions produce different results. Any insights?
This is probably due to the 32-bit size of `int()`. If your values are
higher than `2 ** 31 -1` (i.e. 214783467; `31` because half of the
numbers go to the negatives; `-1` because there is also NULL), the
counter starts to rotate. It means that `2 ** 31` will be NULL and `2 ** 31 + 1` will be `- (2 ** 31 - 1)`. You can encounter this behaviour
also in other software, such as in `numpy`.
Indeed, there is no mention of that in the docs. I will update it
today to include this information. It is true that it's easy to break
your teeth on this.
And that WIKI is wrong. gis.h defines CELL as int thus only thing
guaranteed is that CELL will be able to store values in range from
-32767 to +32767. But by magic of C standard on my system CELL can
store values in range from -2147483647 to +2147483647 (as long as my
compiler uses /usr/include/limits.h). One running GRASS on an ARM CPU
system most likely will see different numbers.
I do agree that it is a good idea to note integer overflow in
r.mapcalc, but I would be careful with "largest number" as it has to
be tied-up to arch/compiler combo.
Māris.
ceturtd., 2024. g. 11. apr., plkst. 12:16 — lietotājs Luís Moreira de
Sousa via grass-user (<grass-user@lists.osgeo.org>) rakstīja:
Thanks Ondrej, that makes sense.
I just found the details in the Wiki [0], perhaps those should be in the manual.
On Thursday, April 11th, 2024 at 9:49 AM, Ondřej Pešek <pesej.ondrek@gmail.com> wrote:
> Hi Luís,
>
> čt 11. 4. 2024 v 8:26 odesílatel Luís Moreira de Sousa via grass-user
> grass-user@lists.osgeo.org napsal:
>
> > I am working with a raster with very large integers. That are stored as 64 bit floats (imported from a GTiff). I wish to extract part of these values with an integer division. The first approach was:
> >
> > r.mapcalc 'short_ints = int(long_ints)/1000000'
> >
> > However, this results in something that is not an integer division. In fact I have no idea what the result is. Eventually, I succeeded with the following formulation:
> >
> > r.mapcalc 'psu_ids = int(ssu_ids/1000000)'
> >
> > I have been going through the manual back and forth and can't understand why these two expressions produce different results. Any insights?
>
>
> This is probably due to the 32-bit size of `int()`. If your values are
> higher than `2 ** 31 -1` (i.e. 214783467; `31` because half of the
> numbers go to the negatives; `-1` because there is also NULL), the
> counter starts to rotate. It means that `2 ** 31` will be NULL and `2 ** 31 + 1` will be `- (2 ** 31 - 1)`. You can encounter this behaviour
> also in other software, such as in `numpy`.
>
> Indeed, there is no mention of that in the docs. I will update it
> today to include this information. It is true that it's easy to break
> your teeth on this.
>
> Cheers,
> Ondra
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user