[GRASS-dev] eval option in r.mapcalc in a python script

GRASS GIS GRASS GIS developer
1

From the r.mapcalc manual, the 'eval' option is described to run more complex expressions. But I wonder how I apply the example below in a python script?

r.mapcalc << EOF
eval(elev_200 = elevation - 200, \
elev_5 = 5 * elevation, \
elev_p = pow(elev_5, 2))
elevation_result = (0.5 * elev_200) + 0.8 * elev_p
EOF

Somewhat related, if I have a large number of maps. For each map I need to compute a value. At the end, the resulting maps need to be summed up. One way is to run a loop, and use r.series to sum up the results. E.g.,

for i in xrange(len(IN)):
  out = 'shi' + str(i)
  out.append(out)
  grass.mapcalc("$out = log($in)", out=out[i], in=in[i])
grass.run_command("r.series", output="endresult", input=out, method="sum")

This generates a lot of intermediate layers, which can potentially take up a lot of space on HD. To avoid that, I could do

grass.mapcalc("$out = 0", out=out)
for i in xrange(len(IN)):
  grass.mapcalc("$out = $out + log($in)", out=out, in=in[i]), overwrite=True

An alternative is to construct one large expression along the lines of expression = "log(in[1]) + log(in[2] + ....". I assume this will be faster? It would result in a unwieldy large expression though, and I remember something about a limit to the length of the expression? So that is why I wonder if using the 'eval' option would be useful.

Any tips for the best approach?

2

On 14/01/16 09:34, Paulo van Breugel wrote:

From the r.mapcalc manual, the 'eval' option is described to run more
complex expressions. But I wonder how I apply the example below in a
python script?

eval(elev_200 = elevation - 200, \
elev_5 = 5 * elevation, \
elev_p = pow(elev_5, 2))
elevation_result = (0.5 * elev_200) + 0.8 * elev_p

Can't you just put the above in a string variable and feed it to the 'expression' parameter of grass.mapcalc (or alternatives grass.run_command('r.mapcalc', ...)?

Somewhat related, if I have a large number of maps. For each map I need
to compute a value. At the end, the resulting maps need to be summed up.
One way is to run a loop, and use r.series to sum up the results. E.g.,

for i in xrange(len(IN)):
     out = 'shi' + str(i)
     out.append(out)
     grass.mapcalc("$out = log($in)", out=out[i], in=in[i])
grass.run_command("r.series", output="endresult", input=out, method="sum")

This generates a lot of intermediate layers, which can potentially take
up a lot of space on HD. To avoid that, I could do

grass.mapcalc("$out = 0", out=out)
for i in xrange(len(IN)):
     grass.mapcalc("$out = $out + log($in)", out=out, in=in[i]),
overwrite=True

ISTR that modifying the contents of a map in-place might be possible, but not recommended.

As a compromose you could do something like this:

grass.mapcalc("$oldout = 0")
for i in xrange(len(IN)):
      grass.mapcalc("$newout = $oldout + log($in)", out=out, in=in[i])
      grass.run_command('g.rename', '$newout,$oldout', overwrite=True)

An alternative is to construct one large expression along the lines of
expression = "log(in[1]) + log(in[2] + ....". I assume this will be
faster? It would result in a unwieldy large expression though, and I
remember something about a limit to the length of the expression? So
that is why I wonder if using the 'eval' option would be useful.

Any tips for the best approach?

Maybe one option, as I could image others to have this kind of issue, would be to add a "modifier" (or similar name) parameter to r.series which would allow to calculate the statistic defined by "method" on a modified version of the maps (e.g. squared, square root, log, sin, cos, etc.).

Moritz

3

On 14-01-16 09:57, Moritz Lennert wrote:

On 14/01/16 09:34, Paulo van Breugel wrote:

From the r.mapcalc manual, the 'eval' option is described to run more
complex expressions. But I wonder how I apply the example below in a
python script?

eval(elev_200 = elevation - 200, \
elev_5 = 5 * elevation, \
elev_p = pow(elev_5, 2))
elevation_result = (0.5 * elev_200) + 0.8 * elev_p

Can't you just put the above in a string variable and feed it to the 'expression' parameter of grass.mapcalc (or alternatives grass.run_command('r.mapcalc', ...)?

Thanks, I'll try that

Somewhat related, if I have a large number of maps. For each map I need
to compute a value. At the end, the resulting maps need to be summed up.
One way is to run a loop, and use r.series to sum up the results. E.g.,

for i in xrange(len(IN)):
     out = 'shi' + str(i)
     out.append(out)
     grass.mapcalc("$out = log($in)", out=out[i], in=in[i])
grass.run_command("r.series", output="endresult", input=out, method="sum")

This generates a lot of intermediate layers, which can potentially take
up a lot of space on HD. To avoid that, I could do

grass.mapcalc("$out = 0", out=out)
for i in xrange(len(IN)):
     grass.mapcalc("$out = $out + log($in)", out=out, in=in[i]),
overwrite=True

ISTR that modifying the contents of a map in-place might be possible, but not recommended.

I think your solution below is perfectly fine, but out of curiosity, and to get a better grasp of the possible pitfalls, what is the reason the above in-place modification is not recommended?

As a compromose you could do something like this:

grass.mapcalc("$oldout = 0")
for i in xrange(len(IN)):
     grass.mapcalc("$newout = $oldout + log($in)", out=out, in=in[i])
     grass.run_command('g.rename', '$newout,$oldout', overwrite=True)

