← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

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

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.volume/main.c
5
5
 * AUTHOR(S):    Dr. James Hinthorne, Central Washington University GIS Laboratory
6
6
 *               December 1988, (revised April 1989) (original contributor)
7
7
 *               Revised Jul 1995 to use new sites API (McCauley)
8
 
 *               GRASS 6 update: Hamish Bowman <hamish_nospam yahoo.com>
 
8
 *               GRASS 6 update: Hamish Bowman <hamish_b yahoo.com>
9
9
 *               Glynn Clements <glynn gclements.plus.com>, Soeren Gebbert <soeren.gebbert gmx.de>
 
10
 *               GRASS 7 update (to use Vlib): Martin Landa <landa.martin gmail.com>
10
11
 * PURPOSE:      
11
12
 *               r.volume is a program to compute the total, and average of cell values
12
13
 *               within regions of a map defined by clumps or patches on a second map
13
14
 *               (or MASK).  It also computes the "volume" by multiplying the total
14
 
 *               within a clump by the area of each cell.  It also outputs the
15
 
 *               "centroid" location of each clump.  Output is to standard out.
 
15
 *               within a clump by the area of each cell. It also outputs the
 
16
 *               "centroid" location of each clump. Output is to standard out.
16
17
 *
17
 
 * COPYRIGHT:    (C) 1999-2006 by the GRASS Development Team
 
18
 * COPYRIGHT:    (C) 1999-2006, 2013 by the GRASS Development Team
18
19
 *
19
20
 *               This program is free software under the GNU General Public
20
21
 *               License (>=v2). Read the file COPYING that comes with GRASS
25
26
#include <stdlib.h>
26
27
#include <string.h>
27
28
#include <grass/gis.h>
28
 
#include <grass/site.h>
 
29
#include <grass/raster.h>
 
30
#include <grass/vector.h>
 
31
#include <grass/dbmi.h>
29
32
#include <grass/glocale.h>
30
33
 
31
 
int centroids(int, int *, int *, int, int);
 
34
#include "local_proto.h"
32
35
 
