← Back to branch summary

~grass/grass/releasebranch_7_0

« back to all changes in this revision

Viewing changes to raster/r.univar2/r.univar_main.c

  • Committer: glynn
  • Date: 2008-08-16 11:51:17 UTC
  • Revision ID: svn-v4:15284696-431f-4ddb-bdfa-cd5b030d7da7:grass/trunk:32813
Rename directories r.univar2 -> r.univar, r.grow2 -> r.grow, r.cats -> r.category
Remove r.proj

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
 
 *
10
 
 *      This program is free software under the GNU General Public
11
 
 *      License (>=v2). Read the file COPYING that comes with GRASS
12
 
 *      for details.
13
 
 *
14
 
 *   This program is a replacement for the r.univar shell script
15
 
 */
16
 
 
17
 
#include <assert.h>
18
 
#include <string.h>
19
 
#include "globals.h"
20
 
 
21
 
param_type param;
22
 
 
23
 
/* ************************************************************************* */
24
 
/* Set up the arguments we are expecting ********************************** */
25
 
/* ************************************************************************* */
26
 
void set_params(void)
27
 
{
28
 
    param.inputfile = G_define_standard_option(G_OPT_R_MAPS);
29
 
 
30
 
    param.percentile = G_define_option();
31
 
    param.percentile->key = "percentile";
32
 
    param.percentile->type = TYPE_INTEGER;
33
 
    param.percentile->required = NO;
34
 
    param.percentile->multiple = YES;
35
 
    param.percentile->options = "0-100";
36
 
    param.percentile->answer = "90";
37
 
    param.percentile->description =
38
 
        _("Percentile to calculate (requires extended statistics flag)");
39
 
 
40
 
    param.shell_style = G_define_flag();
41
 
    param.shell_style->key = 'g';
42
 
    param.shell_style->description =
43
 
        _("Print the stats in shell script style");
44
 
 
45
 
    param.extended = G_define_flag();
46
 
    param.extended->key = 'e';
47
 
    param.extended->description = _("Calculate extended statistics");
48
 
 
49
 
    return;
50
 
}
51
 
 
52
 
static int open_raster(const char *infile);
53
 
static univar_stat *univar_stat_with_percentiles(int map_type, int size);
54
 
static void process_raster(univar_stat * stats, int fd,
55
 
                           const struct Cell_head *region);
56
 
 
57
 
/* *************************************************************** */
58
 
/* **** the main functions for r.univar ************************** */
59
 
/* *************************************************************** */
60
 
int main(int argc, char *argv[])
61
 
{
62
 
    unsigned int rows, cols;    /*  totals  */
63
 
    int rasters;
64
 
 
65
 
    struct Cell_head region;
66
 
    struct GModule *module;
67
 
    univar_stat *stats;
68
 
 
69
 
 
70
 
    G_gisinit(argv[0]);
71
 
 
72
 
    module = G_define_module();
73
 
    module->keywords = _("raster, statistics");
74
 
    module->description =
75
 
        _("Calculates univariate statistics from the non-null cells of a raster map.");
76
 
 
77
 
    /* Define the different options */
78
 
    set_params();
79
 
 
80
 
    if (G_parser(argc, argv))
81
 
        exit(EXIT_FAILURE);
82
 
 
83
 
    G_get_window(&region);
84
 
    rows = region.rows;
85
 
    cols = region.cols;
86
 
 
87
 
    /* count the rasters given */
88
 
    {
89
 
        const char **p;
90
 
 
91
 
        for (p = (const char **)param.inputfile->answers, rasters = 0;
92
 
             *p; p++, rasters++) ;
93
 
    }
94
 
 
95
 
    /* process them all */
96
 
    {
97
 
        size_t cells = rasters * cols * rows;
98
 
        int map_type = param.extended->answer ? -2 : -1;
99
 
        char **p;
100
 
 
101
 
        stats = ((map_type == -1)
102
 
                 ? create_univar_stat_struct(-1, cells, 0)
103
 
                 : 0);
104
 
 
105
 
        for (p = param.inputfile->answers; *p; p++) {
106
 
            int fd = open_raster(*p);
107
 
 
108
 
            if (map_type != -1) {
109
 
                /* NB: map_type must match when doing extended stats */
110
 
                int this_type = G_get_raster_map_type(fd);
111
 
 
112
 
                assert(this_type > -1);
113
 
                if (map_type < -1) {
114
 
                    assert(stats == 0);
115
 
                    map_type = this_type;
116
 
                    stats = univar_stat_with_percentiles(map_type, cells);
117
 
                }
118
 
                else if (this_type != map_type) {
119
 
                    G_fatal_error(_("Raster <%s> type mismatch"), *p);
120
 
                }
121
 
            }
122
 
 
123
 
            process_raster(stats, fd, &region);
124
 
        }
125
 
    }
126
 
 
127
 
    if (!(param.shell_style->answer))
128
 
        G_percent(rows, rows, 2);       /* finish it off */
129
 
 
130
 
    /* create the output */
131
 
    print_stats(stats);
132
 
 
133
 
    /* release memory */
134
 
    free_univar_stat_struct(stats);
135
 
 
136
 
    exit(EXIT_SUCCESS);
137
 
}
138
 
 
139
 
