[GRASS-dev] r.mapcalc bug on max and min?

Dear all,

I would like to rescale a raster map from a -1, 1 interval to 0, 255.

To make an example I've tried to normalize the elevation map:

{{{
$ r.mapcalc "el = elevation / max(elevation)" --o
100%
$ r.info el
+----------------------------------------------------------------------------+
| Map: el Date: Mon Oct 6 10:29:19 2014 |
| Mapset: user1 Login of Creator: pietro |
| Location: nc_basic_spm_grass7 |
| DataBase: /home/pietro/docdat/gis |
| Title: ( el ) |
| Timestamp: none |
|----------------------------------------------------------------------------|
| |
| Type of Map: raster Number of Categories: 0 |
| Data Type: FCELL |
| Rows: 1350 |
| Columns: 1500 |
| Total Cells: 2025000 |
| Projection: Lambert Conformal Conic |
| N: 228500 S: 215000 Res: 10 |
| E: 645000 W: 630000 Res: 10 |
| Range of data: min = 1 max = 1 |
| |
| Data Description: |
| generated by r.mapcalc |
| |
| Comments: |
| elevation / max(elevation) |
| |
+----------------------------------------------------------------------------+

}}}

As you can see the result is a map with all the values that are 1.

Instead If I execute:

{{{
$ r.mapcalc "el = elevation / 156.3299" --o
100%

$ r.info el
+----------------------------------------------------------------------------+
| Map: el Date: Mon Oct 6 10:31:17 2014 |
| Mapset: user1 Login of Creator: pietro |
| Location: nc_basic_spm_grass7 |
| DataBase: /home/pietro/docdat/gis |
| Title: ( el ) |
| Timestamp: none |
|----------------------------------------------------------------------------|
| |
| Type of Map: raster Number of Categories: 0 |
| Data Type: DCELL |
| Rows: 1350 |
| Columns: 1500 |
| Total Cells: 2025000 |
| Projection: Lambert Conformal Conic |
| N: 228500 S: 215000 Res: 10 |
| E: 645000 W: 630000 Res: 10 |
| Range of data: min = 0.355522472489405 max = 0.999999772928615 |
| |
| Data Description: |
| generated by r.mapcalc |
| |
| Comments: |
| elevation / 156.3299 |
| |
+----------------------------------------------------------------------------+
}}}

Now it si working... do you know why this r.mapcalc expression:
"el = elevation / max(elevation)" is wrong?

Now trying to normalize del "el" map:

{{{
$ RAST="el"
$ NEW_MAX=255.
$ NEW_MIN=0.
$ r.mapcalc "s$RAST = int(float($RAST - min($RAST)) / float(max($RAST)
- min($RAST)) * ($NEW_MAX - $NEW_MIN) + $NEW_MIN)" --o
100%
$ r.info sel
+----------------------------------------------------------------------------+
| Map: sel Date: Mon Oct 6 10:42:42 2014 |
| Mapset: user1 Login of Creator: pietro |
| Location: nc_basic_spm_grass7 |
| DataBase: /home/pietro/docdat/gis |
| Title: ( sel ) |
| Timestamp: none |
|----------------------------------------------------------------------------|
| |
| Type of Map: raster Number of Categories: 0 |
| Data Type: CELL |
| Rows: 1350 |
| Columns: 1500 |
| Total Cells: 2025000 |
| Projection: Lambert Conformal Conic |
| N: 228500 S: 215000 Res: 10 |
| E: 645000 W: 630000 Res: 10 |
| Range of data: min = NULL max = NULL |
| |
| Data Description: |
| generated by r.mapcalc |
| |
| Comments: |
| int(float(el - min(el)) / float(max(el) - min(el)) * (255 - 0) + 0) |
| |
+----------------------------------------------------------------------------+
}}}

Again everything is NULL, I'm doing something wrong or it is a bug?

Best regards

Pietro

On Mon, Oct 6, 2014 at 4:51 AM, Pietro <peter.zamb@gmail.com> wrote:

Dear all,

I would like to rescale a raster map from a -1, 1 interval to 0, 255.

To make an example I've tried to normalize the elevation map:

