[GRASS-dev] problem using pygrass

GRASS GIS GRASS GIS developer
1

Hello,

I have a couple of problems when using pygrass for analyzing raster
maps stored in GRASS with the numpy library:

- First: I cannot close the maps, since I always get the following error:

import numpy
import grass.pygrass as pygrass
demx = pygrass.raster.RasterNumpy('dem')
demx.open('r')
dem2 = demx.flatten()
len(dem2)

43860

demx.close()

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/grass-7.0.y14/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 654, in close
    self._write()
  File "/usr/local/grass-7.0.y14/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 597, in _write
    self.tofile(self.filename)
ValueError: 43860 requested and 0 written

- Second: when I change the region (with g.region) and read again some
maps, it uses the same extent as before and I always get a "core dump
error". I print the region in the script and it has in deed changed
but I always get maps from the same size as before...

Could anybody help?

Thank you very much and best regards,

Javier

2

Hi Javier,

On Thu, Feb 27, 2014 at 4:45 PM, Javier Martínez-López
<javi.martinez.lopez@gmail.com> wrote:

I have a couple of problems when using pygrass for analyzing raster
maps stored in GRASS with the numpy library:

The RasterNumpy class was quite buggy, I did some changes last week
that should fix some of them, which version are you using?

- First: I cannot close the maps, since I always get the following error:

import numpy
import grass.pygrass as pygrass
demx = pygrass.raster.RasterNumpy('dem')
demx.open('r')
dem2 = demx.flatten()
len(dem2)

43860

demx.close()

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/grass-7.0.y14/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 654, in close
    self._write()
  File "/usr/local/grass-7.0.y14/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 597, in _write
    self.tofile(self.filename)
ValueError: 43860 requested and 0 written

yes, this bug should be fixed in (r59127), I've tried this code using
the North Carolina mapset:

from grass.pygrass.raster import RasterNumpy
elev = RasterNumpy('elevation')
elev.open()
flt = elev.flatten()
len(flt)
elev.close()

and it works, are you in Linux?

- Second: when I change the region (with g.region) and read again some
maps, it uses the same extent as before and I always get a "core dump
error". I print the region in the script and it has in deed changed
but I always get maps from the same size as before...

in numpy the raster map is read and load into the memory when you open
the map using the current region, so changing the region when your map
is opened will not affect the raster map. Using the other RasterClass
will easily raise some errors.
So it is not a worth idea change a region when one or more raster maps
are opened, the right way to go is: set the region, do something with
the raster map, close then you are free to change the region again.

I would like to introduce the with statement also for the Region
object, to do something like:

{{{
## change resolution
with Region(nsres=100., ewres=100.) as reg:
    # do something with the changed region
    #

## back to the original region
#
}}}

3

Dear Pietro,

thank you very much for your help! I will upgrade the GRASS version
and try again, and will close the maps before changing the region, if
possible. I am using CentOS linux. This is a great tool!

Cheers,

Javier

On Thu, Feb 27, 2014 at 10:06 PM, Pietro <peter.zamb@gmail.com> wrote:

Hi Javier,

On Thu, Feb 27, 2014 at 4:45 PM, Javier Martínez-López
<javi.martinez.lopez@gmail.com> wrote:

I have a couple of problems when using pygrass for analyzing raster
maps stored in GRASS with the numpy library:

The RasterNumpy class was quite buggy, I did some changes last week
that should fix some of them, which version are you using?

- First: I cannot close the maps, since I always get the following error:

import numpy
import grass.pygrass as pygrass
demx = pygrass.raster.RasterNumpy('dem')
demx.open('r')
dem2 = demx.flatten()
len(dem2)

43860

demx.close()

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/grass-7.0.y14/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 654, in close
    self._write()
  File "/usr/local/grass-7.0.y14/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 597, in _write
    self.tofile(self.filename)
ValueError: 43860 requested and 0 written

yes, this bug should be fixed in (r59127), I've tried this code using
the North Carolina mapset:

from grass.pygrass.raster import RasterNumpy
elev = RasterNumpy('elevation')
elev.open()
flt = elev.flatten()
len(flt)
elev.close()

and it works, are you in Linux?

- Second: when I change the region (with g.region) and read again some
maps, it uses the same extent as before and I always get a "core dump
error". I print the region in the script and it has in deed changed
but I always get maps from the same size as before...

in numpy the raster map is read and load into the memory when you open
the map using the current region, so changing the region when your map
is opened will not affect the raster map. Using the other RasterClass
will easily raise some errors.
So it is not a worth idea change a region when one or more raster maps
are opened, the right way to go is: set the region, do something with
the raster map, close then you are free to change the region again.

