[GRASS-dev] r.series threshold patch

[Markus Neteler]

Glynn,

to easier operate on incomplete time series from MODIS (and
others), we would like to suggest attached patch. It
adds a threshold to filter out incomplete pixel series
before calling the aggregation function which saves us
to perform extra runs on counting valid pixels and to
post-filter the aggregated results.

Markus
-------------- next part --------------
Index: raster/r.series/description.html
===================================================================
RCS file: /grassrepository/grass6/raster/r.series/description.html,v
retrieving revision 1.12
diff -u -r1.12 description.html
--- raster/r.series/description.html	16 Jul 2007 07:44:31 -0000	1.12
+++ raster/r.series/description.html	16 Aug 2007 16:08:03 -0000
@@ -28,12 +28,16 @@
 With <EM>-n</EM> flag, any cell for which any of the corresponding input cells are
 NULL is automatically set to NULL (NULL propagation). The aggregate function is not
 called, so all methods behave this way with respect to the <EM>-n</EM> flag.
+With <EM>-n</EM> flag, the parameter <EM>threshold</EM> is ignored.
 <P>
-Without <EM>-n</EM> flag, the complete list of inputs for each cell (including
-NULLs) is passed to the aggregate function. Individual aggregates can
-handle data as they choose. Mostly, they just compute the aggregate
-over the non-NULL values, producing a NULL result only if all inputs
-are NULL.
+Without <EM>-n</EM> flag, the parameter <EM>threshold</EM> is used to decide
+when to call the aggregate function: if the list of inputs for each cell
+contains at least <EM>threshold</EM> non NULLs, then the complete list of
+values (including NULLs) is passed over to the aggregate function.
+<P>
+Individual aggregates can handle data as they choose. Mostly, they just compute
+the aggregate over the non-NULL values, producing a NULL result only if all
+inputs are NULL.
 <p>
 The <EM>min_raster</EM> and <EM>max_raster</EM> methods generate a map with the
 number of the raster map that holds the minimum/maximum value of the
Index: raster/r.series/main.c
===================================================================
RCS file: /grassrepository/grass6/raster/r.series/main.c,v
retrieving revision 2.16
diff -u -r2.16 main.c
--- raster/r.series/main.c	23 Jul 2007 15:47:13 -0000	2.16
+++ raster/r.series/main.c	16 Aug 2007 16:08:03 -0000
@@ -77,7 +77,7 @@
 {
 	struct GModule *module;
 	struct {
-		struct Option *input, *output, *method;
+		struct Option *input, *output, *method, *thresh;
 	} parm;
 	struct {
     /* please, remove before GRASS 7 released */
@@ -95,7 +95,7 @@
 	DCELL *out_buf;
 	DCELL *values;
 	int nrows, ncols;
-	int row, col;
+	int row, col, thresh;
 
 	G_gisinit(argv[0]);
 
@@ -117,6 +117,13 @@
 	parm.method->options    = build_method_list();
 	parm.method->description= _("Aggregate operation") ;
 
+	parm.thresh = G_define_option() ;
+	parm.thresh->key        = "threshold" ;
+	parm.thresh->type       = TYPE_INTEGER ;
+	parm.thresh->required   = NO ;
+	parm.thresh->answer     = "1" ;
+	parm.thresh->description= _("Minimum number of non-NULL observations to perform computation") ;
+
 	/* please, remove before GRASS 7 released */
 	flag.quiet = G_define_flag();
 	flag.quiet->key = 'q';
@@ -177,6 +184,7 @@
 
 	/* process the output map */
 	out_name = parm.output->answer;
+	thresh = atoi(parm.thresh->answer);
 
 	out_fd = G_open_raster_new(out_name, DCELL_TYPE);
 	if (out_fd < 0)
@@ -209,12 +217,12 @@
 				DCELL v = inputs[i].buf[col];
 
 				if (G_is_d_null_value(&v))
-					null = 1;
+					null++;
 
 				values[i] = v;
 			}
 
-			if (null && flag.nulls->answer)
+			if ((null && flag.nulls->answer)||((num_inputs-null) < thresh))
 				G_set_d_null_value(&out_buf[col], 1);
 			else
 				(*method_fn)(&out_buf[col], values, num_inputs);