← Back to branch summary

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

« back to all changes in this revision

Viewing changes to raster/r.univar2/r.univar_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.univar2/r.univar_main.c
1
 
/*
2
 
 * r.univar
3
 
 *
4
 
 *  Calculates univariate statistics from the non-null cells of a GRASS raster map
5
 
 *
6
 
 *   Copyright (C) 2004-2006 by the GRASS Development Team
7
 
 *   Author(s): Hamish Bowman, University of Otago, New Zealand
8
 
 *              Extended stats: Martin Landa
9
 
 *              Zonal stats: Markus Metz
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
 
 *   This program is a replacement for the r.univar shell script
16
 
 */
17
 
 
18
 
#include <assert.h>
19
 
#include <string.h>
20
 
#define MAIN
21
 
#include "globals.h"
22
 
 
23
 
/* local proto */
24
 
void set_params();
25
 
 
26
 
/* ************************************************************************* */
27
 
/* Set up the arguments we are expecting ********************************** */
28
 
/* ************************************************************************* */
29
 
void set_params()
30
 
{
31
 
    param.inputfile = G_define_standard_option(G_OPT_R_MAPS);
32
 
 
33
 
    param.zonefile = G_define_standard_option(G_OPT_R_MAP);
34
 
    param.zonefile->key = "zones";
35
 
    param.zonefile->required = NO;
36
 
    param.zonefile->description =
37
 
        _("Raster map used for zoning, must be of type CELL");
38
 
 
39
 
    param.output_file = G_define_standard_option(G_OPT_F_OUTPUT);
40
 
    param.output_file->required = NO;
41
 
    param.output_file->description =
42
 
        _("Name for output file (if omitted or \"-\" output to stdout)");
43
 
 
44
 
    param.percentile = G_define_option();
45
 
    param.percentile->key = "percentile";
46
 
    param.percentile->type = TYPE_DOUBLE;
47
 
    param.percentile->required = NO;
48
 
    param.percentile->multiple = YES;
49
 
    param.percentile->options = "0-100";
50
 
    param.percentile->answer = "90";
51
 
    param.percentile->description =
52
 
        _("Percentile to calculate (requires extended statistics flag)");
53
 
 
54
 
    param.separator = G_define_standard_option(G_OPT_F_SEP);
55
 
    param.separator->description = _("Special characters: space, comma, tab");
56
 
 
57
 
    param.shell_style = G_define_flag();
58
 
    param.shell_style->key = 'g';
59
 
    param.shell_style->description =
60
 
        _("Print the stats in shell script style");
61
 
 
62
 
    param.extended = G_define_flag();
63
 
    param.extended->key = 'e';
64
 
    param.extended->description = _("Calculate extended statistics");
65
 
 
66
 
    param.table = G_define_flag();
67
 
    param.table->key = 't';
68
 
    param.table->description = _("Table output format instead of standard output format");
69
 
 
70
 
    return;
71
 
}
72
 
 
73
 
static int open_raster(const char *infile);
74
 
static univar_stat *univar_stat_with_percentiles(int map_type);
75
 
static void process_raster(univar_stat * stats, int fd, int fdz,
76
 
                           const struct Cell_head *region);
77
 
 
78
 
/* *************************************************************** */
79
 
/* **** the main functions for r.univar ************************** */
80
 
/* *************************************************************** */
81
 
int main(int argc, char *argv[])
82
 