I would like to introduce the with statement also for the Region
object, to do something like:

{{{
## change resolution
with Region(nsres=100., ewres=100.) as reg:
    # do something with the changed region
    #

## back to the original region
#
}}}

4

On 27 February 2014 22:06, Pietro <peter.zamb@gmail.com> wrote:

Hi Javier,

Hi,

yes, this bug should be fixed in (r59127), I've tried this code using
the North Carolina mapset:

from grass.pygrass.raster import RasterNumpy
elev = RasterNumpy('elevation')
elev.open()
flt = elev.flatten()
len(flt)
elev.close()

and it works, are you in Linux?

I confirm this.
I would also add that it if you want to save the map you have to use:

flt = pygrass.raster.RasterNumpy('elevation_flatten')
flt.open('w',mtype='FCELL')
flt = elev.flatten()
elev.close()
flt.close()

otherwise you obtain an error (segmentation fault)

--
ciao
Luca

http://gis.cri.fmach.it/delucchi/
www.lucadelu.org

5

Hello again,

I keep testing it and I get the following errors when trying to save
the results into a new raster map:

out = pygrass.raster.RasterNumpy('eco_pa22') # eco_pa22 is a raster
map which I created just to fill it in with new values

out.open('w',mtype='FCELL')

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: open() got multiple values for keyword argument 'mtype'

out.open('w')

out.max()
1
out.min()
0

out = np.where(out >= 0,(pmh),(pmh)) # filling the map with the values from pmh

out.max()

RasterNumpy(0.9995180622329483)

out.min()

RasterNumpy(0.0)

out.close()

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 655, in close
    np.memmap._close(self)
AttributeError: type object 'memmap' has no attribute '_close'

out

RasterNumpy([[ 0., 0., 0., ..., 0., 0., 0.],
       [ 0., 0., 0., ..., 0., 0., 0.],
       [ 0., 0., 0., ..., 0., 0., 0.],
       ...,
       [ 0., 0., 0., ..., 0., 0., 0.],
       [ 0., 0., 0., ..., 0., 0., 0.],
       [ 0., 0., 0., ..., 0., 0., 0.]])

Am I doing something wrong? Is there an easy way to get the raster map
out as a tiff file instead to avoid this problem?

Thank you very much!

Cheers,

Javier

On Sat, Mar 1, 2014 at 12:11 AM, Luca Delucchi <lucadeluge@gmail.com> wrote:

On 27 February 2014 22:06, Pietro <peter.zamb@gmail.com> wrote:

Hi Javier,

Hi,

yes, this bug should be fixed in (r59127), I've tried this code using
the North Carolina mapset:

from grass.pygrass.raster import RasterNumpy
elev = RasterNumpy('elevation')
elev.open()
flt = elev.flatten()
len(flt)
elev.close()

and it works, are you in Linux?

I confirm this.
I would also add that it if you want to save the map you have to use:

flt = pygrass.raster.RasterNumpy('elevation_flatten')
flt.open('w',mtype='FCELL')
flt = elev.flatten()
elev.close()
flt.close()

otherwise you obtain an error (segmentation fault)

--
ciao
Luca

http://gis.cri.fmach.it/delucchi/
www.lucadelu.org

6

Hi Javier,

the RasterNumpy class is different from other Raster classes in
pygrass, because inherit from the numpy.memmap array.

On Tue, Mar 4, 2014 at 2:38 PM, Javier Martínez-López
<javi.martinez.lopez@gmail.com> wrote:

out = pygrass.raster.RasterNumpy('eco_pa22')

In the RasterNumpy class you need to specify the raster type with
mtype when you instance the object.
and the open method should be call only after you close the map.

so a working session in the python interpreter could be:

{{{
In [1]: from grass.pygrass.raster import RasterNumpy

In [2]: new = RasterNumpy('newraster', mtype='FCELL')

In [3]: import numpy as np

In [4]: new[:] = np.random.randn(*new.shape)

In [5]: new.min()
Out[5]: RasterNumpy(-5.076531410217285, dtype=float32)

In [6]: new.max()
Out[6]: RasterNumpy(4.842053413391113, dtype=float32)

In [7]: new[:3, :4]
Out[7]:
RasterNumpy([[-0.74486554, -1.80636668, -0.62290537, 1.10075259],
       [ 0.65961069, -0.71186149, 1.31982756, 0.70869118],
       [ 2.11516094, -1.15944469, 0.71266735, 0.52502483]], dtype=float32)

In [8]: new.close()

}}}

