← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to raster/r.series/main.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.series/main.c
18
18
#include <unistd.h>
19
19
 
20
20
#include <grass/gis.h>
 
21
#include <grass/raster.h>
21
22
#include <grass/glocale.h>
22
23
#include <grass/stats.h>
23
24
 
39
40
    {c_stddev, 0, "stddev",     "standard deviation"},
40
41
    {c_range,  0, "range",      "range of values"},
41
42
    {c_sum,    0, "sum",        "sum of values"},
42
 
    {c_thresh, 1, "threshold",  "threshold value"},
43
43
    {c_var,    0, "variance",   "statistical variance"},
44
44
    {c_divr,   1, "diversity",  "number of different values"},
45
45
    {c_reg_m,  0, "slope",      "linear regression slope"},
46
46
    {c_reg_c,  0, "offset",     "linear regression offset"},
47
47
    {c_reg_r2, 0, "detcoeff",   "linear regression coefficient of determination"},
 
48
    {c_reg_t,  0, "tvalue",     "linear regression t-value"},
48
49
    {c_quart1, 0, "quart1",     "first quartile"},
49
50
    {c_quart3, 0, "quart3",     "third quartile"},
50
51
    {c_perc90, 0, "perc90",     "ninetieth percentile"},
56
57
 
57
58
struct input
58
59
{
59
 
    const char *name, *mapset;
 
60
    const char *name;
60
61
    int fd;
61
62
    DCELL *buf;
 
63
    DCELL weight;
62
64
};
63
65
 
64
66
struct output
68
70
    DCELL *buf;
69
71
    stat_func *method_fn;
70
72
    double quantile;
71
 
    double threshold;
72
73
};
73
74
 
74
75
static char *build_method_list(void)
108
109
    struct GModule *module;
109
110
    struct
110
111
    {
111
 
        struct Option *input, *output, *method, *quantile, *threshold, *range;
 
112
        struct Option *input, *file, *output, *method, *weights, *quantile, *range;
112
113
    } parm;
113
114
    struct
114
115
    {
115
 
        /* please, remove before GRASS 7 released */
116
 
        struct Flag *quiet;
117
 
        struct Flag *nulls;
 
116
        struct Flag *nulls, *lazy;
118
117
    } flag;
119
118
    int i;
120
119
    int num_inputs;
121
 
    struct input *inputs;
 
120
    int num_weights;
 
121
    struct input *inputs = NULL;
122
122
    int num_outputs;
123
 
    struct output *outputs;
 
123
    struct output *outputs = NULL;
124
124
    struct History history;
125
 
    DCELL *values, *values_tmp;
 
125
    DCELL *values = NULL, *values_tmp = NULL;
126
126
    int nrows, ncols;
127
127
    int row, col;
128
128
    double lo, hi;
130
130
    G_gisinit(argv[0]);
131
131
 
132
132
    module = G_define_module();
133
 
    module->keywords = _("raster, series");
 
133
    G_add_keyword(_("raster"));
 
134
    G_add_keyword(_("aggregation"));
 
135
    G_add_keyword(_("series"));
134
136
    module->description =
135
137
        _("Makes each output cell value a "
136
138
          "function of the values assigned to the corresponding cells "
137
139
          "in the input raster map layers.");
138
140
 
139
141
    parm.input = G_define_standard_option(G_OPT_R_INPUTS);
 
142
    parm.input->required = NO;
 
143
 
 
144
    parm.file = G_define_standard_option(G_OPT_F_INPUT);
 
145
    parm.file->key = "file";
 
146
    parm.file->description = _("Input file with one raster map name and optional one weight per line, field separator between name and weight is |");
 
147
    parm.file->required = NO;
140
148
 
141
149
    parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
142
150
    parm.output->multiple = YES;
157
165
    parm.quantile->options = "0.0-1.0";
158
166
    parm.quantile->multiple = YES;
159
167
 
160
 
    parm.threshold = G_define_option();
161
 
    parm.threshold->key = "threshold";
162
 
    parm.threshold->type = TYPE_DOUBLE;
163
 
    parm.threshold->required = NO;
164
 
    parm.threshold->description = _("Threshold to calculate for method=threshold");
165
 
    parm.threshold->multiple = YES;
166
 
 
167
 
    /* please, remove before GRASS 7 released */