{{{
$ r.mapcalc "el = elevation / max(elevation)" --o
100%
$ r.info el

+----------------------------------------------------------------------------+
| Map: el Date: Mon Oct 6 10:29:19
2014 |
| Mapset: user1 Login of Creator: pietro
    |
| Location: nc_basic_spm_grass7
    |
| DataBase: /home/pietro/docdat/gis
    |
| Title: ( el )
    |
| Timestamp: none
    |

|----------------------------------------------------------------------------|
|
    |
| Type of Map: raster Number of Categories: 0
   |
| Data Type: FCELL
    |
| Rows: 1350
   |
| Columns: 1500
   |
| Total Cells: 2025000
    |
| Projection: Lambert Conformal Conic
   |
| N: 228500 S: 215000 Res: 10
   |
| E: 645000 W: 630000 Res: 10
   |
| Range of data: min = 1 max = 1
   |
|
    |
| Data Description:
    |
| generated by r.mapcalc
    |
|
    |
| Comments:
    |
| elevation / max(elevation)
    |
|
    |

+----------------------------------------------------------------------------+

}}}

As you can see the result is a map with all the values that are 1.

Instead If I execute:

{{{
$ r.mapcalc "el = elevation / 156.3299" --o
100%

$ r.info el

+----------------------------------------------------------------------------+
| Map: el Date: Mon Oct 6 10:31:17
2014 |
| Mapset: user1 Login of Creator: pietro
    |
| Location: nc_basic_spm_grass7
    |
| DataBase: /home/pietro/docdat/gis
    |
| Title: ( el )
    |
| Timestamp: none
    |

|----------------------------------------------------------------------------|
|
    |
| Type of Map: raster Number of Categories: 0
   |
| Data Type: DCELL
    |
| Rows: 1350
   |
| Columns: 1500
   |
| Total Cells: 2025000
    |
| Projection: Lambert Conformal Conic
   |
| N: 228500 S: 215000 Res: 10
   |
| E: 645000 W: 630000 Res: 10
   |
| Range of data: min = 0.355522472489405 max = 0.999999772928615
   |
|
    |
| Data Description:
    |
| generated by r.mapcalc
    |
|
    |
| Comments:
    |
| elevation / 156.3299
    |
|
    |

+----------------------------------------------------------------------------+
}}}

Now it si working... do you know why this r.mapcalc expression:
"el = elevation / max(elevation)" is wrong?

Now trying to normalize del "el" map:

{{{
$ RAST="el"
$ NEW_MAX=255.
$ NEW_MIN=0.
$ r.mapcalc "s$RAST = int(float($RAST - min($RAST)) / float(max($RAST)
- min($RAST)) * ($NEW_MAX - $NEW_MIN) + $NEW_MIN)" --o
100%
$ r.info sel

+----------------------------------------------------------------------------+
| Map: sel Date: Mon Oct 6 10:42:42
2014 |
| Mapset: user1 Login of Creator: pietro
    |
| Location: nc_basic_spm_grass7
    |
| DataBase: /home/pietro/docdat/gis
    |
| Title: ( sel )
   |
| Timestamp: none
    |

|----------------------------------------------------------------------------|
|
    |
| Type of Map: raster Number of Categories: 0
   |
| Data Type: CELL
   |
| Rows: 1350
   |
| Columns: 1500
   |
| Total Cells: 2025000
    |
| Projection: Lambert Conformal Conic
   |
| N: 228500 S: 215000 Res: 10
   |
| E: 645000 W: 630000 Res: 10
   |
| Range of data: min = NULL max = NULL
   |
|
    |
| Data Description:
    |
| generated by r.mapcalc
    |
|
    |
| Comments:
    |
| int(float(el - min(el)) / float(max(el) - min(el)) * (255 - 0) + 0)
   |
|
    |

+----------------------------------------------------------------------------+
}}}

Again everything is NULL, I'm doing something wrong or it is a bug?

Glynn could perhaps give you some more precise comment, but basically, max
doesn't operate on the entire map but only on the current cell. So for
example, r.macalc "map_max = max(elevation1, elevation2, elevation3)"
creates a map where each cell is the largest value of the corresponding
cells in the 3 elevation maps.

Hope that helps,
Anna

Best regards

Pietro
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

On Mon, Oct 6, 2014 at 2:25 PM, Anna Petrášová <kratochanna@gmail.com> wrote:

On Mon, Oct 6, 2014 at 4:51 AM, Pietro <peter.zamb@gmail.com> wrote:

Again everything is NULL, I'm doing something wrong or it is a bug?

