1
2 /****************************************************************************
3 *
4 * MODULE: r.series
5 * AUTHOR(S): Glynn Clements <glynn gclements.plus.com> (original contributor)
6 * Hamish Bowman <hamish_b yahoo.com>, Jachym Cepicky <jachym les-ejk.cz>,
7 * Martin Wegmann <wegmann biozentrum.uni-wuerzburg.de>
8 * PURPOSE:
9 * COPYRIGHT: (C) 2002-2008 by the GRASS Development Team
10 *
11 * This program is free software under the GNU General Public
12 * License (>=v2). Read the file COPYING that comes with GRASS
13 * for details.
14 *
15 *****************************************************************************/
16 #include <string.h>
17 #include <stdlib.h>
18 #include <unistd.h>
19
20 #include <grass/gis.h>
21 #include <grass/raster.h>
22 #include <grass/glocale.h>
23 #include <grass/stats.h>
24
25 struct menu
26 {
27 stat_func *method; /* routine to compute new value */
28 stat_func_w *method_w; /* routine to compute new value (weighted) */
29 RASTER_MAP_TYPE outtype; /* type of result */
30 char *name; /* method name */
31 char *text; /* menu display - full description */
32 } menu[] = {
33 {c_ave, w_ave, DCELL_TYPE, "average", "average value"},
34 {c_count, w_count, CELL_TYPE, "count", "count of non-NULL cells"},
35 {c_median, w_median, DCELL_TYPE, "median", "median value"},
36 {c_mode, w_mode, -1, "mode", "most frequently occurring value"},
37 {c_min, NULL, -1, "minimum", "lowest value"},
38 {c_minx, NULL, CELL_TYPE, "min_raster", "raster with lowest value"},
39 {c_max, NULL, -1, "maximum", "highest value"},
40 {c_maxx, NULL, CELL_TYPE, "max_raster", "raster with highest value"},
41 {c_stddev, w_stddev, DCELL_TYPE, "stddev", "standard deviation"},
42 {c_range, NULL, -1, "range", "range of values"},
43 {c_sum, w_sum, DCELL_TYPE, "sum", "sum of values"},
44 {c_var, w_var, DCELL_TYPE, "variance", "statistical variance"},
45 {c_divr, NULL, CELL_TYPE, "diversity", "number of different values"},
46 {c_reg_m, w_reg_m, DCELL_TYPE, "slope", "linear regression slope"},
47 {c_reg_c, w_reg_c, DCELL_TYPE, "offset", "linear regression offset"},
48 {c_reg_r2, w_reg_r2, DCELL_TYPE, "detcoeff", "linear regression coefficient of determination"},
49 {c_reg_t, w_reg_t, DCELL_TYPE, "tvalue", "linear regression t-value"},
50 {c_quart1, w_quart1, DCELL_TYPE, "quart1", "first quartile"},
51 {c_quart3, w_quart3, DCELL_TYPE, "quart3", "third quartile"},
52 {c_perc90, w_perc90, DCELL_TYPE, "perc90", "ninetieth percentile"},
53 {c_quant, w_quant, DCELL_TYPE, "quantile", "arbitrary quantile"},
54 {c_skew, w_skew, DCELL_TYPE, "skewness", "skewness"},
55 {c_kurt, w_kurt, DCELL_TYPE, "kurtosis", "kurtosis"},
56 {NULL, NULL, 0, NULL, NULL}
57 };
58
59 struct input
60 {
61 const char *name;
62 int fd;
63 DCELL *buf;
64 DCELL weight;
65 };
66
67 struct output
68 {
69 const char *name;
70 int fd;
71 DCELL *buf;
72 stat_func *method_fn;
73 stat_func_w *method_fn_w;
74 double quantile;
75 };
76
build_method_list(void)77 static char *build_method_list(void)
78 {
79 char *buf = G_malloc(1024);
80 char *p = buf;
81 int i;
82
83 for (i = 0; menu[i].name; i++) {
84 char *q;
85
86 if (i)
87 *p++ = ',';
88 for (q = menu[i].name; *q; p++, q++)
89 *p = *q;
90 }
91 *p = '\0';
92
93 return buf;
94 }
95
find_method(const char * method_name)96 static int find_method(const char *method_name)
97 {
98 int i;
99
100 for (i = 0; menu[i].name; i++)
101 if (strcmp(menu[i].name, method_name) == 0)
102 return i;
103
104 G_fatal_error(_("Unknown method <%s>"), method_name);
105
106 return -1;
107 }
108
main(int argc,char * argv[])109 int main(int argc, char *argv[])
110 {
111 struct GModule *module;
112 struct
113 {
114 struct Option *input, *file, *output, *method, *weights, *quantile, *range;
115 } parm;
116 struct
117 {
118 struct Flag *nulls, *lazy;
119 } flag;
120 int i;
121 int num_inputs;
122 struct input *inputs = NULL;
123 int num_outputs;
124 struct output *outputs = NULL;
125 struct History history;
126 DCELL *values = NULL, *values_tmp = NULL;
127 DCELL(*values_w)[2]; /* list of values and weights */
128 DCELL(*values_w_tmp)[2]; /* list of values and weights */
129 int have_weights;
130 int nrows, ncols;
131 int row, col;
132 double lo, hi;
133 RASTER_MAP_TYPE intype, maptype;
134
135 G_gisinit(argv[0]);
136
137 module = G_define_module();
138 G_add_keyword(_("raster"));
139 G_add_keyword(_("aggregation"));
140 G_add_keyword(_("series"));
141 module->description =
142 _("Makes each output cell value a "
143 "function of the values assigned to the corresponding cells "
144 "in the input raster map layers.");
145
146 parm.input = G_define_standard_option(G_OPT_R_INPUTS);
147 parm.input->required = NO;
148
149 parm.file = G_define_standard_option(G_OPT_F_INPUT);
150 parm.file->key = "file";
151 parm.file->description = _("Input file with one raster map name and optional one weight per line, field separator between name and weight is |");
152 parm.file->required = NO;
153
154 parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
155 parm.output->multiple = YES;
156
157 parm.method = G_define_option();
158 parm.method->key = "method";
159 parm.method->type = TYPE_STRING;
160 parm.method->required = YES;
161 parm.method->options = build_method_list();
162 parm.method->description = _("Aggregate operation");
163 parm.method->multiple = YES;
164
165 parm.quantile = G_define_option();
166 parm.quantile->key = "quantile";
167 parm.quantile->type = TYPE_DOUBLE;
168 parm.quantile->required = NO;
169 parm.quantile->description = _("Quantile to calculate for method=quantile");
170 parm.quantile->options = "0.0-1.0";
171 parm.quantile->multiple = YES;
172
173 parm.weights = G_define_option();
174 parm.weights->key = "weights";
175 parm.weights->type = TYPE_DOUBLE;
176 parm.weights->required = NO;
177 parm.weights->description = _("Weighting factor for each input map, default value is 1.0 for each input map");
178 parm.weights->multiple = YES;
179
180 parm.range = G_define_option();
181 parm.range->key = "range";
182 parm.range->type = TYPE_DOUBLE;
183 parm.range->key_desc = "lo,hi";
184 parm.range->description = _("Ignore values outside this range");
185
186 flag.nulls = G_define_flag();
187 flag.nulls->key = 'n';
188 flag.nulls->description = _("Propagate NULLs");
189
190 flag.lazy = G_define_flag();
191 flag.lazy->key = 'z';
192 flag.lazy->description = _("Do not keep files open");
193
194 if (G_parser(argc, argv))
195 exit(EXIT_FAILURE);
196
197 lo = -1.0 / 0.0; /* -inf */
198 hi = 1.0 / 0.0; /* inf */
199 if (parm.range->answer) {
200 lo = atof(parm.range->answers[0]);
201 hi = atof(parm.range->answers[1]);
202 }
203
204 if (parm.input->answer && parm.file->answer)
205 G_fatal_error(_("%s= and %s= are mutually exclusive"),
206 parm.input->key, parm.file->key);
207
208 if (!parm.input->answer && !parm.file->answer)
209 G_fatal_error(_("Please specify %s= or %s="),
210 parm.input->key, parm.file->key);
211
212 have_weights = 0;
213
214 intype = -1;
215
216 /* process the input maps from the file */
217 if (parm.file->answer) {
218 FILE *in;
219 int max_inputs;
220
221 if (strcmp(parm.file->answer, "-") == 0)
222 in = stdin;
223 else {
224 in = fopen(parm.file->answer, "r");
225 if (!in)
226 G_fatal_error(_("Unable to open input file <%s>"), parm.file->answer);
227 }
228
229 num_inputs = 0;
230 max_inputs = 0;
231
232 for (;;) {
233 char buf[GNAME_MAX + 50]; /* Name and weight*/
234 char tok_buf[GNAME_MAX + 50];
235 char *name;
236 int ntokens;
237 char **tokens;
238 struct input *p;
239 double weight = 1.0;
240
241 if (!G_getl2(buf, sizeof(buf), in))
242 break;
243
244 strcpy(tok_buf, buf);
245 tokens = G_tokenize(tok_buf, "|");
246 ntokens = G_number_of_tokens(tokens);
247
248 name = G_chop(tokens[0]);
249 if (ntokens > 1) {
250 weight = atof(G_chop(tokens[1]));
251
252 if (weight < 0)
253 G_fatal_error(_("Weights must be positive"));
254
255 if (weight != 1)
256 have_weights = 1;
257 }
258
259 /* Ignore empty lines */
260 if (!*name)
261 continue;
262
263 if (num_inputs >= max_inputs) {
264 max_inputs += 100;
265 inputs = G_realloc(inputs, max_inputs * sizeof(struct input));
266 }
267 p = &inputs[num_inputs++];
268
269 p->name = G_store(name);
270 p->weight = weight;
271 G_verbose_message(_("Reading raster map <%s> using weight %f..."), p->name, p->weight);
272 p->fd = Rast_open_old(p->name, "");
273 if (p->fd < 0)
274 G_fatal_error(_("Unable to open input raster <%s>"), p->name);
275 maptype = Rast_get_map_type(p->fd);
276 if (intype == -1)
277 intype = maptype;
278 else {
279 if (intype != maptype)
280 intype = DCELL_TYPE;
281 }
282 if (flag.lazy->answer)
283 Rast_close(p->fd);
284 p->buf = Rast_allocate_d_buf();
285 }
286
287 if (num_inputs < 1)
288 G_fatal_error(_("No raster map name found in input file"));
289
290 fclose(in);
291 }
292 else {
293 int num_weights;
294
295 for (i = 0; parm.input->answers[i]; i++)
296 ;
297 num_inputs = i;
298
299 if (num_inputs < 1)
300 G_fatal_error(_("Raster map not found"));
301
302 /* count weights */
303 num_weights = 0;
304 if (parm.weights->answers) {
305 for (i = 0; parm.weights->answers[i]; i++)
306 ;
307 num_weights = i;
308 }
309
310 if (num_weights && num_weights != num_inputs)
311 G_fatal_error(_("input= and weights= must have the same number of values"));
312
313 inputs = G_malloc(num_inputs * sizeof(struct input));
314
315 for (i = 0; i < num_inputs; i++) {
316 struct input *p = &inputs[i];
317
318 p->name = parm.input->answers[i];
319 p->weight = 1.0;
320
321 if (num_weights) {
322 p->weight = (DCELL)atof(parm.weights->answers[i]);
323
324 if (p->weight < 0)
325 G_fatal_error(_("Weights must be positive"));
326
327 if (p->weight != 1)
328 have_weights = 1;
329 }
330
331 G_verbose_message(_("Reading raster map <%s> using weight %f..."), p->name, p->weight);
332 p->fd = Rast_open_old(p->name, "");
333 if (p->fd < 0)
334 G_fatal_error(_("Unable to open input raster <%s>"), p->name);
335 maptype = Rast_get_map_type(p->fd);
336 if (intype == -1)
337 intype = maptype;
338 else {
339 if (intype != maptype)
340 intype = DCELL_TYPE;
341 }
342 if (flag.lazy->answer)
343 Rast_close(p->fd);
344 p->buf = Rast_allocate_d_buf();
345 }
346 }
347
348 /* process the output maps */
349 for (i = 0; parm.output->answers[i]; i++)
350 ;
351 num_outputs = i;
352
353 for (i = 0; parm.method->answers[i]; i++)
354 ;
355 if (num_outputs != i)
356 G_fatal_error(_("output= and method= must have the same number of values"));
357
358 outputs = G_calloc(num_outputs, sizeof(struct output));
359
360 for (i = 0; i < num_outputs; i++) {
361 struct output *out = &outputs[i];
362 const char *output_name = parm.output->answers[i];
363 const char *method_name = parm.method->answers[i];
364 int method = find_method(method_name);
365
366 out->name = output_name;
367
368 if (have_weights) {
369 if (menu[method].method_w) {
370 out->method_fn = NULL;
371 out->method_fn_w = menu[method].method_w;
372 /* special case mode: the result of a weighed mode
373 * can be stored as type of input
374 * all other weighed versions: result as DCELL_TYPE */
375 if (menu[method].outtype == CELL_TYPE)
376 menu[method].outtype = DCELL_TYPE;
377 }
378 else {
379 G_warning(_("Method %s not compatible with weights, using unweighed version instead"),
380 method_name);
381
382 out->method_fn = menu[method].method;
383 out->method_fn_w = NULL;
384 }
385 }
386 else {
387 out->method_fn = menu[method].method;
388 out->method_fn_w = NULL;
389 }
390
391 out->quantile = (parm.quantile->answer && parm.quantile->answers[i])
392 ? atof(parm.quantile->answers[i])
393 : 0;
394 out->buf = Rast_allocate_d_buf();
395 if (menu[method].outtype == -1)
396 out->fd = Rast_open_new(output_name, intype);
397 else
398 out->fd = Rast_open_new(output_name, menu[method].outtype);
399 }
400
401 /* initialise variables */
402 values = G_malloc(num_inputs * sizeof(DCELL));
403 values_tmp = G_malloc(num_inputs * sizeof(DCELL));
404 values_w = NULL;
405 values_w_tmp = NULL;
406 if (have_weights) {
407 values_w = (DCELL(*)[2]) G_malloc(num_inputs * 2 * sizeof(DCELL));
408 values_w_tmp = (DCELL(*)[2]) G_malloc(num_inputs * 2 * sizeof(DCELL));
409 }
410
411 nrows = Rast_window_rows();
412 ncols = Rast_window_cols();
413
414 /* process the data */
415 G_verbose_message(_("Percent complete..."));
416
417 for (row = 0; row < nrows; row++) {
418 G_percent(row, nrows, 2);
419
420 if (flag.lazy->answer) {
421 /* Open the files only on run time */
422 for (i = 0; i < num_inputs; i++) {
423 inputs[i].fd = Rast_open_old(inputs[i].name, "");
424 Rast_get_d_row(inputs[i].fd, inputs[i].buf, row);
425 Rast_close(inputs[i].fd);
426 }
427 }
428 else {
429 for (i = 0; i < num_inputs; i++)
430 Rast_get_d_row(inputs[i].fd, inputs[i].buf, row);
431 }
432
433 for (col = 0; col < ncols; col++) {
434 int null = 0;
435
436 for (i = 0; i < num_inputs; i++) {
437 DCELL v = inputs[i].buf[col];
438
439 if (Rast_is_d_null_value(&v))
440 null = 1;
441 else if (parm.range->answer && (v < lo || v > hi)) {
442 Rast_set_d_null_value(&v, 1);
443 null = 1;
444 }
445 values[i] = v;
446 if (have_weights) {
447 values_w[i][0] = v;
448 values_w[i][1] = inputs[i].weight;
449 }
450 }
451
452 for (i = 0; i < num_outputs; i++) {
453 struct output *out = &outputs[i];
454
455 if (null && flag.nulls->answer)
456 Rast_set_d_null_value(&out->buf[col], 1);
457 else {
458 if (out->method_fn_w) {
459 memcpy(values_w_tmp, values_w, num_inputs * 2 * sizeof(DCELL));
460 (*out->method_fn_w)(&out->buf[col], values_w_tmp, num_inputs, &out->quantile);
461 }
462 else {
463 memcpy(values_tmp, values, num_inputs * sizeof(DCELL));
464 (*out->method_fn)(&out->buf[col], values_tmp, num_inputs, &out->quantile);
465 }
466 }
467 }
468 }
469
470 for (i = 0; i < num_outputs; i++)
471 Rast_put_d_row(outputs[i].fd, outputs[i].buf);
472 }
473
474 G_percent(row, nrows, 2);
475
476 /* close output maps */
477 for (i = 0; i < num_outputs; i++) {
478 struct output *out = &outputs[i];
479
480 Rast_close(out->fd);
481
482 Rast_short_history(out->name, "raster", &history);
483 Rast_command_history(&history);
484 Rast_write_history(out->name, &history);
485 }
486
487 /* Close input maps */
488 if (!flag.lazy->answer) {
489 for (i = 0; i < num_inputs; i++)
490 Rast_close(inputs[i].fd);
491 }
492
493 exit(EXIT_SUCCESS);
494 }
495