[GRASS-dev] grass, python, Pygrass

[Ivan Marchesini]

Hi Pietro
many thanks for your answer


> Some note on the code:
> 
> 1) maybe to this kind of things you should consider to use the
> r.mapcalc instead of pygrass;

Yes you are true but later, in the code, I need to do some more
complicated stuf, like increasing or decreasing the value of a cell
based on 
1) the value of direction and inclination stored in the same position
into two other maps and  and 
2) the three cells that are in the opposite side respect to the
direction. 

As an example: Suppose I'm considering the cell[i,j].
>From the map of directions (drainage) I know that the direction[i,j] is
towards cell[1,j+1] (towards east). From the inclination map I know that
the inclination[i,j] is 45°.
In my algorithm I would like to estimate a new value of the elevation
for the cell[i,j] based on the contribution of the elevation values
stored on the cells [i-1,j-1], [1,j-1],[i+1,j-1].


> 2) Why are you using RasterNumpy?
>     - RasterNumpy is quite slow because need to export your map on the
> hard disk, load in memory using numpy.memmap... therefore, if
> possible, try to avoid this class...
>     - Avoid to use cycle in python, they are really slow and are
> working only with small region.

ok but it is possible to work with indexed cells [i,j] as I described
before.

> Following the first point the cycle will be:
> 
> {{{
> from grass import script as grass
> 
> grass.run_command("r.mapcalc", expression="{outmap} = if({a}==1, 1,
> 3)".format(a = "tmp_ocs", outmap="visited"), overwrite = True)
> }}}

ok.. thank you. I was stupid not to think at that for this simple
cycle... It happened because I was already thinking about the second
part.. the one of direction and inclination  

> 
> About the second point:
> 
> {{{
> from grass.pygrass.raster import RasterRow
> 
> with RasterRow("tmp_ocs") as ocs:
>     new = RasterRow("visited")
>     new.open("w", "CELL")
>     for row in ocs:
>         rbool = (row == 1)
>         row[rbool] = 1
>         row[~rbool] = 3
>         new.write(row)
>     new.close()
> }}}

yes.. good.. this is also another way. thank you


> In this piece of code you are assign a new name to the raster, but the
> map already exist...
> Therefore as Luca said, should be sufficient to add:
> 
> {{{
> visited.overwrite = True
> visited.name = 'visitati'
> }}}
> 
> let us know.

As I said to luca probably the problem was simply due to 4 spaces in
place of a tab

Many many thanks and please let me know if you think that for the
problem of direction and inclination the use of numpy is correct

Ivan