Glynn could perhaps give you some more precise comment, but basically, max
doesn't operate on the entire map but only on the current cell. So for
example, r.macalc "map_max = max(elevation1, elevation2, elevation3)"
creates a map where each cell is the largest value of the corresponding
cells in the 3 elevation maps.

Ahh, I misunderstood the documentation I thought it was not pixel by
pixel but for the whole raster map...
now the results are explained! :slight_smile:

Perhaps we could add two new functions, like: rmin and rmax that stay
for range min and range max that give this information, do you think
could be useful?

On 06/10/14 15:07, Pietro wrote:

On Mon, Oct 6, 2014 at 2:25 PM, Anna Petrášová <kratochanna@gmail.com> wrote:

On Mon, Oct 6, 2014 at 4:51 AM, Pietro <peter.zamb@gmail.com> wrote:

Again everything is NULL, I'm doing something wrong or it is a bug?

Glynn could perhaps give you some more precise comment, but basically, max
doesn't operate on the entire map but only on the current cell. So for
example, r.macalc "map_max = max(elevation1, elevation2, elevation3)"
creates a map where each cell is the largest value of the corresponding
cells in the 3 elevation maps.

Ahh, I misunderstood the documentation I thought it was not pixel by
pixel but for the whole raster map...
now the results are explained! :slight_smile:

Perhaps we could add two new functions, like: rmin and rmax that stay
for range min and range max that give this information, do you think
could be useful?

You can use either r.info (entire map) or r.univar (current region) for that and then reuse the info in r.mapcalc.

Moritz

On Mon, Oct 6, 2014 at 9:07 AM, Pietro <peter.zamb@gmail.com> wrote:

On Mon, Oct 6, 2014 at 2:25 PM, Anna Petrášová <kratochanna@gmail.com>
wrote:
> On Mon, Oct 6, 2014 at 4:51 AM, Pietro <peter.zamb@gmail.com> wrote:
>> Again everything is NULL, I'm doing something wrong or it is a bug?
>
>
> Glynn could perhaps give you some more precise comment, but basically,
max
> doesn't operate on the entire map but only on the current cell. So for
> example, r.macalc "map_max = max(elevation1, elevation2, elevation3)"
> creates a map where each cell is the largest value of the corresponding
> cells in the 3 elevation maps.

Ahh, I misunderstood the documentation I thought it was not pixel by
pixel but for the whole raster map...
now the results are explained! :slight_smile:

Perhaps we could add two new functions, like: rmin and rmax that stay
for range min and range max that give this information, do you think
could be useful?

Yes, definitely, I was also missing this functionality. I assume it
shouldn't be too difficult since range is stored in the map anyway.

Anna

On Mon, Oct 6, 2014 at 9:39 AM, Anna Petrášová <kratochanna@gmail.com>
wrote:

On Mon, Oct 6, 2014 at 9:07 AM, Pietro <peter.zamb@gmail.com> wrote:

On Mon, Oct 6, 2014 at 2:25 PM, Anna Petrášová <kratochanna@gmail.com>
wrote:
> On Mon, Oct 6, 2014 at 4:51 AM, Pietro <peter.zamb@gmail.com> wrote:
>> Again everything is NULL, I'm doing something wrong or it is a bug?
>
>
> Glynn could perhaps give you some more precise comment, but basically,
max
> doesn't operate on the entire map but only on the current cell. So for
> example, r.macalc "map_max = max(elevation1, elevation2, elevation3)"
> creates a map where each cell is the largest value of the corresponding
> cells in the 3 elevation maps.

Ahh, I misunderstood the documentation I thought it was not pixel by
pixel but for the whole raster map...
now the results are explained! :slight_smile:

Perhaps we could add two new functions, like: rmin and rmax that stay
for range min and range max that give this information, do you think
could be useful?

Yes, definitely, I was also missing this functionality. I assume it
shouldn't be too difficult since range is stored in the map anyway.

This should strictly work on computational region.

Another function which I miss is coalesce which would use first non-NULL
value from the left but this is unrelated.

Vaclav

Anna

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

On Mon, Oct 6, 2014 at 4:05 PM, Vaclav Petras <wenzeslaus@gmail.com> wrote:

On Mon, Oct 6, 2014 at 9:39 AM, Anna Petrášová <kratochanna@gmail.com>

On Mon, Oct 6, 2014 at 9:07 AM, Pietro <peter.zamb@gmail.com> wrote:

Perhaps we could add two new functions, like: rmin and rmax that stay
for range min and range max that give this information, do you think
could be useful?

