[GRASS-dev] graph function in r.mapcalc

Hi,

I am using the graph function in r.mapcalc. The input is the name of the map to be converted and a string with XY values, like:

“newmap = graph(map, 1, x1,y1, x2,y2,… xi,yi)”

Often, X and Y values are available as separate columns or vectors. In such cases, it would be much easier if X and Y values can be given as separate vectors, e.g., something like:

“newmap = graph(map, x=x1,x2,x3,x4,…xi, y=y1,y2,y3,y4,…,yi)”

For me the second option would make it especially much easier to create this graphs within scripts. So I wondering, what are the reasons that makes the current way of doing this better? And would it be possible to implement a new option to give the x and y values as separate vectors?

Cheers,

Paulo

Paulo van Breugel wrote:

I am using the graph function in r.mapcalc. The input is the name of the
map to be converted and a string with XY values, like:

"newmap = graph(map, 1, x1,y1, x2,y2,... xi,yi)"

The second argument ("1") shouldn't be there.

Often, X and Y values are available as separate columns or vectors. In
such cases, it would be much easier if X and Y values can be given as
separate vectors, e.g., something like:

"newmap = graph(map, x=x1,x2,x3,x4,...xi, y=y1,y2,y3,y4,...,yi)"

That exact syntax (i.e. with "x=" and "y=" markers) isn't possible
without fundamentally re-writing r.mapcalc.

However, it would be possible to implement:

  newmap = graph(map, x1,x2,x3,x4, y1,y2,y3,y4)

This would boil down to cloning f_graph() in raster/r.mapcalc/xgraph.c
with

  #define X(j) (argz[2 + 2 * (j) + 0][i])
  #define Y(j) (argz[2 + 2 * (j) + 1][i])

changed to:

  #define X(j) (argz[2 + (j) + 0][i])
  #define Y(j) (argz[2 + (j) + n][i])

In practice, we'd want to re-factor the common code.

A (untested) patch to implement a graph2() function with the above
syntax is attached.

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

(attachments)

graph2.diff (3.71 KB)

Hi Glynn

It works perfectly in my (limited) testing, great! I am impressed.

I am not sure what you mean by “re-factor the common code”, but I hope it will be possible to include this graph2() function as a standard as it makes, i.m.h.o. the function more user-friendly.

Cheers,

Paulo

p.s. I applied the patch to my GRASS7, can I also apply the patch to GRASS64?

···

On Mon, Nov 19, 2012 at 4:46 AM, Glynn Clements <glynn@gclements.plus.com> wrote:

Paulo van Breugel wrote:

I am using the graph function in r.mapcalc. The input is the name of the
map to be converted and a string with XY values, like:

“newmap = graph(map, 1, x1,y1, x2,y2,… xi,yi)”

The second argument (“1”) shouldn’t be there.

Often, X and Y values are available as separate columns or vectors. In
such cases, it would be much easier if X and Y values can be given as
separate vectors, e.g., something like:

“newmap = graph(map, x=x1,x2,x3,x4,…xi, y=y1,y2,y3,y4,…,yi)”

That exact syntax (i.e. with “x=” and “y=” markers) isn’t possible
without fundamentally re-writing r.mapcalc.

However, it would be possible to implement:

newmap = graph(map, x1,x2,x3,x4, y1,y2,y3,y4)

This would boil down to cloning f_graph() in raster/r.mapcalc/xgraph.c
with

#define X(j) (argz[2 + 2 * (j) + 0][i])
#define Y(j) (argz[2 + 2 * (j) + 1][i])

changed to:

#define X(j) (argz[2 + (j) + 0][i])
#define Y(j) (argz[2 + (j) + n][i])

In practice, we’d want to re-factor the common code.

A (untested) patch to implement a graph2() function with the above
syntax is attached.


Glynn Clements <glynn@gclements.plus.com>

Index: raster/r.mapcalc/r.mapcalc.html

— raster/r.mapcalc/r.mapcalc.html (revision 53894)
+++ raster/r.mapcalc/r.mapcalc.html (working copy)
@@ -298,6 +298,8 @@
exp(x,y) x to the power y F
float(x) convert x to single-precision floating point F
graph(x,x1,y1[x2,y2…]) convert the x to a y based on points in a graph F
+graph2(x,x1[,x2,…],y1[,y2…])

  • alternative form of graph() F
    if decision options: *
    if(x) 1 if x not zero, 0 otherwise
    if(x,a) a if x not zero, 0 otherwise
    Index: raster/r.mapcalc/function.c
    ===================================================================
    — raster/r.mapcalc/function.c (revision 53894)
    +++ raster/r.mapcalc/function.c (working copy)
    @@ -74,6 +74,7 @@
    {“nmode”, c_varop, f_nmode},

