[GRASS-user] Grass, Python, mapcalc

Thomas Adams - NOAA Federal thomas.adams at noaa.gov
Fri Dec 28 08:24:53 PST 2012


I may be off base with what you want, but it seems like may you want to use
GRASS r.series, namely…

base=$(basename $input .grib2)
•
•
•
r.series input="`g.mlist rast pattern=${base} .* sep=,`"
output=${base".total method=sum --overwrite


So, if you have a bunch of *.grib2 files, for instance, that have been read
into GRASS that you want to sum; these GRASS raster maps all have the same
base filename root, $base.*

>From the GRASS manual:

*r.series* makes each output cell value a function of the values assigned
to the corresponding cells in the input raster map layers. Following
methods are available:
average: average value
count: count of non-NULL cells
median: median value
mode: most frequently occuring value
minimum: lowest value
maximum: highest value
range: range of values (max - min)
stddev: standard deviation
sum: sum of values
variance: statistical variance
diversity: number of different values
slope: linear regression slope
offset: linear regression offset
detcoeff: linear regression coefficient of determination
min_raster: raster map number with the minimum time-series value
max_raster: raster map number with the maximum time-series value


I hope this helps…

Tom



On Fri, Dec 28, 2012 at 7:36 AM, S. Koukoulas (lists) <
sotkouklistes at gmail.com> wrote:

> 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
>
> _______________________________________________
> grass-user mailing list
> grass-user at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user
>



-- 

Thomas E Adams

Development & Operations Hydrologist
National Weather Service
Ohio River Forecast Center
1901 South State Route 134
Wilmington, OH 45177

http://www.erh.noaa.gov/er/ohrfc/

EMAIL:	thomas.adams at noaa.gov
VOICE:	937-383-0528
FAX:	937-383-0033
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20121228/cc51fd15/attachment-0001.html>


More information about the grass-user mailing list