← Back to branch summary

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

« back to all changes in this revision

Viewing changes to raster/wildfire/r.spreadpath/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/wildfire/r.spreadpath/main.c
1
 
 
2
 
/****************************************************************************
3
 
 *
4
 
 * MODULE:       r.spreadpath
5
 
 * AUTHOR(S):    Jianping Xu 1995: WIldfire SPread Simulation, WiSpS (original contributor)
6
 
 *               Markus Neteler <neteler itc.it>
7
 
 *               Roberto Flor <flor itc.it>, Brad Douglas <rez touchofmadness.com>,
8
 
 *               Glynn Clements <glynn gclements.plus.com>, Jachym Cepicky <jachym les-ejk.cz>
9
 
 * PURPOSE:      This is the main program for tracing out the shortest path(s)
10
 
 *               based on the raster map showing back path cells from which the   
11
 
 *               cumulative costs were determined.
12
 
 * COPYRIGHT:    (C) 2000-2006 by the GRASS Development Team
13
 
 *
14
 
 *               This program is free software under the GNU General Public
15
 
 *               License (>=v2). Read the file COPYING that comes with GRASS
16
 
 *               for details.
17
 
 *
18
 
 *****************************************************************************/
19
 
 
20
 
/**********************************************************************/
21
 
/*                                                                    */
22
 
/*      This is the main program for tracing out the shortest path(s) */
23
 
/*      based on the raster map showing back path cells from which the */
24
 
/*      cumulative costs were determined.                             */
25
 
/*                                                                    */
26
 
 
27
 
/**********************************************************************/
28
 
 
29
 
 
30
 
#include <stdio.h>
31
 
#include <stdlib.h>
32
 
#include <string.h>
33
 
#include <fcntl.h>
34
 
#include <grass/segment.h>
35
 
#include <grass/gis.h>
36
 
#define MAIN
37
 
#include "stash.h"
38
 
#include <grass/glocale.h>
39
 
#include "local_proto.h"
40
 
 
41
 
 
42
 
char *value;
43
 
int nrows, ncols;
44
 
SEGMENT in_row_seg, in_col_seg, out_seg;
45
 
 
46
 
 
47
 
int main(int argc, char **argv)
48
 
