← Back to branch summary

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

« back to all changes in this revision

Viewing changes to raster/r.series.accumulate/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.accumulate/main.c
 
1
/****************************************************************************
 
2
 * 
 
3
 * MODULE:       r.series.accumulate
 
4
 * AUTHOR(S):    Markus Metz
 
5
 *               Soeren Gebbert
 
6
 *               based on r.series
 
7
 * PURPOSE:      Calculates (accumulated) raster value means, growing degree days
 
8
 *               (GDDs) or Winkler indices from several input maps.
 
9
 * COPYRIGHT:    (C) 2012 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
#include <limits.h>
 
20
#include <math.h>
 
21
#include <stdio.h>
 
22
 
 
23
#include <grass/gis.h>
 
24
#include <grass/raster.h>
 
25
#include <grass/glocale.h>
 
26
 
 
27
#define METHOD_GDD 1
 
28
#define METHOD_MEAN 2
 
29
#define METHOD_WINKLER 3
 
30
#define METHOD_BEDD 4
 
31
#define METHOD_HUGLIN 5
 
32
 
 
33
struct map_info
 
34
{
 
35
    const char *name;
 
36
    int fd;
 
37
    DCELL *buf;
 
38
};
 
39
 
 
40
struct map_info_out
 
41
{
 
42
    const char *name;
 
43
    int fd;
 
44
    void *buf;
 
45
};
 
46
 
 
47
int main(int argc, char *argv[])
 
48
{
 
49
    struct GModule *module;
 
50
    struct
 
51
    {
 
52
        struct Option *input, *basemap, *file, *output,
 
53
        *range, *scale, *shift, *lower,
 
54
        *upper, *limits, *method;
 
55
    } parm;
 
56
    struct
 
57
    {
 
58
        struct Flag *nulls, *lazy, *float_output;
 
59
    } flag;
 
60
    int i;
 
61
    int num_inputs, max_inputs;
 
62
    int method;
 
63
    struct map_info *inputs = NULL;
 
64
    struct map_info_out *out = NULL;
 
65
    struct map_info *basemap = NULL;
 
66
    struct map_info *map_lower = NULL;
 
67
    struct map_info *map_upper = NULL;
 
68
    struct History history;
 
69
    int nrows, ncols;
 
70
    int row, col;
 
71
    DCELL lo, hi, tscale, tshift, lower = 10.0, upper = 30.0;
 
72
    DCELL dcell_null;
 
73
    RASTER_MAP_TYPE out_type;
 
74
    int out_size;
 
75
    char *desc = NULL;
 
76
    
 
77
    G_gisinit(argv[0]);
 
78
    
 
79
    module = G_define_module();
 
80
    G_add_keyword(_("raster"));
 
81
    G_add_keyword(_("series"));
 
82
    G_add_keyword(_("accumulation"));
 
83
    module->description =
 
84
    _("Makes each output cell value a accumulation"
 
85
      "function of the values assigned to the corresponding cells "
 
86
      "in the input raster map layers.");
 
87
 
 
88
    parm.basemap = G_define_standard_option(G_OPT_R_INPUT);
 
89
    parm.basemap->key = "basemap";
 
90
    parm.basemap->description = _("Existing map to be added to output");
 
91
    parm.basemap->required = NO;
 
92
 
 
93
    parm.input = G_define_standard_option(G_OPT_R_INPUTS);
 
94
    parm.input->required = NO;
 
95
 
 
96
    parm.file = G_define_standard_option(G_OPT_F_INPUT);
 
97
    parm.file->key = "file";
 
98
    parm.file->description = _("Input file with raster map names, one per line");
 
99
    parm.file->required = NO;
 
100
    
 
101
    parm.output = G_define_standard_option(G_OPT_R_OUTPUT);
 
102
    parm.output->multiple = NO;
 
103
    
 
104
    parm.scale = G_define_option();
 
105
    parm.scale->key = "scale";
 
106
    parm.scale->type = TYPE_DOUBLE;
 
107
    parm.scale->answer = "1.0";
 
108
    parm.scale->required = NO;
 
109
    parm.scale->description = _("Scale factor for input");
 
110
    
 
111
    parm.shift = G_define_option();
 
112
    parm.shift->key = "shift";
 
113
    parm.shift->type = TYPE_DOUBLE;
 
114
    parm.shift->answer = "0.0";
 
115
    parm.shift->required = NO;
 
116
    parm.shift->description = _("Shift factor for input");
 
117
    
 
118
    parm.lower = G_define_standard_option(G_OPT_R_INPUT);
 
119
    parm.lower->key = "lower";
 
120
    parm.lower->required = NO;
 
121
    parm.lower->description = _("The raster map specifying the lower accumulation limit, also called baseline");
 
122
    
 
123
    parm.upper = G_define_standard_option(G_OPT_R_INPUT);
 
124
    parm.upper->key = "upper";
 
125
    parm.upper->required = NO;
 
126
    parm.upper->description = _("The raster map specifying the upper accumulation limit, also called cutoff. Only applied to BEDD computation.");
 
127
    
 
128
    parm.range = G_define_option();
 
129
    parm.range->key = "range";
 
130
    parm.range->type = TYPE_DOUBLE;
 
131
    parm.range->key_desc = "min,max";
 
132
    parm.range->description = _("Ignore values outside this range");
 
133
    
 
134
    parm.limits = G_define_option();
 
135
    parm.limits->key = "limits";
 
136
    parm.limits->type = TYPE_DOUBLE;
 
137
    parm.limits->key_desc = "lower,upper";
 
138
    parm.limits->answer = "10,30";
 
139
    parm.limits->description = _("Use these limits in case lower and/or upper input maps are not defined");
 
140
    
 
141
    parm.method = G_define_option();
 
142
    parm.method->key = "method";
 
143
    parm.method->type = TYPE_STRING;
 
144
    parm.method->multiple = NO;
 
145
    parm.method->required = NO;
 
146
    parm.method->options = "gdd,bedd,huglin,mean";
 
147
    parm.method->answer = "gdd";
 
148
    parm.method->label = "This method will be applied to compute the accumulative values from the input maps";
 
149
    G_asprintf(&desc,
 
150
           "gdd;%s;mean;%s;bedd;%s;huglin;%s",
 
151
           _("Growing Degree Days or Winkler indices"),
 
152
           _("Mean: sum(input maps)/(number of input maps)"),
 
153
           _("Biologically Effective Degree Days"),
 
154
           _("Huglin Heliothermal index"));
 
155
    parm.method->descriptions = desc;
 
156
 
 
157
    flag.nulls = G_define_flag();
 
158
    flag.nulls->key = 'n';
 
159
    flag.nulls->description = _("Propagate NULLs");
 
160
    
 
161
    flag.lazy = G_define_flag();
 
162
    flag.lazy->key = 'z';
 
163
    flag.lazy->description = _("Do not keep files open");
 
164
    
 
165
    flag.float_output = G_define_flag();
 
166
    flag.float_output->key = 'f';
 
167
    flag.float_output->description = _("Create a FCELL map (floating point single precision) as output");
 
168
    
 
169
    if (G_parser(argc, argv))
 
170
        exit(EXIT_FAILURE);
 
171
    
 
172
    lo = -1.0 / 0.0;    /* -inf */
 
