1
2 /****************************************************************************
3 *
4 * MODULE: r.spread
5 * AUTHOR(S): Jianping Xu, 1995 (original contributor)
6 * Center for Remote Sensing and Spatial Analysis (CRSSA)
7 * Rutgers University.
8 * Andreas.Lange <andreas.lange rhein-main.de>, Eric G. Miller <egm2 jps.net>,
9 * Markus Neteler <neteler itc.it>, Roberto Flor <flor itc.it>,
10 * Brad Douglas <rez touchofmadness.com>,
11 * Glynn Clements <glynn gclements.plus.com>, Jachym Cepicky <jachym les-ejk.cz>
12 * PURPOSE:
13 * This is the main program for simulating elliptical spread.
14 *
15 * It
16 * 1) determines the earliest time a phenomenon REACHES to a
17 * map cell, NOT the time that cell is EXHAUSTED.
18 * 3) If a cell is spread barrier, a no-data value is assigned
19 * to it.
20 * COPYRIGHT: (C) 2000-2006 by the GRASS Development Team
21 *
22 * This program is free software under the GNU General Public
23 * License (>=v2). Read the file COPYING that comes with GRASS
24 * for details.
25 *
26 *****************************************************************************/
27
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <math.h>
31 #include <sys/types.h>
32 #include <unistd.h>
33 #include <grass/gis.h>
34 #include <grass/raster.h>
35 #include <grass/glocale.h>
36 #include "cmd_line.h"
37 #include "costHa.h"
38 #include "cell_ptrHa.h"
39 #include "local_proto.h"
40
41 #define DATA(map, r, c) (map)[(r) * ncols + (c)]
42
43 CELL *cell;
44 CELL *x_cell;
45 CELL *y_cell;
46
47 CELL *map_max;
48 CELL *map_dir;
49 CELL *map_base;
50 CELL *map_spotdist;
51 CELL *map_velocity;
52 CELL *map_mois;
53 float *map_out;
54 CELL *map_x_out;
55 CELL *map_y_out;
56 CELL *map_visit;
57
58 char buf[400];
59
60 float zero = 0.0;
61 float neg = -2.0;
62
63 int BARRIER = 0;
64 int max_fd, dir_fd, base_fd, start_fd;
65 int spotdist_fd, velocity_fd, mois_fd;
66 int cum_fd, x_fd, y_fd;
67 int nrows, ncols;
68
69 long heap_len;
70
71 struct Cell_head window;
72
73 struct costHa *heap;
74
75
main(int argc,char * argv[])76 int main(int argc, char *argv[])
77 {
78 int col, row;
79
80 /* to menage start (source) raster map */
81 struct Range start_range;
82 CELL start_range_min, start_range_max;
83 int start_is_time; /* 0 or 1 */
84
85 struct
86 {
87 struct Option *max, *dir, *base, *start,
88 *spotdist, *velocity, *mois,
89 *least, *comp_dens, *init_time,
90 *time_lag, *backdrop, *out, *x_out, *y_out;
91 } parm;
92 struct
93 {
94 /* please, remove display before GRASS 7 released */
95 struct Flag *display, *spotting, *start_is_time;
96 } flag;
97 struct GModule *module;
98
99 /* initialize access to database and create temporary files */
100
101 G_gisinit(argv[0]);
102
103 /* Set description */
104 module = G_define_module();
105 G_add_keyword(_("raster"));
106 G_add_keyword(_("fire"));
107 G_add_keyword(_("spread"));
108 G_add_keyword(_("hazard"));
109 G_add_keyword(_("model"));
110 module->label =
111 _("Simulates elliptically anisotropic spread.");
112 module->description =
113 _("Generates a raster map of the cumulative time of spread, "
114 "given raster maps containing the rates of spread (ROS), "
115 "the ROS directions and the spread origins. "
116 "It optionally produces raster maps to contain backlink UTM "
117 "coordinates for tracing spread paths. "
118 "Usable for fire spread simulations.");
119
120 parm.base = G_define_option();
121 parm.base->key = "base_ros";
122 parm.base->type = TYPE_STRING;
123 parm.base->required = YES;
124 parm.base->gisprompt = "old,cell,raster";
125 parm.base->guisection = _("Input");
126 parm.base->label =
127 _("Raster map containing base ROS (cm/min)");
128 parm.base->description =
129 _("Name of an existing raster map layer in the user's "
130 "current mapset search path containing the ROS values in the directions "
131 "perpendicular to maximum ROSes' (cm/minute). These ROSes are also the ones "
132 "without the effect of directional factors.");
133
134 parm.max = G_define_option();
135 parm.max->key = "max_ros";
136 parm.max->type = TYPE_STRING;
137 parm.max->required = YES;
138 parm.max->gisprompt = "old,cell,raster";
139 parm.max->guisection = _("Input");
140 parm.max->label =
141 _("Raster map containing maximal ROS (cm/min)");
142 parm.max->description =
143 _("Name of an existing raster map layer in the user's current "
144 "mapset search path containing the maximum ROS values (cm/minute).");
145
146 parm.dir = G_define_option();
147 parm.dir->key = "direction_ros";
148 parm.dir->type = TYPE_STRING;
149 parm.dir->required = YES;
150 parm.dir->gisprompt = "old,cell,raster";
151 parm.dir->guisection = _("Input");
152 parm.dir->label =
153 _("Raster map containing directions of maximal ROS (degree)");
154 parm.dir->description =
155 _("Name of an existing raster map layer in the user's "
156 "current mapset search path containing directions of the maximum ROSes, "
157 "clockwise from north (degree)."); /* TODO: clockwise from north? see r.ros */
158
159 parm.start = G_define_option();
160 parm.start->key = "start";
161 parm.start->type = TYPE_STRING;
162 parm.start->required = YES;
163 parm.start->gisprompt = "old,cell,raster";
164 parm.start->guisection = _("Input");
165 parm.start->label =
166 _("Raster map containing starting sources");
167 parm.start->description =
168 _("Name of an existing raster map layer in the "
169 "user's current mapset search path containing starting locations of the "
170 "spread phenomenon. Any positive integers in this map are recognized as "
171 "starting sources (seeds).");
172
173 parm.spotdist = G_define_option();
174 parm.spotdist->key = "spotting_distance";
175 parm.spotdist->type = TYPE_STRING;
176 parm.spotdist->gisprompt = "old,cell,raster";
177 parm.spotdist->guisection = _("Input");
178 parm.spotdist->label =
179 _("Raster map containing maximal spotting distance (m, required with -s)");
180 parm.spotdist->description =
181 _("Name of an existing raster map layer in "
182 "the user's current mapset search path containing the maximum potential "
183 "spotting distances (meters).");
184
185 parm.velocity = G_define_option();
186 parm.velocity->key = "wind_speed";
187 parm.velocity->type = TYPE_STRING;
188 parm.velocity->gisprompt = "old,cell,raster";
189 parm.velocity->guisection = _("Input");
190 parm.velocity->label =
191 _("Raster map containing midflame wind speed (ft/min, required with -s)");
192 parm.velocity->description =
193 _("Name of an existing raster map layer in the "
194 "user's current mapset search path containing wind velocities at half of "
195 "the average flame height (feet/minute).");
196
197 parm.mois = G_define_option();
198 parm.mois->key = "fuel_moisture";
199 parm.mois->type = TYPE_STRING;
200 parm.mois->gisprompt = "old,cell,raster";
201 parm.mois->guisection = _("Input");
202 parm.mois->label =
203 _("Raster map containing fine fuel moisture of the cell receiving a spotting firebrand (%, required with -s)");
204 parm.mois->description =
205 _("Name of an existing raster map layer in the "
206 "user's current mapset search path containing the 1-hour (<.25\") fuel "
207 "moisture (percentage content multiplied by 100).");
208
209 parm.least = G_define_option();
210 parm.least->key = "least_size";
211 parm.least->type = TYPE_STRING;
212 parm.least->key_desc = "odd int";
213 parm.least->options = "3,5,7,9,11,13,15";
214 parm.least->label =
215 _("Basic sampling window size needed to meet certain accuracy (3)"); /* TODO: what is 3 here? default? */
216 parm.least->description =
217 _("An odd integer ranging 3 - 15 indicating "
218 "the basic sampling window size within which all cells will be considered "
219 "to see whether they will be reached by the current spread cell. The default "
220 "number is 3 which means a 3x3 window.");
221
222 parm.comp_dens = G_define_option();
223 parm.comp_dens->key = "comp_dens";
224 parm.comp_dens->type = TYPE_STRING;
225 parm.comp_dens->key_desc = "decimal";
226 parm.comp_dens->label =
227 _("Sampling density for additional computing (range: 0.0 - 1.0 (0.5))"); /* TODO: again, what is 0.5?, TODO: range not set */
228 parm.comp_dens->description =
229 _("A decimal number ranging 0.0 - 1.0 indicating "
230 "additional sampling cells will be considered to see whether they will be "
231 "reached by the current spread cell. The closer to 1.0 the decimal number "
232 "is, the longer the program will run and the higher the simulation accuracy "
233 "will be. The default number is 0.5.");
234
235 parm.init_time = G_define_option();
236 parm.init_time->key = "init_time";
237 parm.init_time->type = TYPE_STRING;
238 parm.init_time->key_desc = "int (>= 0)"; /* TODO: move to ->options */
239 parm.init_time->answer = "0";
240 parm.init_time->label =
241 _("Initial time for current simulation (0) (min)");
242 parm.init_time->description =
243 _("A non-negative number specifying the initial "
244 "time for the current spread simulation (minutes). This is useful when multiple "
245 "phase simulation is conducted. The default time is 0.");
246
247 parm.time_lag = G_define_option();
248 parm.time_lag->key = "lag";
249 parm.time_lag->type = TYPE_STRING;
250 parm.time_lag->key_desc = "int (>= 0)"; /* TODO: move to ->options */
251 parm.time_lag->label =
252 _("Simulating time duration LAG (fill the region) (min)"); /* TODO: what does this mean? */
253 parm.time_lag->description =
254 _("A non-negative integer specifying the simulating "
255 "duration time lag (minutes). The default is infinite, but the program will "
256 "terminate when the current geographic region/mask has been filled. It also "
257 "controls the computational time, the shorter the time lag, the faster the "
258 "program will run.");
259
260 /* TODO: what's this? probably display, so remove */
261 parm.backdrop = G_define_option();
262 parm.backdrop->key = "backdrop";
263 parm.backdrop->type = TYPE_STRING;
264 parm.backdrop->gisprompt = "old,cell,raster";
265 parm.backdrop->label =
266 _("Name of raster map as a display backdrop");
267 parm.backdrop->description =
268 _("Name of an existing raster map layer in the "
269 "user's current mapset search path to be used as the background on which "
270 "the \"live\" movement will be shown.");
271
272 parm.out = G_define_option();
273 parm.out->key = "output";
274 parm.out->type = TYPE_STRING;
275 parm.out->required = YES;
276 parm.out->gisprompt = "new,cell,raster";
277 parm.out->guisection = _("Output");
278 parm.out->label =
279 _("Raster map to contain output spread time (min)");
280 parm.out->description =
281 _("Name of the new raster map layer to contain "
282 "the results of the cumulative spread time needed for a phenomenon to reach "
283 "each cell from the starting sources (minutes).");
284
285 parm.x_out = G_define_option();
286 parm.x_out->key = "x_output";
287 parm.x_out->type = TYPE_STRING;
288 parm.x_out->gisprompt = "new,cell,raster";
289 parm.x_out->guisection = _("Output");
290 parm.x_out->label =
291 _("Name of raster map to contain X back coordinates");
292 parm.x_out->description =
293 _("Name of the new raster map layer to contain "
294 "the results of backlink information in UTM easting coordinates for each "
295 "cell.");
296
297 parm.y_out = G_define_option();
298 parm.y_out->key = "y_output";
299 parm.y_out->type = TYPE_STRING;
300 parm.y_out->gisprompt = "new,cell,raster";
301 parm.y_out->guisection = _("Output");
302 parm.y_out->label =
303 _("Name of raster map to contain Y back coordinates");
304 parm.y_out->description =
305 _("Name of the new raster map layer to contain "
306 "the results of backlink information in UTM northing coordinates for each "
307 "cell.");
308
309 #if 0
310 flag.display = G_define_flag();
311 flag.display->key = 'd';
312 flag.display->label = _("DISPLAY 'live' spread process on screen");
313 flag.display->description =
314 _("Display the 'live' simulation on screen. A graphics window "
315 "must be opened and selected before using this option.");
316 #endif
317
318 flag.spotting = G_define_flag();
319 flag.spotting->key = 's';
320 flag.spotting->description = _("Consider spotting effect (for wildfires)");
321
322 flag.start_is_time = G_define_flag();
323 flag.start_is_time->key = 'i';
324 flag.start_is_time->label = _("Use start raster map values in"
325 " output spread time raster map");
326 flag.start_is_time->description = _("Designed to be used with output"
327 " of previous run of r.spread when computing spread iteratively."
328 " The values in start raster map are considered as time."
329 " Allowed values in raster map are from zero"
330 " to the value of init_time option."
331 " If not enabled, init_time is used in the area of start raster map");
332
333 /* Parse command line */
334 if (G_parser(argc, argv))
335 exit(EXIT_FAILURE);
336
337
338 /* FIXME - allow seed to be specified for repeatability */
339 G_srand48_auto();
340
341 display = 0;
342
343 spotting = flag.spotting->answer;
344
345 max_layer = parm.max->answer;
346 dir_layer = parm.dir->answer;
347 base_layer = parm.base->answer;
348 start_layer = parm.start->answer;
349 backdrop_layer = parm.backdrop->answer;
350 out_layer = parm.out->answer;
351 if (parm.x_out->answer) {
352 x_out = 1;
353 x_out_layer = parm.x_out->answer;
354 }
355 if (parm.y_out->answer) {
356 y_out = 1;
357 y_out_layer = parm.y_out->answer;
358 }
359 if (spotting) {
360 if (!
361 (parm.spotdist->answer && parm.velocity->answer &&
362 parm.mois->answer)) {
363 G_warning
364 ("SPOTTING DISTANCE, fuel MOISTURE, or wind VELOCITY map not given w/ -s");
365 G_usage();
366 exit(EXIT_FAILURE);
367 }
368 else {
369 spotdist_layer = parm.spotdist->answer;
370 velocity_layer = parm.velocity->answer;
371 mois_layer = parm.mois->answer;
372 }
373 }
374 /*Check the given the least sampling size, assign the default if needed */
375 if (parm.least->answer)
376 least = atoi(parm.least->answer);
377 else
378 least = 3;
379 /*Check the given computing density, assign the default if needed */
380 if (parm.comp_dens->answer) {
381 comp_dens = atof(parm.comp_dens->answer);
382 if (comp_dens < 0.0 || comp_dens > 1.0) {
383 G_warning("Illegal computing density <%s>",
384 parm.comp_dens->answer);
385 G_usage();
386 exit(EXIT_FAILURE);
387 }
388 }
389 else {
390 comp_dens = 0.5;
391 }
392 /*Check the given initial time and simulation time lag, assign the default if needed */
393 init_time = atoi(parm.init_time->answer);
394 if (init_time < 0) {
395 G_warning("Illegal initial time <%s>", parm.init_time->answer);
396 G_usage();
397 exit(EXIT_FAILURE);
398 }
399
400 if (parm.time_lag->answer) {
401 time_lag = atoi(parm.time_lag->answer);
402 if (time_lag < 0) {
403 G_warning("Illegal simulating time lag <%s>",
404 parm.time_lag->answer);
405 G_usage();
406 exit(EXIT_FAILURE);
407 }
408 }
409 else {
410 time_lag = 99999;
411 }
412
413 /* Get database window parameters */
414
415 G_get_window(&window);
416
417 /* find number of rows and columns in window */
418
419 nrows = Rast_window_rows();
420 ncols = Rast_window_cols();
421
422 /*transfor measurement unit from meters to centimeters due to ROS unit
423 *if the input ROSs are in m/min units, cancell the following*/
424 window.ns_res = 100 * window.ns_res;
425 window.ew_res = 100 * window.ew_res;
426
427 /* Initialize display screens */
428 #if 0
429 if (display)
430 display_init();
431 #endif
432
433 /* Check if input layers exists in data base */
434
435 if (G_find_raster2(max_layer, "") == NULL)
436 G_fatal_error("Raster map <%s> not found", max_layer);
437
438 if (G_find_raster2(dir_layer, "") == NULL)
439 G_fatal_error(_("Raster map <%s> not found"), dir_layer);
440
441 if (G_find_raster2(base_layer, "") == NULL)
442 G_fatal_error(_("Raster map <%s> not found"), base_layer);
443
444 if (G_find_raster2(start_layer, "") == NULL)
445 G_fatal_error(_("Raster map <%s> not found"), start_layer);
446
447 if (spotting) {
448 if (G_find_raster2(spotdist_layer, "") == NULL)
449 G_fatal_error(_("Raster map <%s> not found"), spotdist_layer);
450
451 if (G_find_raster2(velocity_layer, "") == NULL)
452 G_fatal_error(_("Raster map <%s> not found"), velocity_layer);
453
454 if (G_find_raster2(mois_layer, "") == NULL)
455 G_fatal_error(_("Raster map <%s> not found"), mois_layer);
456 }
457
458 /* Open input cell layers for reading */
459
460 max_fd = Rast_open_old(max_layer, G_find_raster2(max_layer, ""));
461
462 dir_fd = Rast_open_old(dir_layer, G_find_raster2(dir_layer, ""));
463
464 base_fd = Rast_open_old(base_layer, G_find_raster2(base_layer, ""));
465
466 if (spotting) {
467 spotdist_fd =
468 Rast_open_old(spotdist_layer, G_find_raster2(spotdist_layer, ""));
469
470 velocity_fd =
471 Rast_open_old(velocity_layer, G_find_raster2(velocity_layer, ""));
472
473 mois_fd = Rast_open_old(mois_layer, G_find_raster2(mois_layer, ""));
474 }
475
476 /* Allocate memories for a row */
477 cell = Rast_allocate_c_buf();
478 if (x_out)
479 x_cell = Rast_allocate_c_buf();
480 if (y_out)
481 y_cell = Rast_allocate_c_buf();
482
483 /* Allocate memories for a map */
484 map_max = (CELL *) G_calloc(nrows * ncols + 1, sizeof(CELL));
485 map_dir = (CELL *) G_calloc(nrows * ncols + 1, sizeof(CELL));
486 map_base = (CELL *) G_calloc(nrows * ncols + 1, sizeof(CELL));
487 map_visit = (CELL *) G_calloc(nrows * ncols + 1, sizeof(CELL));
488 map_out = (float *)G_calloc(nrows * ncols + 1, sizeof(float));
489 if (spotting) {
490 map_spotdist = (CELL *) G_calloc(nrows * ncols + 1, sizeof(CELL));
491 map_velocity = (CELL *) G_calloc(nrows * ncols + 1, sizeof(CELL));
492 map_mois = (CELL *) G_calloc(nrows * ncols + 1, sizeof(CELL));
493 }
494 if (x_out)
495 map_x_out = (CELL *) G_calloc(nrows * ncols + 1, sizeof(CELL));
496 if (y_out)
497 map_y_out = (CELL *) G_calloc(nrows * ncols + 1, sizeof(CELL));
498
499
500 /* Write the input layers in the map "arrays" */
501
502 G_message(_("Reading inputs..."));
503
504 for (row = 0; row < nrows; row++) {
505 G_percent(row, nrows, 2);
506 Rast_get_c_row(max_fd, cell, row);
507 for (col = 0; col < ncols; col++)
508 DATA(map_max, row, col) = cell[col];
509 Rast_get_c_row(dir_fd, cell, row);
510 for (col = 0; col < ncols; col++)
511 DATA(map_dir, row, col) = cell[col];
512 Rast_get_c_row(base_fd, cell, row);
513 for (col = 0; col < ncols; col++)
514 DATA(map_base, row, col) = cell[col];
515 if (spotting) {
516 Rast_get_c_row(spotdist_fd, cell, row);
517 for (col = 0; col < ncols; col++)
518 DATA(map_spotdist, row, col) = cell[col];
519 Rast_get_c_row(velocity_fd, cell, row);
520 for (col = 0; col < ncols; col++)
521 DATA(map_velocity, row, col) = cell[col];
522 Rast_get_c_row(mois_fd, cell, row);
523 for (col = 0; col < ncols; col++)
524 DATA(map_mois, row, col) = cell[col];
525 }
526 }
527 G_percent(row, nrows, 2);
528
529
530 /* Scan the START layer searching for starting points.
531 * Create an array of starting points (min_heap) ordered by costs.
532 */
533
534 start_fd = Rast_open_old(start_layer, G_find_raster2(start_layer, ""));
535
536 Rast_read_range(start_layer, G_find_file("cell", start_layer, ""),
537 &start_range);
538 Rast_get_range_min_max(&start_range, &start_range_min, &start_range_max);
539
540 start_is_time = flag.start_is_time->answer;
541 /* values higher than init_time are unexpected and may cause segfaults */
542 if (start_is_time && start_range_max > init_time)
543 G_fatal_error(_("Maximum of start raster map is grater than init_time"
544 " (%d > %d)"), start_range_max, init_time);
545 /* values lower then zero does not make sense for time */
546 if (start_is_time && start_range_min < 0)
547 G_fatal_error(_("Minimum of start raster map is less than zero"
548 " (%d < 0)"), start_range_min);
549
550 /* Initialize the heap */
551 heap =
552 (struct costHa *)G_calloc(nrows * ncols + 1, sizeof(struct costHa));
553 heap_len = 0;
554
555 G_message(_("Reading %s..."), start_layer);
556 G_debug(1, "Collecting origins...");
557 collect_ori(start_fd, start_is_time);
558 G_debug(1, "Done");
559
560
561 /* Major computation of spread time */
562 G_debug(1, "Spreading...");
563 spread();
564 G_debug(1, "Done");
565
566
567 /* Open cumulative cost layer (and x, y direction layers) for writing */
568
569 cum_fd = Rast_open_c_new(out_layer);
570 if (x_out)
571 x_fd = Rast_open_c_new(x_out_layer);
572 if (y_out)
573 y_fd = Rast_open_c_new(y_out_layer);
574
575 /* prepare output -- adjust from cm to m */
576 window.ew_res = window.ew_res / 100;
577 window.ns_res = window.ns_res / 100;
578
579 /* copy maps in ram to output maps */
580 ram2out();
581
582 G_free(map_max);
583 G_free(map_dir);
584 G_free(map_base);
585 G_free(map_out);
586 G_free(map_visit);
587 if (x_out)
588 G_free(map_x_out);
589 if (y_out)
590 G_free(map_y_out);
591 if (spotting) {
592 G_free(map_spotdist);
593 G_free(map_mois);
594 G_free(map_velocity);
595 }
596
597 Rast_close(max_fd);
598 Rast_close(dir_fd);
599 Rast_close(base_fd);
600 Rast_close(start_fd);
601 Rast_close(cum_fd);
602 if (x_out)
603 Rast_close(x_fd);
604 if (y_out)
605 Rast_close(y_fd);
606 if (spotting) {
607 Rast_close(spotdist_fd);
608 Rast_close(velocity_fd);
609 Rast_close(mois_fd);
610 }
611
612 /* close graphics */
613 #if 0
614 if (display)
615 display_close();
616 #endif
617
618 exit(EXIT_SUCCESS);
619 }
620