[GRASS-user] Grass, Python, mapcalc

[Paul Meems]

Hi Sotiris,

Thanks for your help.
I will have a look at the r.series.

Meanwhile I've discovered the graphical modeler. It is not working 100% in
all cases but it does generates a Python script with all the necessary
parameters.

This is my script so far, which seems to be working:
for x in range(1, 10, 1):
    print "Working on day %d" % (x)
    # Create radiation for each day:
    glob_rad = 'global' + str(x)
    grass.run_command('r.sun', flags = 's', overwrite = True, elevin =
'w001001 at temp', aspin = 'aspect at temp', slopein = 'slope at temp', glob_rad =
glob_rad, day = x)

    # Combine the radiation file with the previous one:
    if x == 2:
      grass.run_command("r.mapcalculator", overwrite = True, amap =
"global1", bmap = "global2", formula = "A+B",
                      outfile = "suntemp2")
    elif x > 2:
      # Use temp name or else it won't work:
      amap = "suntemp" + str(x-1)
      bmap = "global" + str(x) + "@temp"
      outfile = "suntemp" + str(x)
      grass.run_command("r.mapcalculator", overwrite = True, amap = amap,
bmap = bmap, formula = "A+B",
                          outfile = outfile)

# Finalize:
grass.mapcalc('sun=' + outfile, quiet=False, verbose=True, overwrite=True)
# Clean up:
grass.run_command('g.mremove', flags = 'f', rast = "suntemp*")
grass.run_command('g.mremove', flags = 'f', rast = "global*")
print "Done"

Thanks

Paul

2012/12/28 S. Koukoulas (lists) <sotkouklistes at gmail.com>

> Dear Paul,
>
> having said all that for r.mapcalc, I wonder whether you really need
> mapcalc. In your case,
> couldn't you just use the * r.series * command with Sum (and perhaps
> g.mlist can help you make a list of all the daily maps).
>
> I think that it might be faster as well.
> best,
> sotiris
>
> On 12/28/2012 03:08 PM, Paul Meems wrote:
> > Thanks Sotiris,
> >
> > I've tried your suggestion but converted it to Python because I'm on
> > Windows.
> > This is my script:
> > grass.run_command('g.mremove', rast = ('out*'))
> > grass.mapcalc('sun=global1', quiet=False, verbose=False, overwrite=True)
> > for x in range(2, 10, 1):
> >   exp = 'out=sun+global' + str(x)
> >   print exp
> >   grass.mapcalc(exp, quiet=False, verbose=False, overwrite=True)
> >   grass.run_command('g.rename', flag='o', rast = ('out','sun'))
> >
> > This doesn't result in the expected outcome.
> > sun is equal to global1 and no values are added
> > out is global1+global2 but is no more values are added. But mapcalc is
> > called several times.
> > I've tested this with:
> > r.what --v -f -n input=global1 east_north=156287.225281,383024.737026
> > r.what --v -f -n input=sun east_north=156287.225281,383024.737026
> > r.what --v -f -n input=out east_north=156287.225281,383024.737026
> >
> > Not only am I a new user to GRASS but I'm also still a novice in
> > Python. Most likely it is something in my script.
> >
> > If somebody has a different approach to get the solar radiation of a
> > whole year I'm happy to try that.
> >
> > Thanks,
> >
> > Paul
> >
> >
> >
> >     Date: Fri, 28 Dec 2012 14:36:15 +0200
> >     From: "S. Koukoulas (lists)" <sotkouklistes at gmail.com
> >     <mailto:sotkouklistes at gmail.com>>
> >     To: grass-user at lists.osgeo.org <mailto:grass-user at lists.osgeo.org>
> >     Subject: Re: [GRASS-user] Grass, Python, mapcalc
> >     Message-ID: <50DD923F.806 at gmail.com <mailto:50DD923F.806 at gmail.com>>
> >     Content-Type: text/plain; charset=windows-1252
> >
> >     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
> >
> >
> >
> >
> > _______________________________________________
> > 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/b32839d8/attachment-0001.html>