Am I doing something wrong?

You should not use the open method with the RasterNumpy class when you
instantiate the object, and be careful that you have to set the name
of your raster map before to close the map or when you close the map,
see this example:

{{{
In [9]: new.open() # re-open the closed map

In [10]: new[:3, :4]
Out[10]:
RasterNumpy([[-0.74486554, -1.80636668, -0.62290537, 1.10075259],
       [ 0.65961069, -0.71186149, 1.31982756, 0.70869118],
       [ 2.11516094, -1.15944469, 0.71266735, 0.52502483]], dtype=float32)

In [11]: norm = (new - new.min())/(new.max() - new.min()) # do something

In [12]: norm.min()
Out[12]: RasterNumpy(0.0, dtype=float32)

In [13]: norm.max()
Out[13]: RasterNumpy(1.0, dtype=float32)

In [14]: norm.name # the name is not set

In [15]: norm.name = 'norm' # set the name

In [16]: norm.close() # then close

In [17]: norm.close('norm') # give the name and close the raster,
only from: r59192 [0]

}}}

Is there an easy way to get the raster map
out as a tiff file instead to avoid this problem?

Why not using r.out.tiff [1]?

Thank you for testing.

Pietro

[0] http://trac.osgeo.org/grass/changeset/59192
[1] http://grass.osgeo.org/grass70/manuals/r.out.tiff.html

7

On Tue, Mar 4, 2014 at 10:25 AM, Pietro <peter.zamb@gmail.com> wrote:

> Is there an easy way to get the raster map
> out as a tiff file instead to avoid this problem?

Why not using r.out.tiff [1]?