{“graph”, c_graph, f_graph},

  • {“graph2”, c_graph, f_graph2},

{“rand”, c_binop, f_rand},

Index: raster/r.mapcalc/xgraph.c

— raster/r.mapcalc/xgraph.c (revision 53894)
+++ raster/r.mapcalc/xgraph.c (working copy)
@@ -90,6 +90,9 @@

break;
}
+#undef X
+#undef Y
+#undef x

continue;

@@ -99,3 +102,79 @@

return 0;
}
+
+int f_graph2(int argc, const int *argt, void **args)
+{

  • DCELL **argz = (DCELL **) args;
  • DCELL *res = argz[0];
  • int n = (argc - 1) / 2;
  • int i, j;
  • if (argc < 3)
  • return E_ARG_LO;
  • if (argc % 2 == 0)
  • return E_ARG_NUM;
  • if (argt[0] != DCELL_TYPE)
  • return E_RES_TYPE;
  • for (i = 1; i <= argc; i++)
  • if (argt[i] != DCELL_TYPE)
  • return E_ARG_TYPE;
  • for (i = 0; i < columns; i++) {
    +#define X(j) (argz[2 + (j) + 0][i])
    +#define Y(j) (argz[2 + (j) + n][i])
    +#define x (argz[1][i])
  • if (IS_NULL_D(&x))
  • goto null;
  • for (j = 0; j < n; j++)
  • if (IS_NULL_D(&X(j)))
  • goto null;
  • for (j = 0; j < n - 1; j++)
  • if (X(j + 1) <= X(j))
  • goto null;
  • if (x <= X(0)) {
  • if (IS_NULL_D(&Y(0)))
  • goto null;
  • res[i] = Y(0);
  • continue;
  • }
  • if (x >= X(n - 1)) {
  • if (IS_NULL_D(&Y(n - 1)))
  • goto null;
  • res[i] = Y(n - 1);
  • continue;
  • }
  • for (j = 0; j < n - 1; j++) {
  • if (x > X(j + 1))
  • continue;
  • if (IS_NULL_D(&Y(j)) || IS_NULL_D(&Y(j + 1)))
  • goto null;
  • res[i] =
  • Y(j) + (x - X(j)) * (Y(j + 1) - Y(j)) / (X(j + 1) - X(j));
  • break;
  • }
    +#undef X
    +#undef Y
    +#undef x
  • continue;
  • null:
  • SET_NULL_D(&res[i]);
  • }
  • return 0;
    +}

Index: raster/r.mapcalc/r3.mapcalc.html

— raster/r.mapcalc/r3.mapcalc.html (revision 53894)
+++ raster/r.mapcalc/r3.mapcalc.html (working copy)
@@ -202,6 +202,8 @@
exp(x,y) x to the power y F
float(x) convert x to single-precision floating point F
graph(x,x1,y1[x2,y2…]) convert the x to a y based on points in a graph F
+graph2(x,x1[,x2,…],y1[,y2…])

  • alternative form of graph() F
    if decision options: *
    if(x) 1 if x not zero, 0 otherwise
    if(x,a) a if x not zero, 0 otherwise
    Index: raster/r.mapcalc/func_proto.h
    ===================================================================
    — raster/r.mapcalc/func_proto.h (revision 53894)
    +++ raster/r.mapcalc/func_proto.h (working copy)
    @@ -74,6 +74,7 @@
    extern args_t c_isnull;

extern func_t f_graph;
+extern func_t f_graph2;
extern args_t c_graph;

extern func_t f_min;

Paulo van Breugel wrote:

It works perfectly in my (limited) testing, great! I am impressed.

I am not sure what you mean by "re-factor the common code", but I hope it
will be possible to include this graph2() function as a standard as it
makes, i.m.h.o. the function more user-friendly.

Can you add a "wish" to the bug tracker?

p.s. I applied the patch to my GRASS7, can I also apply the patch to
GRASS64?

I'm not sure that the patch will apply to 6.4, but there's no
fundamental reason the same change cannot be made.

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

On Mon, Nov 19, 2012 at 7:07 PM, Glynn Clements <glynn@gclements.plus.com>wrote:

Paulo van Breugel wrote:

> It works perfectly in my (limited) testing, great! I am impressed.
>
> I am not sure what you mean by "re-factor the common code", but I hope it
> will be possible to include this graph2() function as a standard as it
> makes, i.m.h.o. the function more user-friendly.

Can you add a "wish" to the bug tracker?

Done

> p.s. I applied the patch to my GRASS7, can I also apply the patch to
> GRASS64?

I'm not sure that the patch will apply to 6.4, but there's no
fundamental reason the same change cannot be made.

OK, thanks, I will try out

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