168
 
    flag.quiet = G_define_flag();
169
 
    flag.quiet->key = 'q';
170
 
    flag.quiet->description = _("Run quietly");
171
 
 
 
168
    parm.weights = G_define_option();
 
169
    parm.weights->key = "weights";
 
170
    parm.weights->type = TYPE_DOUBLE;
 
171
    parm.weights->required = NO;
 
172
    parm.weights->description = _("Weighting factor for each input map, default value is 1.0 for each input map");
 
173
    parm.weights->multiple = YES;
 
174
    
172
175
    parm.range = G_define_option();
173
176
    parm.range->key = "range";
174
177
    parm.range->type = TYPE_DOUBLE;
179
182
    flag.nulls->key = 'n';
180
183
    flag.nulls->description = _("Propagate NULLs");
181
184
 
 
185
    flag.lazy = G_define_flag();
 
186
    flag.lazy->key = 'z';
 
187
    flag.lazy->description = _("Do not keep files open");
 
188
 
182
189
    if (G_parser(argc, argv))
183
190
        exit(EXIT_FAILURE);
184
191
 
185
 
    /* please, remove before GRASS 7 released */
186
 
    if (flag.quiet->answer) {
187
 
        putenv("GRASS_VERBOSE=0");
188
 
        G_warning(_("The '-q' flag is superseded and will be removed "
189
 
                    "in future. Please use '--quiet' instead."));
190
 
    }
191
 
 
192
192
    if (parm.range->answer) {
193
193
        lo = atof(parm.range->answers[0]);
194
194
        hi = atof(parm.range->answers[1]);
195
195
    }
196
 
 
197
 
    /* process the input maps */
198
 
    for (i = 0; parm.input->answers[i]; i++) ;
199
 
    num_inputs = i;
200
 
 
201
 
    if (num_inputs < 1)
202
 
        G_fatal_error(_("No input raster map(s) specified"));
203
 
 
204
 
    inputs = G_malloc(num_inputs * sizeof(struct input));