173
    hi =  1.0 / 0.0;    /* inf */
 
174
    
 
175
    method = METHOD_GDD;
 
176
    if (G_strncasecmp(parm.method->answer, "gdd", 3) == 0)
 
177
        method = METHOD_GDD;
 
178
    else if (G_strncasecmp(parm.method->answer, "mean", 4) == 0)
 
179
        method = METHOD_MEAN;
 
180
    else if (G_strncasecmp(parm.method->answer, "bedd", 4) == 0)
 
181
        method = METHOD_BEDD;
 
182
    else if (G_strncasecmp(parm.method->answer, "huglin", 7) == 0)
 
183
        method = METHOD_HUGLIN;
 
184
 
 
185
    if (parm.range->answer) {
 
186
        lo = atof(parm.range->answers[0]);
 
187
        hi = atof(parm.range->answers[1]);
 
188
    }
 
189
    
 
190
    if (parm.limits->answer) {
 
191
        lower = atof(parm.limits->answers[0]);
 
192
        upper = atof(parm.limits->answers[1]);
 
193
    }
 
194
    
 
195
    if (parm.scale->answer)
 
196
        tscale = atof(parm.scale->answer);
 
197
    else
 
198
        tscale = 1.;
 
199
    
 
200
    if (parm.shift->answer)
 
201
        tshift = atof(parm.shift->answer);
 
202
    else
 
203
        tshift = 0.;
 
204
    
 
205
    if (parm.input->answer && parm.file->answer)
 
206
        G_fatal_error(_("%s= and %s= are mutually exclusive"),
 
207
                        parm.input->key, parm.file->key);
 
208
    
 
209
    if (!parm.input->answer && !parm.file->answer)
 
210
        G_fatal_error(_("Please specify %s= or %s="),
 
211
                        parm.input->key, parm.file->key);
 
212
    
 
213
    max_inputs = 0;
 
214
    
 
215
    /* process the input maps from the file */
 