An alternative is to construct one large expression along the lines of
expression = "log(in[1]) + log(in[2] + ....". I assume this will be
faster? It would result in a unwieldy large expression though, and I
remember something about a limit to the length of the expression? So
that is why I wonder if using the 'eval' option would be useful.

Any tips for the best approach?

Maybe one option, as I could image others to have this kind of issue, would be to add a "modifier" (or similar name) parameter to r.series which would allow to calculate the statistic defined by "method" on a modified version of the maps (e.g. squared, square root, log, sin, cos, etc.).

The log in the example above was for illustrative purposes, the real expression is a bit more complex. Nevertheless, I think it would be very useful to have such a modifier, really good idea!

Moritz

4

On 14/01/16 13:05, Paulo van Breugel wrote:

grass.mapcalc("$out = 0", out=out)
for i in xrange(len(IN)):
     grass.mapcalc("$out = $out + log($in)", out=out, in=in[i]),
overwrite=True

ISTR that modifying the contents of a map in-place might be possible,
but not recommended.

I think your solution below is perfectly fine, but out of curiosity, and
to get a better grasp of the possible pitfalls, what is the reason the
above in-place modification is not recommended?

See Glynn's comment (to you :wink: ) here http://lists.osgeo.org/pipermail/grass-dev/2015-December/077600.html

Maybe one option, as I could image others to have this kind of issue,
would be to add a "modifier" (or similar name) parameter to r.series
which would allow to calculate the statistic defined by "method" on a
modified version of the maps (e.g. squared, square root, log, sin,
cos, etc.).

The log in the example above was for illustrative purposes, the real
expression is a bit more complex. Nevertheless, I think it would be very
useful to have such a modifier, really good idea!

Try the attached patch, which is just a brute-force, proof-of-concept. It would be nicer to solve this with function pointers, but I'm not at ease with that, so I'll leave that to others.

If you would like to see this implemented, I suggest filing an enhancement ticket. You can attach this patch if you want to.

Moritz

(attachments)

r_series_modifier.diff (1.75 KB)

5

I'm sorry, I haven't had time to update r.mapcalc documentation to
clarify things in a nice way. Here's part of discussion on this topic:
https://lists.osgeo.org/pipermail/grass-dev/2015-November/077491.html

Modifying rasters in-place is not possible also because raster library
does not support it. It would also posses a problem for r.mapcalc, as
it supports neighbour operator - what should be an outcome if you
request data from -1 row (the one you just modified)?

2016-01-14 14:05 GMT+02:00 Paulo van Breugel <p.vanbreugel@gmail.com>:

ISTR that modifying the contents of a map in-place might be possible, but
not recommended.

I think your solution below is perfectly fine, but out of curiosity, and to
get a better grasp of the possible pitfalls, what is the reason the above
in-place modification is not recommended?

6

On 14-01-16 13:54, Moritz Lennert wrote:

On 14/01/16 13:05, Paulo van Breugel wrote:

grass.mapcalc("$out = 0", out=out)
for i in xrange(len(IN)):
     grass.mapcalc("$out = $out + log($in)", out=out, in=in[i]),
overwrite=True

ISTR that modifying the contents of a map in-place might be possible,
but not recommended.

I think your solution below is perfectly fine, but out of curiosity, and
to get a better grasp of the possible pitfalls, what is the reason the
above in-place modification is not recommended?

See Glynn's comment (to you :wink: ) here http://lists.osgeo.org/pipermail/grass-dev/2015-December/077600.html

OK, that is embarrassing (and sorry Glynn, I do try to pay attention). I did remember having been point out to this before, but forgot the details. But yes, I should have know better than to fall into the same pitfall.

Maybe one option, as I could image others to have this kind of issue,
would be to add a "modifier" (or similar name) parameter to r.series
which would allow to calculate the statistic defined by "method" on a
modified version of the maps (e.g. squared, square root, log, sin,
cos, etc.).

The log in the example above was for illustrative purposes, the real
expression is a bit more complex. Nevertheless, I think it would be very
useful to have such a modifier, really good idea!

Try the attached patch, which is just a brute-force, proof-of-concept. It would be nicer to solve this with function pointers, but I'm not at ease with that, so I'll leave that to others.

If you would like to see this implemented, I suggest filing an enhancement ticket. You can attach this patch if you want to.

Just create a ticket, many thanks for the patch and helping out on this

Moritz

7

Paulo van Breugel wrote:

From the r.mapcalc manual, the 'eval' option is described to run more
complex expressions. But I wonder how I apply the example below in a
python script?

r.mapcalc << EOF
eval(elev_200 = elevation - 200, \
elev_5 = 5 * elevation, \
elev_p = pow(elev_5, 2))
elevation_result = (0.5 * elev_200) + 0.8 * elev_p
EOF

Use grass.script.mapcalc().

The first parameter is a template string, which is passed to
string.Template() (see the Python documentation for the string module)
then instantiated using any additional keyword arguments. This can be
used if you need to be able to substitute parameters into the
expression.

For complex expressions, you can use any of Python's features to
construct the expression. E.g.

        assignments = [
            "elev_200 = elevation - 200",
            "elev_5 = 5 * elevation",
            "elev_p = pow(elev_5, 2)"]
        result = "elevation_result = (0.5 * elev_200) + 0.8 * elev_p"
        expr = "eval(" + ",".join(assignments) + ");" + result
        gscript.mapcalc(expr)

This may be useful if you need to be able to dynamically manipulate
the structure of the expression rather than just parameters.

--
Glynn Clements <glynn@gclements.plus.com>