static int open_raster(const char *infile)
140
 
{
141
 
    char *mapset;
142
 
    int fd;
143
 
 
144
 
    mapset = G_find_cell2(infile, "");
145
 
    if (mapset == NULL) {
146
 
        G_fatal_error(_("Raster map <%s> not found"), infile);
147
 
    }
148
 
 
149
 
    fd = G_open_cell_old(infile, mapset);
150
 
    if (fd < 0)
151
 
        G_fatal_error(_("Unable to open raster map <%s>"), infile);
152
 
 
153
 
    /* . */
154
 
    return fd;
155
 
}
156
 
 
157
 
static univar_stat *univar_stat_with_percentiles(int map_type, int size)
158
 
{
159
 
    univar_stat *stats;
160
 
    int i;
161
 
 
162
 
    i = 0;
163
 
    while (param.percentile->answers[i])
164
 
        i++;
165
 
    stats = create_univar_stat_struct(map_type, size, i);
166
 
    for (i = 0; i < stats->n_perc; i++) {
167
 
        sscanf(param.percentile->answers[i], "%i", &stats->perc[i]);
168
 
    }
169
 
 
170
 
    /* . */
171
 
    return stats;
172
 
}
173
 
 
174
 
static void
175
 
process_raster(univar_stat * stats, int fd, const struct Cell_head *region)
176
 
{
177
 
    /* use G_window_rows(), G_window_cols() here? */
178
 
    const int rows = region->rows;
179
 
    const int cols = region->cols;
180
 
    int first = (stats->n < 1);
181
 
 
182
 
    const RASTER_MAP_TYPE map_type = G_get_raster_map_type(fd);
183
 
    void *nextp
184
 
        = ((!param.extended->answer) ? 0
185
 
           : (map_type == DCELL_TYPE) ? (void *)stats->dcell_array
186
 
           : (map_type == FCELL_TYPE) ? (void *)stats->fcell_array
187
 
           : (void *)stats->cell_array);
188
 
    const size_t value_sz = G_raster_size(map_type);
189
 
    unsigned int row;
190
 
    void *raster_row;
191
 
 
192
 
    raster_row = G_calloc(cols, value_sz);
193
 
 
194
 
    for (row = 0; row < rows; row++) {
195
 
        void *ptr;
196
 
        unsigned int col;
197
 
 
198
 
        if (G_get_raster_row(fd, raster_row, row, map_type) < 0)
199
 
            G_fatal_error(_("Reading row %d"), row);
200
 
 
201
 
        ptr = raster_row;
202
 
 
203
 
        for (col = 0; col < cols; col++) {
204
 
 
205
 
            if (G_is_null_value(ptr, map_type)) {
206
 
                ptr = G_incr_void_ptr(ptr, value_sz);
207
 
                continue;
208
 
            }
209
 
 
210
 
            if (nextp) {
211
 
                /* put the value into stats->XXXcell_array */
212
 
                memcpy(nextp, ptr, value_sz);
213
 
                nextp = G_incr_void_ptr(nextp, value_sz);
214
 
            }
215
 
 
216
 
            {
217
 
                double val = ((map_type == DCELL_TYPE) ? *((DCELL *) ptr)
218
 
                              : (map_type == FCELL_TYPE) ? *((FCELL *) ptr)
219
 
                              : *((CELL *) ptr));
220
 
 
221
 
                stats->sum += val;
222
 
                stats->sumsq += val * val;
223
 
                stats->sum_abs += fabs(val);
224
 
 
225
 
                if (first) {
226
 
                    stats->max = val;
227
 
                    stats->min = val;
228
 
                    first = FALSE;
229
 
                }
230
 
                else {
231
 
                    if (val > stats->max)
232
 
                        stats->max = val;
233
 
                    if (val < stats->min)
234
 
                        stats->min = val;
235
 
                }
236
 
            }
237
 
 
238
 
            ptr = G_incr_void_ptr(ptr, value_sz);
239
 
            stats->n++;
240
 
        }
241
 
        if (!(param.shell_style->answer))
242
 
            G_percent(row, rows, 2);
243
 
    }
244
 
}