[GRASS-user] Grass, Python, mapcalc

[S. Koukoulas (lists)]

To my experience with r.mapcalc you cannot do what yo do with
programming languages... e.g. i=i+1 will not work directly. Instead you
could try a more indirect way, something like (here is an example with
bash shell script):

#g.mremove rast=out* (to remove previous stuff, if needed)
r.mapcalc sun=global1
for a in ‘seq 2 10‘; do
r.mapcalc out = sun + global$a
g.rename --o rast=out,sun
done

basically the main idea is to use an intermediate map and then before
the loop ends, pass it (with g.rename) to the map variable you want.
There might be a more elegant way with python, but hopefully this shell
script example might help you.
HTH,
sotiris


On 12/28/2012 01:34 PM, Paul Meems wrote:
> Thanks again Daniel for your help.
>
> When I do the calculation in the GUI instead of in the Python shell I
> get these results:
> r.mapcalc sun = global1 at temp+ global2 at temp
> (Fri Dec 28 12:26:10 2012) Command finished (1 sec)
> r.mapcalc sun = sun at temp+ global3 at temp
> ERROR: Unable to close raster map
>
> When I Google on the error message I find this post from 2010
> http://lists.osgeo.org/pipermail/grass-user/2010-October/058463.html
> saying:
> >/ Why am I obtaining this error? and How can I avoid it?
> /
> Because you're using the same map as both input and output, which
> won't work. Either use a temporary name for the intermediate map, or
> replace the second grass.mapcalc call with r.null.
> I'll have a look in using a temp file. This will make my script more
> complicated ;)
>
> Thanks,
>
> Paul
>
> *Paul Meems *
> Release manager, configuration manager
> and forum moderator of MapWindow GIS.
> www.mapwindow.org <http://www.mapwindow.org/>
>
> Owner of MapWindow.nl - Support for
> Dutch speaking users.
> www.mapwindow.nl <http://www.mapwindow.nl/>
>
> *
> *
>
>
>
> 2012/12/28 Daniel Lee <lee at isi-solutions.org
> <mailto:lee at isi-solutions.org>>
>
>     Hi there,
>
>         I'm now trying to use r.sun and r.mapcalc with Python to
>         create my annual solar radiation map.
>         I'm using the Python shell inside GRASS.
>
>         For testing purposes I start with 9 days.
>         The r.sun part seems to be working:
>         for x in range(1, 10, 1):
>         print "Working on day %d" % (x)
>         glob_rad = 'global' + str(x)
>         grass.run_command('r.sun', flags = 's', elevin =
>         'w001001 at temp', aspin = 'aspect at temp', slopein = 'slope at temp',
>         glob_rad = glob_rad, day = x)
>         I do have a question about parsing the --overwrite flag. How
>         to do that. Adding flags='s, --overwrite' or flags='s,
>         -overwrite' gives a compile error.
>
>
>     I believe you overwrite by setting:
>     flags="s, o"
>
>
>         Now I try to combine the results of r.sun into 1 raster.
>         This is working:
>         for x in range(1, 10, 1):
>         if x == 2:
>         exp = 'sun = global1 at temp + global2 at temp'
>         print exp
>         grass.mapcalc(exp, quiet=False, verbose=False, overwrite=True)
>
>         This is also working:
>         for x in range(1, 10, 1):
>         if x == 2:
>         exp = 'sun = global1 at temp + global2 at temp'
>         print exp
>         grass.mapcalc(exp, quiet=False, verbose=False, overwrite=True)
>         elif x > 2:
>         exp = 'sun = sun at temp + global' + str(x) + '@temp'
>         print exp
>         and produces this output:
>         sun = sun at temp + global3 at temp
>         sun = sun at temp + global4 at temp
>         sun = sun at temp + global5 at temp
>         sun = sun at temp + global6 at temp
>         sun = sun at temp + global7 at temp
>         sun = sun at temp + global8 at temp
>         sun = sun at temp + global9 at temp
>         But when I execute this using
>         grass.mapcalc(exp, quiet=False, verbose=False, overwrite=True)
>         Grass crashes completely.
>
>         Has this something to do with having the output file is the
>         input file as well?
>         If this is not allowed how can I combine all output files from
>         r.sun into 1 raster?
>
>
>     GRASS doesn't care what you're combining, when you use the map
>     calculator it just sees numbers so that shouldn't be a problem.
>     I don't know why the map calculator to overwrite the input map
>     would be a problem either. I'm not familiar with your syntax
>     though. What happens when you just enter the expression on the
>     "normal" command line?
>
>     r.mapcalc "sun = sun + global3"
>
>     That should work.
>
>     Best,
>     Daniel
>
>
>
>
> _______________________________________________
> grass-user mailing list
> grass-user at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user