216
    if (parm.file->answer) {
 
217
        FILE *in;
 
218
        
 
219
        in = fopen(parm.file->answer, "r");
 
220
        if (!in)
 
221
            G_fatal_error(_("Unable to open input file <%s>"),
 
222
                          parm.file->answer);
 
223
            
 
224
            num_inputs = 0;
 
225
        
 
226
        for (;;) {
 
227
            char buf[GNAME_MAX];
 
228
            char *name;
 
229
            struct map_info *p;
 
230
            
 
231
            if (!G_getl2(buf, sizeof(buf), in))
 
232
                break;
 
233
            
 
234
            name = G_chop(buf);
 
235
            
 
236
            /* Ignore empty lines */
 
237
            if (!*name)
 
238
                continue;
 
239
            
 
240
            if (num_inputs >= max_inputs) {
 
241
                max_inputs += 100;
 
242
                inputs =
 
243
                G_realloc(inputs, max_inputs * sizeof(struct map_info));
 
244
            }
 
245
            p = &inputs[num_inputs++];
 
246
            
 
247
            p->name = G_store(name);
 
248
            G_verbose_message(_("Reading raster map <%s>..."), p->name);
 
249
            p->buf = Rast_allocate_d_buf();
 
250
            if (!flag.lazy->answer)
 
251
                p->fd = Rast_open_old(p->name, "");
 
252
        }
 
253
        
 
254
        if (num_inputs < 1)
 
255
            G_fatal_error(_("No raster map name found in input file"));
 
256
        
 
257
        fclose(in);
 
258
    }
 
259
    else {
 
260
        for (i = 0; parm.input->answers[i]; i++)
 
261
            ;
 
262
        num_inputs = i;
 
263
        
 
264
        if (num_inputs < 1)
 
265
            G_fatal_error(_("Raster map not found"));
 
266
        
 
267
        inputs = G_malloc(num_inputs * sizeof(struct map_info));
 
268
        
 
269
        for (i = 0; i < num_inputs; i++) {
 
270
            struct map_info *p = &inputs[i];
 
271
            
 
272
            p->name = parm.input->answers[i];
 
273
            G_verbose_message(_("Reading raster map <%s>..."), p->name);
 
274
            p->buf = Rast_allocate_d_buf();
 
275
            if (!flag.lazy->answer)
 
276
                p->fd = Rast_open_old(p->name, "");
 
277
        }
 
278
        max_inputs = num_inputs;
 
279
    }
 
280
    
 
281
    if (parm.basemap->answer) {
 
282
        basemap = G_malloc(1 * sizeof(struct map_info));
 
283
        basemap->name = parm.basemap->answer;
 
284
        G_verbose_message(_("Reading raster map <%s>..."), basemap->name);
 
285
        basemap->buf = Rast_allocate_d_buf();
 
286
        basemap->fd = Rast_open_old(basemap->name, "");
 
287
    }
 
288
    
 
289
    if (parm.lower->answer) {
 
290
        map_lower = G_malloc(1 * sizeof(struct map_info));
 
291
        map_lower->name = parm.lower->answer;
 
292
        G_verbose_message(_("Reading raster map <%s>..."), map_lower->name);
 
293
        map_lower->buf = Rast_allocate_d_buf();
 
294
        map_lower->fd = Rast_open_old(map_lower->name, "");
 
295
    }
 
296
    
 
297
    if (parm.upper->answer) {
 
298
        map_upper = G_malloc(1 * sizeof(struct map_info));
 
299
        map_upper->name = parm.upper->answer;
 
300
        G_verbose_message(_("Reading raster map <%s>..."), map_upper->name);
 
301
        map_upper->buf = Rast_allocate_d_buf();
 
302
        map_upper->fd = Rast_open_old(map_upper->name, "");
 
303
    }
 
304
    
 
305
    out = G_calloc(1, sizeof(struct map_info_out));
 
306
    out->name = parm.output->answer;
 
307
    out_type = flag.float_output->answer ? FCELL_TYPE : DCELL_TYPE;
 
308
    out->buf = Rast_allocate_buf(out_type);
 
309
    out_size = Rast_cell_size(out_type);
 
310
    out->fd = Rast_open_new(out->name, out_type);
 
311
    
 
312
    nrows = Rast_window_rows();
 
313
    ncols = Rast_window_cols();
 
314
    
 
315
    Rast_set_d_null_value(&dcell_null, 1);
 
316
    
 
317
    /* process the data */
 
318
    G_verbose_message(_("Percent complete..."));
 