205
 
 
206
 
    for (i = 0; i < num_inputs; i++) {
207
 
        struct input *p = &inputs[i];
208
 
 
209
 
        p->name = parm.input->answers[i];
210
 
        p->mapset = G_find_cell2(p->name, "");
211
 
        if (!p->mapset)
212
 
            G_fatal_error(_("Raster map <%s> not found"), p->name);
213
 
        else
214
 
            G_message(_("Reading raster map <%s>..."), p->name);
215
 
        p->fd = G_open_cell_old(p->name, p->mapset);
216
 
        if (p->fd < 0)
217
 
            G_fatal_error(_("Unable to open raster map <%s> in mapset <%s>"),
218
 
                          p->name, p->mapset);
219
 
        p->buf = G_allocate_d_raster_buf();
 
196
    
 
197
    if (parm.input->answer && parm.file->answer)
 
198
        G_fatal_error(_("%s= and %s= are mutually exclusive"),
 
199
                        parm.input->key, parm.file->key);
 
200
 
 
201
    if (!parm.input->answer && !parm.file->answer)
 
202
        G_fatal_error(_("Please specify %s= or %s="),
 
203
                        parm.input->key, parm.file->key);
 
204
 
 
205
 
 
206
    /* process the input maps from the file */
 
207
    if (parm.file->answer) {
 
208
        FILE *in;
 
209
        int max_inputs;
 
210
 
 
211
        if (strcmp(parm.file->answer, "-") == 0)
 
212
            in = stdin;
 
213
        else {
 
214
            in = fopen(parm.file->answer, "r");
 
215
            if (!in)
 
216
                G_fatal_error(_("Unable to open input file <%s>"), parm.file->answer);
 
217
        }
 
218
    
 
219
        num_inputs = 0;
 
220
        max_inputs = 0;
 
221
 
 
222
        for (;;) {
 
223
            char buf[GNAME_MAX + 50]; /* Name and weight*/
 
224
            char tok_buf[GNAME_MAX + 50];
 
225
            char *name;
 
226
            int ntokens;
 
227
            char **tokens;
 
228
            struct input *p;
 
229
            double weight = 1.0;
 
230
 
 
231
            if (!G_getl2(buf, sizeof(buf), in))
 
232
                break;
 
233
 
 
234
            strcpy(tok_buf, buf);
 
235
            tokens = G_tokenize(tok_buf, "|");
 
236
            ntokens = G_number_of_tokens(tokens);
 
237
 
 
238
            if(ntokens > 1) {
 
239
                name = G_chop(tokens[0]);
 
240
                weight = atof(G_chop(tokens[1]));
 
241
            } else {
 
242
                name = G_chop(buf);
 
243
            }
 
244
 
 
245
            /* Ignore empty lines */
 
246
            if (!*name)
 
247
                continue;
 
248
 
 
249
            if (num_inputs >= max_inputs) {
 
250
                max_inputs += 100;
 
251
                inputs = G_realloc(inputs, max_inputs * sizeof(struct input));
 
252
            }
 
253
            p = &inputs[num_inputs++];
 
254
 
 
255
            p->name = G_store(name);
 
256
            p->weight = weight;
 
257
            G_verbose_message(_("Reading raster map <%s> using weight %f..."), p->name, p->weight);
 
258
            p->buf = Rast_allocate_d_buf();
 
259
            if (!flag.lazy->answer)
 
260
                p->fd = Rast_open_old(p->name, "");
 
261
        }
 
262
 
 
263
        if (num_inputs < 1)
 
264
            G_fatal_error(_("No raster map name found in input file"));
 
265
        
 
266
        fclose(in);
 
267
    }
 
268
    else {
 
269
        for (i = 0; parm.input->answers[i]; i++)
 
270
            ;
 
271
        num_inputs = i;
 
272
 
 
273
        if (num_inputs < 1)
 
274
            G_fatal_error(_("Raster map not found"));
 
275
 
 
276
        /* count weights */
 
277
        if(parm.weights->answers) {
 
278
            for (i = 0; parm.weights->answers[i]; i++)
 
279
                    ;
 
280
            num_weights = i;
 
281
        } else {
 
282
            num_weights = 0;
 
283
        }
 
284
    
 
285
        if (num_weights && num_weights != num_inputs)
 
286
                G_fatal_error(_("input= and weights= must have the same number of values"));
 
287
        
 
288
        inputs = G_malloc(num_inputs * sizeof(struct input));
 
289
 
 
290
        for (i = 0; i < num_inputs; i++) {
 
291
            struct input *p = &inputs[i];
 
292
 
 
293
            p->name = parm.input->answers[i];
 
294
 
 
295
            if(num_weights)
 
296
                p->weight = (DCELL)atof(parm.weights->answers[i]);
 
297
            else
 
298
                p->weight = 1.0;
 
299
 
 
300
            G_verbose_message(_("Reading raster map <%s> using weight %f..."), p->name, p->weight);
 
301
            p->buf = Rast_allocate_d_buf();
 
302
            if (!flag.lazy->answer)
 
303
                p->fd = Rast_open_old(p->name, "");
 
304
        }
220
305
    }
221
306
 
222
307
    /* process the output maps */
223
 
    for (i = 0; parm.output->answers[i]; i++) ;
 
308
    for (i = 0; parm.output->answers[i]; i++)
 
309
        ;
224
310
    num_outputs = i;
225
311
 
226
 
    for (i = 0; parm.method->answers[i]; i++) ;
 
312
    for (i = 0; parm.method->answers[i]; i++)
 
313
        ;
227
314
    if (num_outputs != i)
228
315
        G_fatal_error(_("output= and method= must have the same number of values"));
229
316
 
240
327
        out->quantile = (parm.quantile->answer && parm.quantile->answers[i])
241
328
            ? atof(parm.quantile->answers[i])
242
329
            : 0;
243
 
        out->threshold = (parm.threshold->answer && parm.threshold->answers[i])
244
 
            ? atof(parm.threshold->answers[i])
245
 
            : 0;
246
 
        out->buf = G_allocate_d_raster_buf();
247
 
        out->fd = G_open_raster_new(
248
 
            output_name, menu[method].is_int ? CELL_TYPE : DCELL_TYPE);
249
 
        if (out->fd < 0)
250
 
            G_fatal_error(_("Unable to create raster map <%s>"), out->name);
 
330
        out->buf = Rast_allocate_d_buf();
 
331
        out->fd = Rast_open_new(output_name,
 
332
                                menu[method].is_int ? CELL_TYPE : DCELL_TYPE);
251
333
    }
