[GRASS-user] Grass, Python, mapcalc

Daniel Lee lee at isi-solutions.org
Fri Dec 28 03:48:35 PST 2012


Hmm, okay. That post is pretty old though and there have been a lot of
improvements to GRASS since. I don't think it applies any more. Check this
out:

GRASS 6.4.3svn (EPSG25832_UTM32N):~ > r.mapcalc test=1
 100%
GRASS 6.4.3svn (EPSG25832_UTM32N):~ > r.mapcalc test=test+3
 100%
GRASS 6.4.3svn (EPSG25832_UTM32N):~ > r.mapcalc "test = test + 3"
 100%

This worked for me. Make sure that you don't forget to use quote marks if
you have spaces in the middle of your expression. I think that's your
problem.

Daniel


2012/12/28 Paul Meems <bontepaarden at gmail.com>

> 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.htmlsaying:
>
> >* 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
>
> Owner of MapWindow.nl - Support for
> Dutch speaking users.
> www.mapwindow.nl
>
> *
> *
>
>
> 2012/12/28 Daniel Lee <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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20121228/aaf44911/attachment.html>


More information about the grass-user mailing list