319
    
 
320
    for (row = 0; row < nrows; row++) {
 
321
        G_percent(row, nrows, 2);
 
322
        
 
323
        if (basemap)
 
324
            Rast_get_d_row(basemap->fd, basemap->buf, row);
 
325
        if (map_lower)
 
326
            Rast_get_d_row(map_lower->fd, map_lower->buf, row);
 
327
        if (map_upper)
 
328
            Rast_get_d_row(map_upper->fd, map_upper->buf, row);
 
329
        
 
330
        if (flag.lazy->answer) {
 
331
            /* Open the files only on run time */
 
332
            for (i = 0; i < num_inputs; i++) {
 
333
                inputs[i].fd = Rast_open_old(inputs[i].name, "");
 
334
                Rast_get_d_row(inputs[i].fd, inputs[i].buf, row);
 
335
                Rast_close(inputs[i].fd);
 
336
            }
 
337
        }
 
338
        else {
 
339
            for (i = 0; i < num_inputs; i++)
 
340
                Rast_get_d_row(inputs[i].fd, inputs[i].buf, row);
 
341
        }
 
342
        
 
343
        #pragma omp for schedule (static) private (col)
 
344
        for (col = 0; col < ncols; col++) {
 
345
            int null = 0, non_null = 0;
 
346
            DCELL min, max, avg, value;
 
347
            
 
348
            if (map_lower)
 
349
                lower = map_lower->buf[col];
 
350
            if (map_upper)
 
351
                upper = map_upper->buf[col];
 
352
            
 
353
            if (upper <= lower)
 
354
                G_fatal_error(_("'%s'=%f must be > '%s'=%f"), parm.upper->key, upper,
 
355
                              parm.lower->key, lower);
 
356
                
 
357
            min = dcell_null;
 
358
            max = dcell_null;
 
359
            avg = 0;
 
360
            
 
361
            for (i = 0; i < num_inputs; i++) {
 
362
                DCELL v = inputs[i].buf[col];
 
363
                
 
364
                if (Rast_is_d_null_value(&v))
 
365
                    null = 1;
 
366
                else {
 
367
                    v = v * tscale + tshift;
 
368
                    if (parm.range->answer && (v < lo || v > hi)) {
 
369
                        null = 1;
 
370
                    }
 
371
                    else  {
 
372
                        avg += v;
 
373
                        if (min > v || Rast_is_d_null_value(&min))
 
374
                            min = v;
 
375
                        if (max < v || Rast_is_d_null_value(&max))
 
376
                            max = v;
 
377
                        non_null++;
 
378
                    }
 
379
                }
 
380
            }
 
381
 
 
382
            value = dcell_null;
 
383
            if (!non_null || (null && flag.nulls->answer)) {
 
384
                if (basemap)
 
385
                    value = basemap->buf[col];
 
386
            }
 
387
            else {
 
388
                /* Compute mean or index */
 
389
                avg /= non_null;
 
390
 
 
391
                switch(method) {
 
392
                    case METHOD_HUGLIN:
 
393
                        avg = (avg + max) / 2;
 
394
                        break;
 
395
                    case METHOD_BEDD:
 
396
                        if(avg > upper)
 
397
                            avg = upper;
 
398
                        break;
 
399
                    case METHOD_MEAN:
 
400
                        value = avg;
 
401
                        break;
 
402
                    default:
 
403
                        /* Winkler or GDD index computation is the default */
 
404
                        break;
 
405
                }
 
406
                if (method != METHOD_MEAN) {
 
407
                    value = avg - lower;
 
408
 
 
409
                    if (value < 0.)
 
410
                        value = 0.;
 
411
                }
 
412
 
 
413
                if (basemap)
 
414
                    value += basemap->buf[col];
 
415
            }
 
416
            Rast_set_d_value((char *)out->buf + col * out_size, value, out_type);
 
417
        }
 
418
        Rast_put_row(out->fd, out->buf, out_type);
 
419
    }
 
420
    
 
421
    G_percent(row, nrows, 2);
 
422
    
 
423
    /* close output map */
 
424
    Rast_close(out->fd);
 
425
    
 
426
    Rast_short_history(out->name, "raster", &history);
 
427
    Rast_command_history(&history);
 
428
    Rast_write_history(out->name, &history);
 
429
    
 
430
    /* Close input maps */
 
431
    if (basemap)
 
432
        Rast_close(basemap->fd);
 
433
    if (map_lower)
 
434
        Rast_close(map_lower->fd);
 
435
    if (map_upper)
 
436
        Rast_close(map_upper->fd);
 
437
    
 
438
    if (!flag.lazy->answer) {
 
439
        for (i = 0; i < num_inputs; i++)
 
440
            Rast_close(inputs[i].fd);
 
441
    }
 
442
    
 
443
    if (method == METHOD_GDD) {
 
444
        struct Colors colors;
 
445
        
 
446
        Rast_init_colors(&colors);
 
447
        Rast_make_colors(&colors, "gdd", 0, 6000);
 
448
        Rast_write_colors(out->name, G_mapset(), &colors);
 
449
    }
 
450
    
 
451
    exit(EXIT_SUCCESS);
 
452
}