Sometimes (don't ask me when) is better to use `r.out.gdal` [0] with
`format=GTiff`.

[0] http://grass.osgeo.org/grass70/manuals/r.out.gdal.html

8

Thank you Pietro! I will try again! However, there seems to be a bug
(at least in version 59147) since I am getting NaN values reading some
maps with RasterNumpy, which are perfectly fine when reading them
using grass70 version 57982... They are FCELL! Did you also notice
that? Cheers, Javier

On Tue, Mar 4, 2014 at 4:25 PM, Pietro <peter.zamb@gmail.com> wrote:

Hi Javier,

the RasterNumpy class is different from other Raster classes in
pygrass, because inherit from the numpy.memmap array.

On Tue, Mar 4, 2014 at 2:38 PM, Javier Martínez-López
<javi.martinez.lopez@gmail.com> wrote:

out = pygrass.raster.RasterNumpy('eco_pa22')

In the RasterNumpy class you need to specify the raster type with
mtype when you instance the object.
and the open method should be call only after you close the map.

so a working session in the python interpreter could be:

{{{
In [1]: from grass.pygrass.raster import RasterNumpy

In [2]: new = RasterNumpy('newraster', mtype='FCELL')

In [3]: import numpy as np

In [4]: new[:] = np.random.randn(*new.shape)

In [5]: new.min()
Out[5]: RasterNumpy(-5.076531410217285, dtype=float32)

In [6]: new.max()
Out[6]: RasterNumpy(4.842053413391113, dtype=float32)

In [7]: new[:3, :4]
Out[7]:
RasterNumpy([[-0.74486554, -1.80636668, -0.62290537, 1.10075259],
       [ 0.65961069, -0.71186149, 1.31982756, 0.70869118],
       [ 2.11516094, -1.15944469, 0.71266735, 0.52502483]], dtype=float32)

In [8]: new.close()

}}}

Am I doing something wrong?

You should not use the open method with the RasterNumpy class when you
instantiate the object, and be careful that you have to set the name
of your raster map before to close the map or when you close the map,
see this example:

{{{
In [9]: new.open() # re-open the closed map

In [10]: new[:3, :4]
Out[10]:
RasterNumpy([[-0.74486554, -1.80636668, -0.62290537, 1.10075259],
       [ 0.65961069, -0.71186149, 1.31982756, 0.70869118],
       [ 2.11516094, -1.15944469, 0.71266735, 0.52502483]], dtype=float32)

In [11]: norm = (new - new.min())/(new.max() - new.min()) # do something

In [12]: norm.min()
Out[12]: RasterNumpy(0.0, dtype=float32)

In [13]: norm.max()
Out[13]: RasterNumpy(1.0, dtype=float32)

In [14]: norm.name # the name is not set

In [15]: norm.name = 'norm' # set the name

In [16]: norm.close() # then close

In [17]: norm.close('norm') # give the name and close the raster,
only from: r59192 [0]

}}}

Is there an easy way to get the raster map
out as a tiff file instead to avoid this problem?

Why not using r.out.tiff [1]?

Thank you for testing.

Pietro

[0] http://trac.osgeo.org/grass/changeset/59192
[1] http://grass.osgeo.org/grass70/manuals/r.out.tiff.html

9

Hi Javier,

On Tue, Mar 4, 2014 at 5:26 PM, Javier Martínez-López
<javi.martinez.lopez@gmail.com> wrote:

Thank you Pietro! I will try again! However, there seems to be a bug
(at least in version 59147) since I am getting NaN values reading some
maps with RasterNumpy, which are perfectly fine when reading them
using grass70 version 57982... They are FCELL! Did you also notice
that?

Please check that the region is set correctly... all the cells are
Nan? check with:

np.nanmax(raster)

Are you able to reproduce the problem with the North Carolina dataset?
If yes please send me the procedure... and I will try to fix it.

All the best

Pietro

10

Dear Pietro,

not all values are nan (consistently 5567600 out of 9476388 cells),
while importing other CELL maps (same region) works fine.

When I try this:

np.nanmax(epr22)

I get the follwing error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/majavie/anaconda/lib/python2.7/site-packages/numpy/lib/nanfunctions.py",
line 322, in nanmax
    a, mask = _replace_nan(a, -np.inf)
  File "/home/majavie/anaconda/lib/python2.7/site-packages/numpy/lib/nanfunctions.py",
line 63, in _replace_nan
    np.copyto(a, val, where=mask)
TypeError: Cannot cast array data from dtype('int32') to dtype('bool')
according to the rule 'safe'

I will check it with the NC dataset and let you know... Thank you and
cheers, Javier

On Tue, Mar 4, 2014 at 7:52 PM, Pietro <peter.zamb@gmail.com> wrote:

Hi Javier,

On Tue, Mar 4, 2014 at 5:26 PM, Javier Martínez-López
<javi.martinez.lopez@gmail.com> wrote:

Thank you Pietro! I will try again! However, there seems to be a bug
(at least in version 59147) since I am getting NaN values reading some
maps with RasterNumpy, which are perfectly fine when reading them
using grass70 version 57982... They are FCELL! Did you also notice
that?

Please check that the region is set correctly... all the cells are
Nan? check with:

np.nanmax(raster)

Are you able to reproduce the problem with the North Carolina dataset?
If yes please send me the procedure... and I will try to fix it.

All the best

Pietro

11

... something is wrong with my installation... I cannot even open any
map from the NC dataset using build 59147, getting the following error
after interrupting the idle process [1], while doing it with version
57982 works fine... any idea or advice? Here is the built
configuration that was used for the 59147 version:

./configure --prefix=/usr/local/grass-7.0.2 --with-libs=/usr/lib64
--with-cxx --without-ffmpeg --with-gdal=/usr/bin/gdal-config
--with-geos=/usr/bin/geos-config --without-odbc --with-sqlite
--without-postgres --without-mysql --with-nls --with-python
--with-cairo --with-wxwidgets=/usr/bin/wx-config --without-fftw
--with-freetype --with-freetype-includes=/usr/include/freetype2
--enable-largefile --with-pthread --with-lapack

... and the one from version 57982:

./configure --prefix=/usr/local --with-libs=/usr/lib64 --with-cxx
--without-ffmpeg --with-gdal=/usr/bin/gdal-config
--with-geos=/usr/bin/geos-config --without-odbc --with-sqlite
--without-postgres --without-mysql --with-nls --with-python
--with-cairo --with-wxwidgets=/usr/bin/wx-config --without-fftw
--with-freetype --with-freetype-includes=/usr/include/freetype2
--enable-largefile --with-pthread

the only difference is "--with-lapack" in the upper one...is there
something wrong with it?

Thank you and cheers,

Javier

[1]

import os
import numpy as np
import grass.pygrass as pygrass
lu = pygrass.raster.RasterNumpy('landuse')
lu.shape

(474, 526)

lu.open('r')

^CTraceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/grass-7.0.2/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 634, in open
    self._read()
  File "/usr/local/grass-7.0.2/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 597, in _read
    self[i] = rst.get_row(i, buff)
  File "/usr/local/grass-7.0.2/grass-7.0.svn/etc/python/grass/pygrass/raster/__init__.py",
line 531, in __array_finalize__
    self.filename = grasscore.tempfile()
  File "/usr/local/grass-7.0.2/grass-7.0.svn/etc/python/grass/script/core.py",
line 680, in tempfile
    return read_command("g.tempfile", flags=flags, pid=os.getpid()).strip()
  File "/usr/local/grass-7.0.2/grass-7.0.svn/etc/python/grass/script/core.py",
line 412, in read_command
    return ps.communicate()[0]
  File "/home/majavie/anaconda/lib/python2.7/subprocess.py", line 790,
in communicate
    stdout = _eintr_retry_call(self.stdout.read)
  File "/home/majavie/anaconda/lib/python2.7/subprocess.py", line 476,
in _eintr_retry_call
    return func(*args)
KeyboardInterrupt

On Tue, Mar 4, 2014 at 7:52 PM, Pietro <peter.zamb@gmail.com> wrote:

Hi Javier,

On Tue, Mar 4, 2014 at 5:26 PM, Javier Martínez-López
<javi.martinez.lopez@gmail.com> wrote:

Thank you Pietro! I will try again! However, there seems to be a bug
(at least in version 59147) since I am getting NaN values reading some
maps with RasterNumpy, which are perfectly fine when reading them
using grass70 version 57982... They are FCELL! Did you also notice
that?

Please check that the region is set correctly... all the cells are
Nan? check with:

np.nanmax(raster)

Are you able to reproduce the problem with the North Carolina dataset?
If yes please send me the procedure... and I will try to fix it.

All the best

Pietro

12

Hi Javier,

On Wed, Mar 5, 2014 at 10:33 AM, Javier Martínez-López
<javi.martinez.lopez@gmail.com> wrote:

... something is wrong with my installation... I cannot even open any
map from the NC dataset using build 59147, getting the following error
after interrupting the idle process [1], while doing it with version
57982 works fine... any idea or advice? Here is the built
configuration that was used for the 59147 version:

[cut]

the only difference is "--with-lapack" in the upper one...is there
something wrong with it?

I don't think so, I'm using this flag too.

Try to build everything from scranch, with:

make distclean
./configure blabla
make -j4

I've run your code using the current svn version (r59193) and it is
working on my laptop.
The only thing is that the region in my case is set to the elevation
raster map, so the shape is different:

In [1]: import numpy as np

In [2]: import grass.pygrass as pygrass

In [3]: lu = pygrass.raster.RasterNumpy('landuse')

In [4]: lu.shape
Out[6]: (1350, 1500)

In [5]: lu.open()

In [6]: lu[:4, :3]
Out[6]:
RasterNumpy([[5, 5, 5],
       [5, 5, 5],
       [5, 5, 5],
       [5, 3, 3]], dtype=int32)

In [7]: lu.close()

13

Dear Pietro,

the issue is apparently not related to the GRASS build but to the fact
that when there is a mask on the region, the areas outside it are
converted to "nan" in the numpy matrix when the map is FCELL, and to
the value "-2147483648" when the map is CELL... at least for me! I
need to apply a mask because I intend to set up an iterative process
in different regions, but I guess it can be just solved by filtering
the matrix before the analysis. I am uploading the code to a github
repo [1], just in case anyone is interested or it could serve as an
example.

Thank you and cheers,

Javier

[1] https://github.com/javimarlop/eHab_grasspy

On Wed, Mar 5, 2014 at 12:15 PM, Pietro <peter.zamb@gmail.com> wrote:

Hi Javier,

On Wed, Mar 5, 2014 at 10:33 AM, Javier Martínez-López
<javi.martinez.lopez@gmail.com> wrote:

... something is wrong with my installation... I cannot even open any
map from the NC dataset using build 59147, getting the following error
after interrupting the idle process [1], while doing it with version
57982 works fine... any idea or advice? Here is the built
configuration that was used for the 59147 version:

[cut]

the only difference is "--with-lapack" in the upper one...is there
something wrong with it?

I don't think so, I'm using this flag too.

Try to build everything from scranch, with:

make distclean
./configure blabla
make -j4

I've run your code using the current svn version (r59193) and it is
working on my laptop.
The only thing is that the region in my case is set to the elevation
raster map, so the shape is different:

In [1]: import numpy as np

In [2]: import grass.pygrass as pygrass

In [3]: lu = pygrass.raster.RasterNumpy('landuse')

In [4]: lu.shape
Out[6]: (1350, 1500)

In [5]: lu.open()

In [6]: lu[:4, :3]
Out[6]:
RasterNumpy([[5, 5, 5],
       [5, 5, 5],
       [5, 5, 5],
       [5, 3, 3]], dtype=int32)

In [7]: lu.close()