{
83
 
    unsigned int rows, cols;    /*  totals  */
84
 
    int rasters;
85
 
 
86
 
    struct Cell_head region;
87
 
    struct GModule *module;
88
 
    univar_stat *stats;
89
 
    char **p, *z;
90
 
    int fd, fdz, cell_type, min, max;
91
 
    struct Range zone_range;
92
 
    char *mapset, *name;
93
 
 
94
 
    G_gisinit(argv[0]);
95
 
 
96
 
    module = G_define_module();
97
 
    module->keywords = _("raster, statistics");
98
 
    module->description =
99
 
        _("Calculates univariate statistics from the non-null cells of a raster map.");
100
 
 
101
 
    /* Define the different options */
102
 
    set_params();
103
 
 
104
 
    if (G_parser(argc, argv))
105
 
        exit(EXIT_FAILURE);
106
 
 
107
 
    name = param.output_file->answer;
108
 
    if (name != NULL && strcmp(name, "-") != 0) {
109
 
        if (NULL == freopen(name, "w", stdout)) {
110
 
            G_fatal_error(_("Unable to open file <%s> for writing"), name);
111
 
        }
112
 
    }
113
 
 
114
 
    G_get_window(&region);
115
 
    rows = region.rows;
116
 
    cols = region.cols;
117
 
 
118
 
    /* table field separator */
119
 
    zone_info.sep = param.separator->answer;
120
 
    if (strcmp(zone_info.sep, "\\t") == 0)
121
 
        zone_info.sep = "\t";
122
 
    if (strcmp(zone_info.sep, "tab") == 0)
123
 
        zone_info.sep = "\t";
124
 
    if (strcmp(zone_info.sep, "space") == 0)
125
 
        zone_info.sep = " ";
126
 
    if (strcmp(zone_info.sep, "comma") == 0)
127
 
        zone_info.sep = ",";
128
 
 
129
 
    zone_info.min = 0.0 / 0.0;  /* set to nan as default */
130
 
    zone_info.max = 0.0 / 0.0;  /* set to nan as default */
131
 
    zone_info.n_zones = 0;
132
 
 
133
 
    fdz = -1;
134
 
    
135
 
    /* open zoning raster */
136
 
    if ((z = param.zonefile->answer)) {
137
 
        mapset = G_find_cell2(z, "");
138
 
 
139
 
        fdz = open_raster(z);
140
 
        
141
 
        cell_type = G_get_raster_map_type(fdz);
142
 
        if (cell_type != CELL_TYPE)
143
 
            G_fatal_error("Zoning raster must be of type CELL");
144
 
 
145
 
        if (G_read_range(z, mapset, &zone_range) == -1)
146
 
            G_fatal_error("Can not read range for zoning raster");
147
 
        if (G_get_range_min_max(&zone_range, &min, &max))
148
 
            G_fatal_error("Can not read range for zoning raster");
149
 
        if (G_read_raster_cats(z, mapset, &(zone_info.cats)))
150
 
            G_warning("no category support for zoning raster");
151
 
 
152
 
        zone_info.min = min;
153
 
        zone_info.max = max;
154
 
        zone_info.n_zones = max - min + 1;
155
 
    }
156
 
 
157
 
    /* count the input rasters given */
158
 
    for (p = (char **)param.inputfile->answers, rasters = 0;
159
 
         *p; p++, rasters++) ;
160
 
 
161
 
    /* process all input rasters */
162
 
    int map_type = param.extended->answer ? -2 : -1;
163
 
 
164
 
    stats = ((map_type == -1)
165
 
             ? create_univar_stat_struct(-1, 0)
166
 
             : 0);
167
 
 
168
 
    for (p = param.inputfile->answers; *p; p++) {
169
 
        fd = open_raster(*p);
170
 
 
171
 
        if (map_type != -1) {
172
 
            /* NB: map_type must match when doing extended stats */
173
 
            int this_type = G_get_raster_map_type(fd);
174
 
 
175
 
            assert(this_type > -1);
176
 
            if (map_type < -1) {
177
 
                /* extended stats */
178
 
                assert(stats == 0);
179
 
                map_type = this_type;
180
 
                stats = univar_stat_with_percentiles(map_type);
181
 
            }
182
 
            else if (this_type != map_type) {
183
 
                G_fatal_error(_("Raster <%s> type mismatch"), *p);
184
 
            }
185
 
        }
186
 
 
187
 
        process_raster(stats, fd, fdz, &region);
188
 
 
189
 
        /* close input raster */
190
 
        G_close_cell(fd);
191
 
    }
192
 
 
193
 
    /* close zoning raster */
194
 
    if (z)
195
 
        G_close_cell(fdz);
196
 
 
197
 
    /* create the output */
198
 
    if (param.table->answer)
199
 
        print_stats_table(stats);
200
 
    else
201
 
        print_stats(stats);
202
 
        
203
 
    /* release memory */
204
 
    free_univar_stat_struct(stats);
205
 
 
206
 
    exit(EXIT_SUCCESS);
207
 
}
208
 
 
209
 
static int open_raster(const char *infile)
210
 
{
211
 
    char *mapset;
212
 
    int fd;
213
 
 
214
 
    mapset = G_find_cell2(infile, "");
215
 
    if (mapset == NULL) {
216
 
        G_fatal_error(_("Raster map <%s> not found"), infile);
217
 
    }
218
 
 
219
 
    fd = G_open_cell_old(infile, mapset);
220
 
    if (fd < 0)
221
 
        G_fatal_error(_("Unable to open raster map <%s>"), infile);
222
 
 
223
 
    /* . */
224
 
    return fd;
225
 
}
226
 
 
227
 
static univar_stat *univar_stat_with_percentiles(int map_type)
228
 
{
229
 
    univar_stat *stats;
230
 
    int i, j;
231
 
    int n_zones = zone_info.n_zones;
232
 
 
233
 
    if (n_zones == 0)
234
 
        n_zones = 1;
235
 
 
236
 
    i = 0;
237
 
    while (param.percentile->answers[i])
238
 
        i++;
239
 
    stats = create_univar_stat_struct(map_type, i);
240
 
    for (i = 0; i < n_zones; i++) {
241
 
        for (j = 0; j < stats[i].n_perc; j++) {
242
 
            sscanf(param.percentile->answers[j], "%lf", &(stats[i].perc[j]));
243
 
        }
244
 
    }
245
 
 
246
 
    /* . */
247
 
    return stats;
248
 
}
249
 
 
250
 
static void
251
 
process_raster(univar_stat * stats, int fd, int fdz, const struct Cell_head *region)
252
 