33
36
int main(int argc, char *argv[])
34
37
{
35
 
    /* variables   */
36
 
    CELL *data_buf, *clump_buf;
 
38
    /* variables */
 
39
    DCELL *data_buf;
 
40
    CELL *clump_buf;
37
41
    CELL i, max;
 
42
 
38
43
    int row, col, rows, cols;
39
44
    int out_mode, use_MASK, *n, *e;
40
45
    long int *count;
41
46
    int fd_data, fd_clump;
42
 
    char buf[200], datamap[100], clumpmap[100], site_list[100];
43
 
    char *curr_mapset, *data_mapset, *clump_mapset;
 
47
 
 
48
    const char *datamap, *clumpmap, *centroidsmap;
 
49
    
44
50
    double avg, vol, total_vol, east, north, *sum;
 
51
 
45
52
    struct Cell_head window;
46
 
    struct Map_info *fd_sites = NULL;
47
 
    Site *mysite;
48
 
    Site_head site_info;
 
53
 
 
54
    struct Map_info *fd_centroids;
 
55
    struct line_pnts *Points;
 
56
    struct line_cats *Cats;
 
57
    struct field_info *Fi;
 
58
 
 
59
    char buf[DB_SQL_MAX];
 
60
    dbString sql;
 
61
    dbDriver *driver;
 
62
 
49
63
    struct GModule *module;
50
 
    struct Option *opt1, *opt2, *opt3;
51
 
    struct Flag *flag1;
 
64
    struct {
 
65
        struct Option *input, *clump, *centroids, *output;
 
66
    } opt;
 
67
    struct {
 
68
        struct Flag *report;
 
69
    } flag;
52
70
 
53
 
    /* Initialize GIS */
 
71
    /* define parameters and flags */
54
72
    G_gisinit(argv[0]);
55
73
 
56
74
    module = G_define_module();
57
 
    module->keywords = _("raster");
58
 
    module->description =
59
 
        _("Calculates the volume of data \"clumps\", "
60
 
          "and (optionally) produces a GRASS vector points map "
61
 
          "containing the calculated centroids of these clumps.");
62
 
 
63
 
    opt1 = G_define_option();
64
 
    opt1->key = "data";
65
 
    opt1->type = TYPE_STRING;
66
 
    opt1->required = YES;
67
 
    opt1->gisprompt = "old,cell,raster";
68
 
    opt1->description =
69
 
        _("Existing raster map representing data that will be summed within clumps");
70
 
 
71
 
    opt2 = G_define_option();
72
 
    opt2->key = "clump";
73
 
    opt2->type = TYPE_STRING;
74
 
    opt2->required = NO;
75
 
    opt2->gisprompt = "old,cell,raster";
76
 
    opt2->description =
77
 
        _("Existing raster map, preferably the output of r.clump");
78
 
 
79
 
    opt3 = G_define_option();
80
 
    opt3->key = "centroids";
81
 
    opt3->type = TYPE_STRING;
82
 
    opt3->required = NO;
83
 
    opt3->gisprompt = "new,vector,vector";
84
 
    opt3->description = _("Vector points map to contain clump centroids");
85
 
 
86
 
    flag1 = G_define_flag();
87
 
    flag1->key = 'f';
88
 
    flag1->description = _("Generate unformatted report");
 
75
    G_add_keyword(_("raster"));
 
76
    G_add_keyword(_("volume"));
 
77
    G_add_keyword(_("clumps"));
 
78
    module->label =
 
79
        _("Calculates the volume of data \"clumps\".");
 
80
    module->description = _("Optionally produces a GRASS vector points map "
 
81
                            "containing the calculated centroids of these clumps.");
 
82
 
 
83
    opt.input = G_define_standard_option(G_OPT_R_INPUT);
 
84
    opt.input->description =
 
85
        _("Name of input raster map representing data that will be summed within clumps");
 
86
 
 
87
    opt.clump = G_define_standard_option(G_OPT_R_INPUT);
 
88
    opt.clump->key = "clump";
 
89
    opt.clump->required = NO;
 
90
    opt.clump->label =
 
91
        _("Name of input clump raster map");
 
92
    opt.clump->description = _("Preferably the output of r.clump. "
 
93
                               "If no clump map is given than MASK is used.");
 
94
 
 
95
    opt.centroids = G_define_standard_option(G_OPT_V_OUTPUT);
 
96
    opt.centroids->key = "centroids";
 
97
    opt.centroids->required = NO;
 
98
    opt.centroids->description = _("Name for output vector points map to contain clump centroids");
 
99
 
 
100
    opt.output = G_define_standard_option(G_OPT_F_OUTPUT);
 
101
    opt.output->required = NO;
 
102
    opt.output->label =
 
103
        _("Name for output file to hold the report");
 
104
    opt.output->description =
 
105
        _("If no output file given report is printed to standard output");
 
106
 
 
107
    flag.report = G_define_flag();
 
108
    flag.report->key = 'f';
 
109
    flag.report->description = _("Generate unformatted report (items separated by colon)");
89
110
 
90
111
    if (G_parser(argc, argv))
91
112
        exit(EXIT_FAILURE);
92
113
 
93
 
    /* get current window */
94
 
    G_get_window(&window);
95
 
 
96
 
    /* initialize */
97
 
    *site_list = 0;
98
 
    out_mode = 1;               /* assume full output text */
99
 
    mysite = G_site_new_struct(CELL_TYPE, 2, 0, 4);
100
 
 
101
114
    /* get arguments */
102
 
    strcpy(datamap, opt1->answer);
103
 
 
104
 
    if (opt2->answer != NULL)
105
 
        strcpy(clumpmap, opt2->answer);
106
 
    else
107
 
        clumpmap[0] = '\0';
108
 
 
109
 
    if (opt3->answer != NULL)
110
 
        strcpy(site_list, opt3->answer);
111
 
    else
112
 
        site_list[0] = '\0';
113
 
 
114
 
    out_mode = (!flag1->answer);
115
 
 
116
 
    if (*datamap == 0)
117
 
        G_fatal_error(_("No data map specified"));
 
115
    datamap = opt.input->answer;
 
116
    
 
117
    clumpmap = NULL;
 
118
    if (opt.clump->answer)
 
119
        clumpmap = opt.clump->answer;
 
120
    
 
121
    centroidsmap = NULL;
 
122
    fd_centroids = NULL;
 
123
    Points = NULL;
 
124
    Cats = NULL;
 
125
    driver = NULL;
 
126
    if (opt.centroids->answer) {
 
127
        centroidsmap = opt.centroids->answer;
 
128
        fd_centroids = G_malloc(sizeof(struct Map_info));
 
129
    }
 
130
    
 
131
    out_mode = (!flag.report->answer);
118
132
 
119
133
    /*
120
 
     * See if MASK or a separate "clumpmap" layer is to be used-- it must(!)
121
 
     * be one of those two choices.
 
134
     * see if MASK or a separate "clumpmap" raster map is to be used
 
135
     * -- it must(!) be one of those two choices.
122
136
     */
123
137
    use_MASK = 0;
124
 
    if (*clumpmap == '\0') {
125
 
        strcpy(clumpmap, "MASK");
 
138
    if (!clumpmap) {
 
139
        clumpmap = "MASK";
126
140
        use_MASK = 1;
127
 
        if (G_find_cell(clumpmap, G_mapset()) == NULL)
128
 
            G_fatal_error(_("No clump map specified and MASK not set."));
129
 
    }
130
 
    curr_mapset = G_mapset();
131
 
    data_mapset = G_find_cell2(datamap, "");
132
 
    if (!data_mapset)
133
 
        G_fatal_error(_("Unable to find data map"));
134
 
 
135
 
    fd_data = G_open_cell_old(datamap, data_mapset);
136
 
    if (use_MASK)
137
 
        clump_mapset = curr_mapset;
138
 
    else
139
 
        clump_mapset = G_find_cell2(clumpmap, "");
140
 
    if (!clump_mapset)
141
 
        G_fatal_error(_("Unable to find clump map"));
142
 
 
143
 
    fd_clump = G_open_cell_old(clumpmap, clump_mapset);
144
 
 
145
 
    /* initialize sites file (for centroids) if needed */
146
 
    if (*site_list) {
147
 
        fd_sites = G_fopen_sites_new(site_list);
148
 
        if (fd_sites == NULL)
149
 
            G_fatal_error(_("Unable to open centroids vector points map"));
 
141
        if (!G_find_raster2(clumpmap, G_mapset()))
 
142
            G_fatal_error(_("No MASK found. If no clump map is given than the MASK is required. "
 
143
                            "You need to define a clump raster map or create a MASK by r.mask command."));
 
144
        G_important_message(_("No clump map given, using MASK"));
 
145
    }
 
146
    
 
147
    /* open input and clump raster maps */
 
148
    fd_data = Rast_open_old(datamap, "");
 
149
    fd_clump = Rast_open_old(clumpmap, use_MASK ? G_mapset() : "");
 
150
    
 
151
    /* initialize vector map (for centroids) if needed */
 
152
    if (centroidsmap) {
 
153
        if (Vect_open_new(fd_centroids, centroidsmap, WITHOUT_Z) < 0)
 
154
            G_fatal_error(_("Unable to create vector map <%s>"), centroidsmap);
 
155
        
 
156
        Points = Vect_new_line_struct();
 
157
        Cats = Vect_new_cats_struct();
 
158
        
 
159
        /* initialize data structures */
 
160
        Vect_append_point(Points, 0., 0., 0.);
 
161
        Vect_cat_set(Cats, 1, 1);
 
162
    }
 
163
    
 
164
    /* initialize output file */
 
165
    if (opt.output->answer && strcmp(opt.output->answer, "-") != 0) {
 
166
        if (freopen(opt.output->answer, "w", stdout) == NULL) {
 
167
            perror(opt.output->answer);
 
168
            exit(EXIT_FAILURE);
 
169
        }
150
170
    }
151
171
 
152
172
    /* initialize data accumulation arrays */
153
 
    max = G_number_of_cats(clumpmap, clump_mapset);
 
173
    max = Rast_get_max_c_cat(clumpmap, use_MASK ? G_mapset() : "");
154
174
 
155
175
    sum = (double *)G_malloc((max + 1) * sizeof(double));
156
176
    count = (long int *)G_malloc((max + 1) * sizeof(long int));
157
177
 
158
 
    for (i = 0; i <= max; i++) {
159
 
        sum[i] = 0;
160
 
        count[i] = 0;
161
 
    }
162
 
 
163
 
    data_buf = G_allocate_cell_buf();
164
 
    clump_buf = G_allocate_cell_buf();
165
 
 
 
178
    G_zero(sum, (max + 1) * sizeof(double));
 
179
    G_zero(count, (max + 1) * sizeof(long int));
 
180
    
 
181
    data_buf = Rast_allocate_d_buf();
 
182
    clump_buf = Rast_allocate_c_buf();
 
183
    
166
184
    /* get window size */
 
185
    G_get_window(&window);
167
186
    rows = window.rows;
168
187
    cols = window.cols;
169
188
 
170
 
    if (fd_data < 0 || fd_clump < 0)
171
 
        G_fatal_error(_("Data or Clump file not open"));
172
 
 
173
189
    /* now get the data -- first pass */
174
 
    G_message("Complete ...");
175
190
    for (row = 0; row < rows; row++) {
176
191
        G_percent(row, rows, 2);
177
 
        G_get_map_row(fd_data, data_buf, row);
178
 
        G_get_map_row(fd_clump, clump_buf, row);
 
192
        Rast_get_d_row(fd_data, data_buf, row);
 
193
        Rast_get_c_row(fd_clump, clump_buf, row);
179
194
        for (col = 0; col < cols; col++) {
180
195
            i = clump_buf[col];
181
 
            if (i > max) {
182
 
                sprintf(buf,
183
 
                        "Row=%d Col=%d Cat=%d in clump map [%s]; max=%d.\n",
184
 
                        row, col, i, clumpmap, max);
185
 
                strcat(buf, "Cat value > max returned by G_number_of_cats.");
186
 
                G_fatal_error(buf);
187
 
            }
188
 
            if (i < 1)
 
196
            if (i > max)
 
197
                G_fatal_error(_("Invalid category value %d (max=%d): row=%d col=%d"),
 
198
                              i, max, row, col);
 
199
            if (i < 1) {
 
200
                G_debug(3, "row=%d col=%d: zero or negs ignored", row, col);
189
201
                continue;       /* ignore zeros and negs */
 
202
            }
 
203
            if (Rast_is_d_null_value(&data_buf[col])) {
 
204
                G_debug(3, "row=%d col=%d: NULL ignored", row, col);
 
205
                continue;
 
206
            }
 
207
            
 
208
            sum[i] += data_buf[col];
190
209
            count[i]++;
191
 
            sum[i] += data_buf[col];
192
210
        }
193
211
    }
194
 
    G_percent(row, rows, 2);
 
212
    G_percent(1, 1, 1);
 
213
    
195
214
    /* free some buffer space */
196
215
    G_free(data_buf);
197
216
    G_free(clump_buf);
202
221
 
203
222
    i = centroids(fd_clump, e, n, 1, max);
204
223
 
 
224
    /* close raster maps */
 
225
    Rast_close(fd_data);
 
226
    Rast_close(fd_clump);
 
227
    
205
228
    /* got everything, now do output */
206
 
    if (*site_list) {
207
 
        char desc[GNAME_MAX * 2 + 40];
208
 
 
209
 
        site_info.form = NULL;
210
 
        site_info.time = NULL;
211
 
        site_info.stime = NULL;
212
 
        sprintf(desc, "from %s on map %s using clumps from %s",
213
 
                argv[0], datamap, clumpmap);
214
 
        site_info.desc = G_store(desc);
215
 
        site_info.name = G_store(site_list);
216
 
        site_info.labels =
217
 
            G_store("centroid east|centroid north|#cat vol avg t n");
218
 
        G_site_put_head(fd_sites, &site_info);
 
229
    if (centroidsmap) {
 
230
        G_message(_("Creating vector point map <%s>..."), centroidsmap);
 
231
        /* set comment */
 
232
        sprintf(buf, _("From '%s' on raster map <%s> using clumps from <%s>"),
 
233
                argv[0], datamap, clumpmap);
 
234
        Vect_set_comment(fd_centroids, buf);
 
235
 
 
236
        /* create attribute table */        
 
237
        Fi = Vect_default_field_info(fd_centroids, 1, NULL, GV_1TABLE);
 
238
        
 
239
        driver = db_start_driver_open_database(Fi->driver,
 
240
                                               Vect_subst_var(Fi->database, fd_centroids));
 
241
        if (driver == NULL) {
 
242
            G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
 
243
                          Vect_subst_var(Fi->database, fd_centroids), Fi->driver);
 
244
        }
 
245
        db_set_error_handler_driver(driver);
 
246
        
 
247
        db_begin_transaction(driver);
 
248
        
 
249
        db_init_string(&sql);
 
250
        sprintf(buf, "create table %s (cat integer, volume double precision, "
 
251
                "average double precision, sum double precision, count integer)",
 
252
                Fi->table);
 
253
        db_set_string(&sql, buf);
 
254
        Vect_map_add_dblink(fd_centroids, 1, NULL, Fi->table, GV_KEY_COLUMN, Fi->database,
 
255
                            Fi->driver);
 
256
 
 
257
        G_debug(3, "%s", db_get_string(&sql));
 
258
        if (db_execute_immediate(driver, &sql) != DB_OK) {
 
259
            G_fatal_error(_("Unable to create table: %s"), db_get_string(&sql));
 
260
        }
219
261
    }
 
262
 
 
263
    /* print header */
220
264
    if (out_mode) {
221
 
        fprintf(stdout, "Volume report on data from %s", datamap);
222
 
        fprintf(stdout, " using clumps on %s map\n\n", clumpmap);
223
 
        fprintf(stdout,
224
 
                " Cat    Average   Data   # Cells        Centroid             Total\n");
225
 
        fprintf(stdout,
226
 
                "Number  in clump  Total  in clump   Easting   Northing       Volume\n\n");
 
265
        fprintf(stdout, _("\nVolume report on data from <%s> using clumps on <%s> raster map"),
 
266
                datamap, clumpmap);
 
267
        fprintf(stdout, "\n\n");
 
268
        fprintf(stdout,
 
269
                _("Category   Average   Data   # Cells        Centroid             Total\n"));
 
270
        fprintf(stdout,
 
271
                _("Number     in clump  Total  in clump   Easting     Northing     Volume"));
 
272
        fprintf(stdout, "\n%s\n", SEP);
227
273
    }
228
274
    total_vol = 0.0;
229
275
 
 
276
    /* print output, write centroids */
230
277
    for (i = 1; i <= max; i++) {
231
278
        if (count[i]) {
232
279
            avg = sum[i] / (double)count[i];
234
281
            total_vol += vol;
235
282
            east = window.west + (e[i] + 0.5) * window.ew_res;
236
283
            north = window.north - (n[i] + 0.5) * window.ns_res;
237
 
            if (*site_list) {
238
 
                mysite->east = east;
239
 
                mysite->north = north;
240
 
                mysite->ccat = i;
241
 
                mysite->dbl_att[0] = vol;
242
 
                mysite->dbl_att[1] = avg;
243
 
                mysite->dbl_att[2] = sum[i];
244
 
                mysite->dbl_att[3] = (double)count[i];
245
 
                /*       "%-1.2f|%-1.2f|#%5d v=%-1.2f a=%-1.2f t=%-1.0f n=%ld\n", */
246
 
                /* east, north, i, vol, avg, sum[i], count[i]); */
247
 
                G_site_put(fd_sites, mysite);
 
284
            if (fd_centroids) { /* write centroids if requested */
 
285
                Points->x[0] = east;
 
286
                Points->y[0] = north;
 
287
                Cats->cat[0] = i;
 
288
                Vect_write_line(fd_centroids, GV_POINT, Points, Cats);
 
289
        
 
290
                sprintf(buf, "insert into %s values (%d, %f, %f, %f, %ld)",
 
291
                        Fi->table, i, vol, avg, sum[i], count[i]);
 
292
                db_set_string(&sql, buf);
 
293
 
 
294
                if (db_execute_immediate(driver, &sql) != DB_OK)
 
295
                    G_fatal_error(_("Cannot insert new row: %s"),
 
296
                                  db_get_string(&sql));
248
297
            }
249
298
            if (out_mode)
250
299
                fprintf(stdout,
251
 
                        "%5d%10.2f%10.0f %7ld  %10.2f  %10.2f %16.2f\n", i,
 
300
                        "%8d%10.2f%10.0f %7ld  %10.2f  %10.2f %16.2f\n", i,
252
301
                        avg, sum[i], count[i], east, north, vol);
253
302
            else
254
303
                fprintf(stdout, "%d:%.2f:%.0f:%ld:%.2f:%.2f:%.2f\n",
255
304
                        i, avg, sum[i], count[i], east, north, vol);
256
305
        }
257
306
    }
258
 
    if (total_vol > 0.0 && out_mode)
259
 
        fprintf(stdout, "%58s %14.2f", "Total Volume =", total_vol);
260
 
    fprintf(stdout, "\n");
 
307
 
 
308
    /* write centroid attributes and close the map*/
 
309
    if (fd_centroids) {
 
310
        db_commit_transaction(driver);
 
311
        Vect_close(fd_centroids);
 
312
    }
 
313
    
 
314
    /* print total value */
 
315
    if (total_vol > 0.0 && out_mode) {
 
316
        fprintf(stdout, "%s\n", SEP);
 
317
        fprintf(stdout, "%60s = %14.2f", _("Total Volume"), total_vol);
 
318
        fprintf(stdout, "\n");
 
319
    }
 
320
 
261
321
    exit(EXIT_SUCCESS);
262
 
}                               /* end of main() */
 
322
}