Yes, definitely, I was also missing this functionality.

...

This should strictly work on computational region.

Yes (r.univar contains the respective code).

Markus

On Mon, Oct 6, 2014 at 10:52 AM, Markus Neteler <neteler@osgeo.org> wrote:

On Mon, Oct 6, 2014 at 4:05 PM, Vaclav Petras <wenzeslaus@gmail.com>
wrote:
> On Mon, Oct 6, 2014 at 9:39 AM, Anna Petrášová <kratochanna@gmail.com>
>> On Mon, Oct 6, 2014 at 9:07 AM, Pietro <peter.zamb@gmail.com> wrote:
>>> Perhaps we could add two new functions, like: rmin and rmax that stay
>>> for range min and range max that give this information, do you think
>>> could be useful?
>>
>> Yes, definitely, I was also missing this functionality.
...
> This should strictly work on computational region.

Yes (r.univar contains the respective code).

This would mean moving this code to the library.

Markus

On 06/10/14 17:19, Vaclav Petras wrote:

On Mon, Oct 6, 2014 at 10:52 AM, Markus Neteler <neteler@osgeo.org
<mailto:neteler@osgeo.org>> wrote:

    On Mon, Oct 6, 2014 at 4:05 PM, Vaclav Petras <wenzeslaus@gmail.com
    <mailto:wenzeslaus@gmail.com>> wrote:
    > On Mon, Oct 6, 2014 at 9:39 AM, Anna Petrášová <kratochanna@gmail.com <mailto:kratochanna@gmail.com>>
    >> On Mon, Oct 6, 2014 at 9:07 AM, Pietro <peter.zamb@gmail.com <mailto:peter.zamb@gmail.com>> wrote:
    >>> Perhaps we could add two new functions, like: rmin and rmax that stay
    >>> for range min and range max that give this information, do you think
    >>> could be useful?
    >>
    >> Yes, definitely, I was also missing this functionality.
    ...
    > This should strictly work on computational region.

    Yes (r.univar contains the respective code).

This would mean moving this code to the library.

I'm not going to stand in the way of those who really want to implement this, but why is this necessary ?

eval(r.univar -g $map)
r.mapcalc new=old/$max

(and equivalent in other programming languages) works perfectly. Just thinking that there are other priorities and that we are going further from the KISS principle...

Moritz

On Mon, Oct 6, 2014 at 5:53 PM, Moritz Lennert
<mlennert@club.worldonline.be> wrote:

I'm not going to stand in the way of those who really want to implement
this, but why is this necessary ?

eval(r.univar -g $map)
r.mapcalc new=old/$max

(and equivalent in other programming languages) works perfectly. Just
thinking that there are other priorities and that we are going further from
the KISS principle...

I can agree with your point, but it is not the topic of this message,
I open this topic because it is to a C/python level.
So please don't focus on the dummy task, It is only an example to test
the existing libraries, and I#m not able to understand where the
problem.is.

Best regards

Pietro

Pietro wrote:

Perhaps we could add two new functions, like: rmin and rmax that stay
for range min and range max that give this information

It's not possible to implement these as functions within the current
structure of r.mapcalc.

r.mapcalc's functions take row buffers as inputs, and return (i.e.
populate) a row buffer as output.

There's no way to write a function which takes a map name as input
(r.mapcalc doesn't have any data types except numbers). A map name
automatically evaluates to the map's contents.

IOW, for a function call such as "sqrt(input)", the function is called
for each row, with the current row of data as its input. The function
has no information on where those values came from (e.g. whether
directly from a map or from the result of some other function).

It would be more realistic to implement it as a modifier, e.g.
min#input, similar to how colour-table lookups are performed. But that
creates issues of its own: what's the type? DCELL? Or the same as the
input? Is it the value for the map as a whole or do we need to invoke
r.univar (or embed r.univar into r.mapcalc or libraster)? Are we going
to add similar modifiers for all of the other useful statistics? For
arbitrary quantiles?

do you think could be useful?

Not really. Just use r.info (or r.univar if you want the min/max for
the current region) and substitute the result into the expression.

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

On Wed, Oct 8, 2014 at 12:38 PM, Glynn Clements
<glynn@gclements.plus.com> wrote:

do you think could be useful?

Not really. Just use r.info (or r.univar if you want the min/max for
the current region) and substitute the result into the expression.

Yes, I did in this way.
Thank you.

Pietro