{
253
 
    /* use G_window_rows(), G_window_cols() here? */
254
 
    const int rows = region->rows;
255
 
    const int cols = region->cols;
256
 
 
257
 
    const RASTER_MAP_TYPE map_type = G_get_raster_map_type(fd);
258
 
    const size_t value_sz = G_raster_size(map_type);
259
 
    unsigned int row;
260
 
    void *raster_row;
261
 
    CELL *zoneraster_row;
262
 
    int n_zones = zone_info.n_zones;
263
 
    
264
 
    raster_row = G_allocate_raster_buf(map_type);
265
 
    if (n_zones)
266
 
        zoneraster_row = G_allocate_c_raster_buf();
267
 
 
268
 
    for (row = 0; row < rows; row++) {
269
 
        void *ptr;
270
 
        CELL *zptr;
271
 
        unsigned int col;
272
 
 
273
 
        if (G_get_raster_row(fd, raster_row, row, map_type) < 0)
274
 
            G_fatal_error(_("Reading row %d"), row);
275
 
        if (n_zones) {
276
 
            if (G_get_c_raster_row(fdz, zoneraster_row, row) < 0)
277
 
                G_fatal_error(_("Reading row %d"), row);
278
 
            zptr = zoneraster_row;
279
 
        }
280
 
 
281
 
        ptr = raster_row;
282
 
 
283
 
        for (col = 0; col < cols; col++) {
284
 
            double val;
285
 
            int zone = 0;
286
 
 
287
 
            if (n_zones) {
288
 
                /* skip NULL cells in zone map */
289
 
                if (G_is_c_null_value(zptr)) {
290
 
                    ptr = G_incr_void_ptr(ptr, value_sz);
291
 
                    zptr++;
292
 
                    continue;
293
 
                }
294
 
                zone = *zptr - zone_info.min;
295
 
            }
296
 
 
297
 
            /* count NULL cells in input map */
298
 
            stats[zone].size++;
299
 
            
300
 
            /* can't do stats with NULL cells in input map */
301
 
            if (G_is_null_value(ptr, map_type)) {
302
 
                ptr = G_incr_void_ptr(ptr, value_sz);
303
 
                if (n_zones)
304
 
                    zptr++;
305
 
                continue;
306
 
            }
307
 
 
308
 
            if (param.extended->answer) {
309
 
                /* check allocated memory */
310
 
                if (stats[zone].n >= stats[zone].n_alloc) {
311
 
                    stats[zone].n_alloc += 1000;
312
 
                    size_t msize;
313
 
                    switch (map_type) {
314
 
                        case DCELL_TYPE:
315
 
                            msize = stats[zone].n_alloc * sizeof(DCELL);
316
 
                            stats[zone].dcell_array =
317
 
                                (DCELL *)G_realloc((void *)stats[zone].dcell_array, msize);
318
 
                            stats[zone].nextp = (void *)&(stats[zone].dcell_array[stats[zone].n]);
319
 
                            break;
320
 
                        case FCELL_TYPE:
321
 
                            msize = stats[zone].n_alloc * sizeof(FCELL);
322
 
                            stats[zone].fcell_array =
323
 
                                (FCELL *)G_realloc((void *)stats[zone].fcell_array, msize);
324
 
                            stats[zone].nextp = (void *)&(stats[zone].fcell_array[stats[zone].n]);
325
 
                            break;
326
 
                        case CELL_TYPE:
327
 
                            msize = stats[zone].n_alloc * sizeof(CELL);
328
 
                            stats[zone].cell_array =
329
 
                                (CELL *)G_realloc((void *)stats[zone].cell_array, msize);
330
 
                            stats[zone].nextp = (void *)&(stats[zone].cell_array[stats[zone].n]);
331
 
                            break;
332
 
                        default:
333
 
                            break;
334
 
                    }
335
 
                }
336
 
                /* put the value into stats->XXXcell_array */
337
 
                memcpy(stats[zone].nextp, ptr, value_sz);
338
 
                stats[zone].nextp = G_incr_void_ptr(stats[zone].nextp, value_sz);
339
 
            }
340
 
 
341
 
            val = ((map_type == DCELL_TYPE) ? *((DCELL *) ptr)
342
 
                          : (map_type == FCELL_TYPE) ? *((FCELL *) ptr)
343
 
                          : *((CELL *) ptr));
344
 
 
345
 
            stats[zone].sum += val;
346
 
            stats[zone].sumsq += val * val;
347
 
            stats[zone].sum_abs += fabs(val);
348
 
 
349
 
            if (stats[zone].first) {
350
 
                stats[zone].max = val;
351
 
                stats[zone].min = val;
352
 
                stats[zone].first = FALSE;
353
 
            }
354
 
            else {
355
 
                if (val > stats[zone].max)
356
 
                    stats[zone].max = val;
357
 
                if (val < stats[zone].min)
358
 
                    stats[zone].min = val;
359
 
            }
360
 
 
361
 
            ptr = G_incr_void_ptr(ptr, value_sz);
362
 
            if (n_zones)
363
 
                zptr++;
364
 
            stats[zone].n++;
365
 
        }
366
 
        if (!(param.shell_style->answer))
367
 
            G_percent(row, rows, 2);
368
 
    }
369
 
    if (!(param.shell_style->answer))
370
 
        G_percent(rows, rows, 2);       /* finish it off */
371
 
 
372
 
}