[GRASS-dev] grass, python, Pygrass

[Pietro]

Ciao Ivan,

On Wed, Jun 26, 2013 at 5:40 AM, Ivan Marchesini
<ivan.marchesini at gmail.com> wrote:
> I have this small code:
>
>     grass.run_command("v.to.rast", input=ocs, type="line", output="tmp_ocs", use="val", overwrite=True)
>     grass.mapcalc("$outmap = $a", a = dtm, outmap="demcut", overwrite = true)
>     grass.run_command("g.remove", rast="visitati")
>
>     np_ocs = raster.RasterNumpy("tmp_ocs")
>     np_ocs.open()
>     visited=np_ocs*0
>     visited.open()
>
>     for i in range(1,new.shape[0]-1):
>         for j in range(1,new.shape[1]-1):
>             if np_ocs[i,j] == 1:
>                 visited[i,j] = 1
>             else:
>                 visited[i,j] = 3
>
>     visited.name='visitati'
>     visited.close()
>     sys.exit("stop")

Some note on the code:

1) maybe to this kind of things you should consider to use the
r.mapcalc instead of pygrass;
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.

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)
}}}

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()
}}}

In this way only the cycle on the rows is done in python, all the
other cycle are done in C trough numpy.

But to do something so simple I suggest to use r.mapcalc.


> When I run it I obtain a lot of:
>
> WARNING: rast=visitati,visitati: files could be the same, no rename
>          possible
>
> I cannot understand where the warnings arise... and why..

I think the warning raise, here:

{{{
visited.name='visitati'
}}}

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.

Pietro