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

[Glynn Clements]

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 
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

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

[Paulo van Breugel]

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?
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

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

[Moritz Lennert]

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
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

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

[Moritz Lennert]

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 ;-) ) 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
Index: raster/r.series/main.c
===
--- raster/r.series/main.c	(révision 67579)
+++ raster/r.series/main.c	(copie de travail)
@@ -16,6 +16,7 @@
 #include 
 #include 
 #include 
+#include 
 
 #include 
 #include 
@@ -109,7 +110,7 @@
 struct GModule *module;
 struct
 {
-	struct Option *input, *file, *output, *method, *weights, *quantile, *range;
+	struct Option *input, *file, *output, *method, *weights, *quantile, *range, *modifier;
 } parm;
 struct
 {
@@ -178,6 +179,13 @@
 parm.range->key_desc = "lo,hi";
 parm.range->description = _("Ignore values outside this range");
 
+parm.modifier = G_define_option(); 
+parm.modifier->key = "modifier";
+parm.modifier->type = TYPE_STRING;
+parm.modifier->options = "log,square,square_root";
+parm.modifier->multiple = NO;
+parm.modifier->description = _("Function to apply to map value before applying method");
+
 flag.nulls = G_define_flag();
 flag.nulls->key = 'n';
 flag.nulls->description = _("Propagate NULLs");
@@ -370,7 +378,16 @@
 		Rast_set_d_null_value(, 1);
 		null = 1;
 		}
-		values[i] = v * inputs[i].weight;
+		if (parm.modifier->answer) { 
+		if (strcmp(parm.modifier->answer, "log") == 0) 
+		 	values[i] = log(v) * inputs[i].weight;
+		if (strcmp(parm.modifier->answer, "square") == 0) 
+		 	values[i] = pow(v, 2) * inputs[i].weight;
+		if (strcmp(parm.modifier->answer, "square_root") == 0) 
+		 	values[i] = sqrt(v) * inputs[i].weight;
+		}
+		else
+	 	values[i] = v * inputs[i].weight;
 	}
 
 	for (i = 0; i < num_outputs; i++) {
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

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

[Paulo van Breugel]

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


___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

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

[Paulo van Breugel]

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 ;-) ) 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


___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

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

[Maris Nartiss]

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 :
>> 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?
___
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev