[GRASS-dev] casting row buffer in memory

Trying to circumvent the r.mapcalc min/max limit, I found out another
problem when I'm trying to cast a buffer from one type to another.

The code is:

{{{
from __future__ import print_function
import numpy as np

from grass.pygrass.raster import RasterRow
from grass.pygrass.messages import get_msgr

def scale(mn, mx, newmn=0, newmx=1, dtype=np.float32):
    print("from (%r, %r) to (%r, %r): %r" % (mn, mx, newmn, newmx, dtype))
    def do(row):
        return dtype((row - mn) / (mx - mn) * (newmx - newmn) + newmn)
    return do

def rescale(input, output, new_min=0, new_max=255, mtype='CELL'):
    dtype = np.int32 if mtype == 'CELL' else np.float32 if mtype ==
'FCELL' else np.float64
    msgr = get_msgr()
    with RasterRow(input, mode='r') as inp:
        scaler = scale(inp.info.min, inp.info.max, new_min, new_max, dtype)
        with RasterRow(output, mode='w', mtype=mtype, overwrite=True) as out:
            binp = inp[0]
            rows = inp._rows
            for i in range(rows):
                msgr.percent(i, rows, 1)
                out.put_row(scaler(inp.get_row(i, binp)))

rescale('el', 'cel255', 0, 255, mtype='CELL')
rescale('el', 'fel255', 0, 255, mtype='FCELL')
rescale('el', 'del255', 0, 255, mtype='DCELL')
}}}

and these are the results

{{{
$ python2 cmd.py
from (0.3555224724894053, 0.9999997729286152) to (0, 255): <type 'numpy.int32'>
from (0.3555224724894053, 0.9999997729286152) to (0, 255): <type
'numpy.float32'>
from (0.3555224724894053, 0.9999997729286152) to (0, 255): <type
'numpy.float64'>

$ r.info cel255
[cut]
    Range of data: min = -2147483405 max = 2147176055

$ r.info fel255
[cut]
    Range of data: min = -3.402784e+38 max = 3.402818e+38

$ r.info del255
[cut]
     Range of data: min = 0 max = 255
}}}

So even if I'm explicitly casting the result in the scaler function
with dtype(...), the result is not properly casted, do you have an
idea on how could I solve this problem?

All the best

Pietro

On 06/10/14 14:54, Pietro wrote:

Trying to circumvent the r.mapcalc min/max limit,

Just out of curiosity: why don't you just use r.recode or r.rescale ?

Moritz

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

On 06/10/14 14:54, Pietro wrote:

Trying to circumvent the r.mapcalc min/max limit,

Just out of curiosity: why don't you just use r.recode or r.rescale ?

r.rescale was my first option, but it just reclassify in two
categories 1 and 255.

I didn't thought about r.recode...but I think that easier combining r.info.

On 06/10/14 16:14, Pietro wrote:

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

On 06/10/14 14:54, Pietro wrote:

Trying to circumvent the r.mapcalc min/max limit,

Just out of curiosity: why don't you just use r.recode or r.rescale ?

r.rescale was my first option, but it just reclassify in two
categories 1 and 255.

r.rescale assumes that input is integer, so if you feed it values between 0 and 1 then yes, but using the NC data as example and feeding it the elevation data directly:

r.rescale elevation out=elevation_0_255 to=0,255

gives the expected result.

I didn't thought about r.recode...but I think that easier combining r.info.

Why would you need r.info for that ?

Dividing by the max and then rescaling should give the same result as rescaling from the original values... All you are doing is dividing all values by a constant before rescaling.

Moritz

Il 06/ott/2014 17:04 “Moritz Lennert” <mlennert@club.worldonline.be> ha scritto:

On 06/10/14 16:14, Pietro wrote:

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

On 06/10/14 14:54, Pietro wrote:

Trying to circumvent the r.mapcalc min/max limit,

Just out of curiosity: why don’t you just use r.recode or r.rescale ?

r.rescale was my first option, but it just reclassify in two
categories 1 and 255.

r.rescale assumes that input is integer, so if you feed it values between 0 and 1 then yes, but using the NC data as example and feeding it the elevation data directly:

r.rescale elevation out=elevation_0_255 to=0,255

gives the expected result.

To make the example reproducible I’ve scaled the elevation map to have all the values between 0 and 1.

So as you said r.rescale is not an option.

I didn’t thought about r.recode…but I think that easier combining r.info.

Why would you need r.info for that ?

Sorry I was in harry to catch the train… I was too succinct… I mean as you suggested combine r.mapcalc and r.info to rescale…

Another option is to modify r.rescale to work with CELL,FCELL, and DCELL.
I don’t see any particular reason to limit the module to the cell case…

Dividing by the max and then rescaling should give the same
result as rescaling from the original values… All you are
doing is dividing all values by a constant before rescaling.

Yes, this was just an example using the north Carolina data set to make the problem reproducible…

But it is true I can also multiply my map to 10000 and then use r.rescale. I think should work.

