[GRASS-dev] graph function in r.mapcalc
Paulo van Breugel
p.vanbreugel at gmail.com
Mon Nov 19 00:40:18 PST 2012
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 at 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 at 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;
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20121119/c1ced91d/attachment-0001.html>
More information about the grass-dev
mailing list