[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