{
49
 
    int n, verbose = 1,
50
 
        backrow, backcol,
51
 
        col, row,
52
 
        len, flag,
53
 
        srows, scols,
54
 
        backrow_fd, backcol_fd, path_fd, in_row_fd, in_col_fd, out_fd;
55
 
    char *current_mapset,
56
 
        *search_mapset,
57
 
        *path_mapset,
58
 
        *backrow_mapset,
59
 
        *backcol_mapset, *in_row_file, *in_col_file, *out_file;
60
 
    CELL *cell;
61
 
    POINT *PRES_PT, *PRESENT_PT, *OLD_PT;
62
 
    struct Cell_head window;
63
 
    double east, north;
64
 
    struct Option *opt1, *opt2, *opt3, *opt4;
65
 
    struct Flag *flag1;
66
 
    struct GModule *module;
67
 
 
68
 
    G_gisinit(argv[0]);
69
 
 
70
 
    /* Set description */
71
 
    module = G_define_module();
72
 
    module->keywords = _("raster, fire");
73
 
    module->description =
74
 
        _("Recursively traces the least cost path backwards to "
75
 
          "cells from which the cumulative cost was determined.");
76
 
 
77
 
    opt1 = G_define_option();
78
 
    opt1->key = "x_input";
79
 
    opt1->type = TYPE_STRING;
80
 
    opt1->required = YES;
81
 
    opt1->gisprompt = "old,cell,raster";
82
 
    opt1->description =
83
 
        _("Name of raster map containing back-path easting information");
84
 
 
85
 
    opt2 = G_define_option();
86
 
    opt2->key = "y_input";
87
 
    opt2->type = TYPE_STRING;
88
 
    opt2->required = YES;
89
 
    opt2->gisprompt = "old,cell,raster";
90
 
    opt2->description =
91
 
        _("Name of raster map containing back-path northing information");
92
 
 
93
 
    opt3 = G_define_option();
94
 
    opt3->key = "coordinate";
95
 
    opt3->type = TYPE_STRING;
96
 
    opt3->multiple = YES;
97
 
    opt3->key_desc = "x,y";
98
 
    opt3->description =
99
 
        _("The map E and N grid coordinates of starting points");
100
 
 
101
 
    opt4 = G_define_option();
102
 
    opt4->key = "output";
103
 
    opt4->type = TYPE_STRING;
104
 
    opt4->required = YES;
105
 
    opt4->gisprompt = "new,cell,raster";
106
 
    opt4->description = _("Name of spread path raster map");
107
 
 
108
 
    flag1 = G_define_flag();
109
 
    flag1->key = 'v';
110
 
    flag1->description = _("Run verbosely");
111
 
 
112
 
    /*   Do command line parsing    */
113
 
    if (G_parser(argc, argv))
114
 
        exit(EXIT_FAILURE);
115
 
 
116
 
    current_mapset = G_mapset();
117
 
    in_row_file = G_tempfile();
118
 
    in_col_file = G_tempfile();
119
 
    out_file = G_tempfile();
120
 
 
121
 
    /*  Get database window parameters      */
122
 
    if (G_get_window(&window) < 0)
123
 
        G_fatal_error("can't read current window parameters");
124
 
 
125
 
    verbose = flag1->answer;
126
 
 
127
 
    /*  Check if backrow layer exists in data base  */
128
 
    search_mapset = "";
129
 
 
130
 
    strcpy(backrow_layer, opt2->answer);
131
 
    strcpy(backcol_layer, opt1->answer);
132
 
 
133
 
    backrow_mapset = G_find_cell(backrow_layer, search_mapset);
134
 
    backcol_mapset = G_find_cell(backcol_layer, search_mapset);
135
 
 
136
 
    if (backrow_mapset == NULL)
137
 
        G_fatal_error("%s - not found", backrow_layer);
138
 
 
139
 
    if (backcol_mapset == NULL)
140
 
        G_fatal_error("%s - not found", backcol_layer);
141
 
 
142
 
    search_mapset = "";
143
 
 
144
 
    strcpy(path_layer, opt4->answer);
145
 
 
146
 
    path_mapset = G_find_cell(path_layer, search_mapset);
147
 
 
148
 
    /*  find number of rows and cols in window    */
149
 
    nrows = G_window_rows();
150
 
    ncols = G_window_cols();
151
 
 
152
 
    cell = G_allocate_cell_buf();
153
 
 
154
 
    /*  Open back cell layers for reading  */
155
 
    backrow_fd = G_open_cell_old(backrow_layer, backrow_mapset);
156
 
    if (backrow_fd < 0)
157
 
        G_fatal_error("%s - can't open raster map", backrow_layer);
158
 
 
159
 
    backcol_fd = G_open_cell_old(backcol_layer, backcol_mapset);
160
 
    if (backcol_fd < 0)
161
 
        G_fatal_error("%s - can't open raster map", backcol_layer);
162
 
 
163
 
    /*   Parameters for map submatrices   */
164
 
    len = sizeof(CELL);
165
 
 
166
 
    srows = nrows / 4 + 1;
167
 
    scols = ncols / 4 + 1;
168
 
 
169
 
    if (verbose)
170
 
        G_message
171
 
            ("\nReading the input map -%s- and -%s- and creating some temporary files...",
172
 
             backrow_layer, backcol_layer);
173
 
 
174
 
    /* Create segmented files for back cell and output layers  */
175
 
    in_row_fd = creat(in_row_file, 0666);
176
 
    segment_format(in_row_fd, nrows, ncols, srows, scols, len);
177
 
    close(in_row_fd);
178
 
    in_col_fd = creat(in_col_file, 0666);
179
 
    segment_format(in_col_fd, nrows, ncols, srows, scols, len);
180
 
    close(in_col_fd);
181
 
 
182
 
    out_fd = creat(out_file, 0666);
183
 
    segment_format(out_fd, nrows, ncols, srows, scols, len);
184
 
    close(out_fd);
185
 
 
186
 
    /*   Open initialize and segment all files  */
187
 
    in_row_fd = open(in_row_file, 2);
188
 
    segment_init(&in_row_seg, in_row_fd, 4);
189
 
    in_col_fd = open(in_col_file, 2);
190
 
    segment_init(&in_col_seg, in_col_fd, 4);
191
 
 
192
 
    out_fd = open(out_file, 2);
193
 
    segment_init(&out_seg, out_fd, 4);
194
 
 
195
 
    /*   Write the back cell layers in the segmented files, and  
196
 
     *   Change UTM coordinates to ROWs and COLUMNs */
197
 
    for (row = 0; row < nrows; row++) {
198
 
        if (G_get_map_row(backrow_fd, cell, row) < 0)
199
 
            G_fatal_error("unable to get map row %d", row);
200
 
 
201
 
        for (col = 0; col < ncols; col++)
202
 
            if (cell[col] > 0)
203
 
                cell[col] =
204
 
                    (window.north - cell[col]) / window.ns_res /* - 0.5 */ ;
205
 
            else
206
 
                cell[col] = -1;
207
 
        segment_put_row(&in_row_seg, cell, row);
208
 
        if (G_get_map_row(backcol_fd, cell, row) < 0)
209
 
            G_fatal_error("unable to get map row %d", row);
210
 
 
211
 
        for (col = 0; col < ncols; col++)
212
 
            if (cell[col] > 0)
213
 
                cell[col] =
214
 
                    (cell[col] - window.west) / window.ew_res /* - 0.5 */ ;
215
 
        segment_put_row(&in_col_seg, cell, row);
216
 
    }
217
 
 
218
 
    /* Convert easting and northing from the command line to row and col */
219
 
    if (opt3->answer) {
220
 
        for (n = 0; opt3->answers[n] != NULL; n += 2) {
221
 
            G_scan_easting(opt3->answers[n], &east, G_projection());
222
 
            G_scan_northing(opt3->answers[n + 1], &north, G_projection());
223
 
            row = (window.north - north) / window.ns_res;
224
 
            col = (east - window.west) / window.ew_res;
225
 
            /* ignore pt outside window */
226
 
            if (east < window.west || east > window.east ||
227
 
                north < window.south || north > window.north) {
228
 
                G_warning("Ignoring point outside window: ");
229
 
                G_warning("   %.4f,%.4f", east, north);
230
 
                continue;
231
 
            }
232
 
 
233
 
            value = (char *)&backrow;
234
 
            segment_get(&in_row_seg, value, row, col);
235
 
            /* ignore pt in no-data area */
236
 
            if (backrow < 0) {
237
 
                G_warning("Ignoring point in NO-DATA area :");
238
 
                G_warning("   %.4f,%.4f", east, north);
239
 
                continue;
240
 
            }
241
 
            value = (char *)&backcol;
242
 
            segment_get(&in_col_seg, value, row, col);
243
 
 
244
 
            insert(&PRESENT_PT, row, col, backrow, backcol);
245
 
        }
246
 
    }
247
 
 
248
 
    /*  Set flag according to input */
249
 
    if (path_mapset != NULL) {
250
 
        if (head_start_pt == NULL)
251
 
            /*output layer exists and start pts are not given on cmd line */
252
 
            flag = 1;
253
 
 
254
 
        /* output layer exists and starting pts are given on cmd line */
255
 
        else
256
 
            flag = 2;
257
 
    }
258
 
    else
259
 
        flag = 3;               /* output layer does not previously exist */
260
 
 
261
 
 
262
 
    /*  Check if specified output layer name is legal   */
263
 
    if (flag == 3) {
264
 
        if (G_legal_filename(path_layer) < 0)
265
 
            G_fatal_error("%s - illegal name", path_layer);
266
 
    }
267
 
 
268
 
    /* If the output layer containing the starting positions */
269
 
    /* create a linked list of of them  */
270
 
    if (flag == 1) {
271
 
        path_fd = G_open_cell_old(path_layer, path_mapset);
272
 
        if (path_fd < 0)
273
 
            G_fatal_error("%s -can't open raster map", path_layer);
274
 
 
275
 
        /*  Search for the marked starting pts and make list    */
276
 
        for (row = 0; row < nrows; row++) {
277
 
            if (G_get_map_row(path_fd, cell, row) < 0)
278
 
                G_fatal_error("unable to get map row %d", row);
279
 
 
280
 
            for (col = 0; col < ncols; col++) {
281
 
                if (cell[col] > 0) {
282
 
                    value = (char *)&backrow;
283
 
                    segment_get(&in_row_seg, value, row, col);
284
 
                    /* ignore pt in no-data area */
285
 
                    if (backrow < 0) {
286
 
                        G_warning("Ignoring point in NO-DATA area:");
287
 
                        G_warning("   %.4f,%.4f\n",
288
 
                                  window.west + window.ew_res * (col + 0.5),
289
 
                                  window.north - window.ns_res * (row + 0.5));
290
 
                        continue;
291
 
                    }
292
 
                    value = (char *)&backcol;
293
 
                    segment_get(&in_col_seg, value, row, col);
294
 
                    insert(&PRESENT_PT, row, col, backrow, backcol);
295
 
                }
296
 
            }                   /* loop over cols */
297
 
        }                       /* loop over rows */
298
 
 
299
 
        G_close_cell(path_fd);
300
 
    }
301
 
 
302
 
    /* loop over the starting points to find the least cost paths */
303
 
    if (verbose)
304
 
        G_message("\nFinding the least cost paths ...");
305
 
 
306
 
    PRES_PT = head_start_pt;
307
 
    while (PRES_PT != NULL) {
308
 
        path_finder(PRES_PT->row, PRES_PT->col, PRES_PT->backrow,
309
 
                    PRES_PT->backcol);
310
 
 
311
 
        OLD_PT = PRES_PT;
312
 
        PRES_PT = NEXT_PT;
313
 
        G_free(OLD_PT);
314
 
    }
315
 
 
316
 
    /* Write pending updates by segment_put() to outputmap */
317
 
    segment_flush(&out_seg);
318
 
 
319
 
    if (verbose)
320
 
        G_message("\nWriting the output map  -%s-...", path_layer);
321
 
 
322
 
    path_fd = G_open_cell_new(path_layer);
323
 
    for (row = 0; row < nrows; row++) {
324
 
        segment_get_row(&out_seg, cell, row);
325
 
        if (G_put_raster_row(path_fd, cell, CELL_TYPE) < 0)
326
 
            G_fatal_error("unable to write map row %d", row);
327
 
    }
328
 
 
329
 
    if (verbose)
330
 
        G_message("finished.");
331
 
 
332
 
    segment_release(&in_row_seg);       /* release memory  */
333
 
    segment_release(&in_col_seg);
334
 
    segment_release(&out_seg);
335
 
 
336
 
    close(in_row_fd);           /* close all files */
337
 
    close(in_col_fd);
338
 
    close(out_fd);
339
 
 
340
 
    G_close_cell(path_fd);
341
 
    G_close_cell(backrow_fd);
342
 
    G_close_cell(backcol_fd);
343
 
 
344
 
    unlink(in_row_file);        /* remove submatrix files  */
345
 
    unlink(in_col_file);
346
 
    unlink(out_file);
347
 
 
348
 
    exit(EXIT_SUCCESS);
349
 
}