On 06/10/14 17:31, Pietro wrote:

Il 06/ott/2014 17:04 "Moritz Lennert" <mlennert@club.worldonline.be
<mailto:mlennert@club.worldonline.be>> ha scritto:
>
> On 06/10/14 16:14, Pietro wrote:
>>
>> On Mon, Oct 6, 2014 at 3:17 PM, Moritz Lennert
>> <mlennert@club.worldonline.be <mailto:mlennert@club.worldonline.be>>
wrote:
>>>
>>> On 06/10/14 14:54, Pietro wrote:
>>>>
>>>> Trying to circumvent the r.mapcalc min/max limit,
>>>
>>> Just out of curiosity: why don't you just use r.recode or r.rescale ?
>>
>> r.rescale was my first option, but it just reclassify in two
>> categories 1 and 255.
>
> r.rescale assumes that input is integer, so if you feed it values
between 0 and 1 then yes, but using the NC data as example and feeding
it the elevation data directly:
>
> r.rescale elevation out=elevation_0_255 to=0,255
>
> gives the expected result.

To make the example reproducible I've scaled the elevation map to have
all the values between 0 and 1.

So as you said r.rescale is not an option.

>> I didn't thought about r.recode...but I think that easier combining
r.info <http://r.info>.
>
> Why would you need r.info <http://r.info> for that ?

Sorry I was in harry to catch the train... I was too succinct... I mean
as you suggested combine r.mapcalc and r.info <http://r.info> to rescale...

Another option is to modify r.rescale to work with CELL,FCELL, and DCELL.
I don't see any particular reason to limit the module to the cell case...

> Dividing by the max and then rescaling should give the same
> result as rescaling from the original values... All you are
> doing is dividing all values by a constant before rescaling.

Yes, this was just an example using the north Carolina data set to make
the problem reproducible...

But it is true I can also multiply my map to 10000 and then use
r.rescale. I think should work.

Or just use r.recode which is meant for FCELL and DCELL. A part from a small bug that made using -1.0:1.0:0:255 difficult (#2053), it works perfectly.

Moritz

Pietro wrote:

So even if I'm explicitly casting the result in the scaler function
with dtype(...), the result is not properly casted, do you have an
idea on how could I solve this problem?

This is a bug in pygrass:

  $ r.info -g foo | fgrep datatype
  datatype=CELL
  $ python
  Python 2.7.7 (default, Sep 26 2014, 00:18:15)
  [GCC 4.6.3] on linux2
  Type "help", "copyright", "credits" or "license" for more information.
  > from grass.pygrass.raster import RasterRow
  > out = RasterRow('foo', mode='w', mtype='DCELL', overwrite=True)
  > out.mtype
  u'CELL'

If a map with the chosen name already exists, its type is used in
place of the one specified.

Looking at the code for grass.pygrass.raster.RasterRow.open(), the
overall structure is rather flawed. It primarily considers whether the
map already exists, then secondarily considers the mode (r/w).

It should be the other way around.

Whether you're reading an existing map or creating a new map is
fundamental. Whether the map exists is secondary; whether it "ought"
to exist and what to do if it does/doesn't is entirely different
depending upon whether you're reading or writing.

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

On Wed, Oct 8, 2014 at 11:24 AM, Glynn Clements
<glynn@gclements.plus.com> wrote:

Pietro wrote:

So even if I'm explicitly casting the result in the scaler function
with dtype(...), the result is not properly casted, do you have an
idea on how could I solve this problem?

This is a bug in pygrass:

        $ r.info -g foo | fgrep datatype
        datatype=CELL
        $ python
        Python 2.7.7 (default, Sep 26 2014, 00:18:15)
        [GCC 4.6.3] on linux2
        Type "help", "copyright", "credits" or "license" for more information.
        > from grass.pygrass.raster import RasterRow
        > out = RasterRow('foo', mode='w', mtype='DCELL', overwrite=True)
        > out.mtype
        u'CELL'

If a map with the chosen name already exists, its type is used in
place of the one specified.

Thank you to find out this bug, I didn't noticed it.

Looking at the code for grass.pygrass.raster.RasterRow.open(), the
overall structure is rather flawed. It primarily considers whether the
map already exists, then secondarily considers the mode (r/w).

It should be the other way around.

Right. Thank you for your code review!
I will try to fix this on the next days.

All the best.

Pietro

Pietro wrote:

> Looking at the code for grass.pygrass.raster.RasterRow.open(), the
> overall structure is rather flawed. It primarily considers whether the
> map already exists, then secondarily considers the mode (r/w).
>
> It should be the other way around.

Right. Thank you for your code review!
I will try to fix this on the next days.

Having thought about it some more, I'd suggest giving serious
consideration to splitting the RasterRow class into separate classes
for input and output.

The two are sufficiently different that I'm not entirely sure that a
common base class wouldn't be empty. Most of the libraster functions
which accept an open map descriptor ("fd") either only work on input
(old) maps or only work on output (new) maps.

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