252
334
 
253
335
    /* initialise variables */
254
336
    values = G_malloc(num_inputs * sizeof(DCELL));
255
337
    values_tmp = G_malloc(num_inputs * sizeof(DCELL));
256
338
 
257
 
    nrows = G_window_rows();
258
 
    ncols = G_window_cols();
 
339
    nrows = Rast_window_rows();
 
340
    ncols = Rast_window_cols();
259
341
 
260
342
    /* process the data */
261
343
    G_verbose_message(_("Percent complete..."));
263
345
    for (row = 0; row < nrows; row++) {
264
346
        G_percent(row, nrows, 2);
265
347
 
266
 
        for (i = 0; i < num_inputs; i++)
267
 
            G_get_d_raster_row(inputs[i].fd, inputs[i].buf, row);
 
348
        if (flag.lazy->answer) {
 
349
            /* Open the files only on run time */
 
350
            for (i = 0; i < num_inputs; i++) {
 
351
                inputs[i].fd = Rast_open_old(inputs[i].name, "");
 
352
                Rast_get_d_row(inputs[i].fd, inputs[i].buf, row);
 
353
                Rast_close(inputs[i].fd);
 
354
            }
 
355
        }
 
356
        else {
 
357
            for (i = 0; i < num_inputs; i++)
 
358
                Rast_get_d_row(inputs[i].fd, inputs[i].buf, row);
 
359
        }
268
360
 
269
361
        for (col = 0; col < ncols; col++) {
270
362
            int null = 0;
272
364
            for (i = 0; i < num_inputs; i++) {
273
365
                DCELL v = inputs[i].buf[col];
274
366
 
275
 
                if (G_is_d_null_value(&v))
 
367
                if (Rast_is_d_null_value(&v))
276
368
                    null = 1;
277
369
                else if (parm.range->answer && (v < lo || v > hi)) {
278
 
                    G_set_d_null_value(&v, 1);
 
370
                    Rast_set_d_null_value(&v, 1);
279
371
                    null = 1;
280
372
                }
281
 
 
282
 
                values[i] = v;
 
373
                values[i] = v * inputs[i].weight;
283
374
            }
284
375
 
285
376
            for (i = 0; i < num_outputs; i++) {
286
377
                struct output *out = &outputs[i];
287
378
 
288
379
                if (null && flag.nulls->answer)
289
 
                    G_set_d_null_value(&out->buf[col], 1);
 
380
                    Rast_set_d_null_value(&out->buf[col], 1);
290
381
                else {
291
382
                    memcpy(values_tmp, values, num_inputs * sizeof(DCELL));
292
 
                    (*out->method_fn)(&out->buf[col], values_tmp, num_inputs, &out->threshold);
 
383
                    (*out->method_fn)(&out->buf[col], values_tmp, num_inputs, &out->quantile);
293
384
                }
294
385
            }
295
386
        }
296
387
 
297
388
        for (i = 0; i < num_outputs; i++)
298
 
            G_put_d_raster_row(outputs[i].fd, outputs[i].buf);
 
389
            Rast_put_d_row(outputs[i].fd, outputs[i].buf);
299
390
    }
300
391
 
301
392
    G_percent(row, nrows, 2);
302
393
 
303
 
    /* close maps */
 
394
    /* close output maps */
304
395
    for (i = 0; i < num_outputs; i++) {
305
396
        struct output *out = &outputs[i];
306
397
 
307
 
        G_close_cell(out->fd);
308
 
 
309
 
        G_short_history(out->name, "raster", &history);
310
 
        G_command_history(&history);
311
 
        G_write_history(out->name, &history);
312
 
    }
313
 
 
314
 
    for (i = 0; i < num_inputs; i++)
315
 
        G_close_cell(inputs[i].fd);
 
398
        Rast_close(out->fd);
 
399
 
 
400
        Rast_short_history(out->name, "raster", &history);
 
401
        Rast_command_history(&history);
 
402
        Rast_write_history(out->name, &history);
 
403
    }
 
404
 
 
405
    /* Close input maps */
 
406
    if (!flag.lazy->answer) {
 
407
        for (i = 0; i < num_inputs; i++)
 
408
            Rast_close(inputs[i].fd);
 
409
    }
316
410
 
317
411
    exit(EXIT_SUCCESS);
318
412
}