1
2 /****************************************************************************
3 *
4 * MODULE: r.ros
5 * AUTHOR(S): Jianping Xu, Rutgers University, 1993
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: Generates rate of spread raster map layers for wildfire modeling
10 *
11 * TODO: (re)move following documentation
12 *
13 * This raster module creates three raster map layers:
14 * 1. Base (perpendicular) rate of spread (ROS);
15 * 2. Maximum (forward) ROS; and
16 * 3. Direction of the Maximum ROS.
17 * The calculation of the two ROS values for each raster
18 * cell is based on the Fortran code by Pat Andrews (1983)
19 * of the Northern Forest Fire Laboratory, USDA Forest
20 * Service. These three raster map layers are expected to
21 * be the inputs for a separate GRASS raster module
22 * 'r.spread'.
23 *
24 * 'r.ros' can be run in two standard GRASS modes:
25 * interactive and command line. For an interactive run,
26 * type in
27 * r.ros
28 * and follow the prompts; for a command line mode,
29 * type in
30 * r.ros [-v] model=name [moisture_1h=name]
31 * [moisture_10h=name]
32 * [moisture_100h=name]
33 * moisture_live=name [velocity=name]
34 * [direction=name] [slope=name]
35 * [aspect=name] output=name
36 * where:
37 * Flag:
38 * Raster Maps:
39 * model 1-13: the standard fuel models,
40 * all other numbers: same as barriers;
41 * moisture_1h 100*moisture_content;
42 * moisture_10h 100*moisture_content;
43 * moisture_100h 100*moisture_content;
44 * moisture_live 100*moisture_content;
45 * velocity ft/minute;
46 * direction degree;
47 * slope degree;
48 * aspect degree starting from East, anti-clockwise;
49 * output
50 * for Base ROS cm/min (technically not ft/min);
51 * for Max ROS cm/min (technically not ft/min);
52 * for Direction degree.
53 *
54 * Note that the name given as output will be used as
55 * PREFIX for several actual raster maps. For example, if
56 * my_ros is given as the output name, 'r.ros' will
57 * actually produce three raster maps named my_ros.base,
58 * my_ros.max, and my_ros.maxdir, or even my_ros.spotdist
59 * respectively.
60 *
61 * COPYRIGHT: (C) 2000-2009 by the GRASS Development Team
62 *
63 * This program is free software under the GNU General Public
64 * License (>=v2). Read the file COPYING that comes with GRASS
65 * for details.
66 *
67 *****************************************************************************/
68
69 #include <stdlib.h>
70 #include <stdio.h>
71 #include <math.h>
72 #include <grass/gis.h>
73 #include <grass/raster.h>
74 #include <grass/glocale.h>
75 #include "local_proto.h"
76
77
78 #define DATA(map, r, c) (map)[(r) * ncols + (c)]
79
80 /*measurements of the 13 fuel models, input of Rothermel equation (1972) */
81 float WO[4][14] =
82 { {0, 0.034, 0.092, 0.138, 0.230, 0.046, 0.069, 0.052, 0.069, 0.134,
83 0.138, 0.069, 0.184, 0.322},
84 {0, 0, 0.046, 0, 0.184, 0.023, 0.115, 0.086, 0.046, 0.019, 0.092, 0.207,
85 0.644, 1.058},
86 {0, 0, 0.023, 0, 0.092, 0, 0.092, 0.069, 0.115, 0.007, 0.230, 0.253, 0.759,
87 1.288},
88 {0, 0, 0.023, 0, 0.230, 0.092, 0, 0.017, 0, 0, 0.092, 0, 0}
89 };
90
91 /*ovendry fuel loading, lb./ft.^2 */
92 float DELTA[] = { 0, 1.0, 1.0, 2.5, 6.0, 2.0, 2.5, 2.5,
93 0.2, 0.2, 1.0, 1.0, 2.3, 3.0
94 }; /*fuel depth, ft. */
95 float SIGMA[4][14] =
96 { {0, 3500, 3000, 1500, 2000, 2000, 1750, 1750, 2000, 2500, 2000, 1500,
97 1500, 1500},
98 {0, 0, 109, 0, 109, 109, 109, 109, 109, 109, 109, 109, 109, 109},
99 {0, 0, 30, 0, 30, 0, 30, 30, 30, 30, 30, 30, 30, 30},
100 {0, 0, 1500, 0, 1500, 1500, 0, 1500, 0, 0, 1500, 0, 0, 0}
101 };
102
103 /*fuel particale surface-area-to-volume ratio, 1/ft. */
104 float MX[] = { 0, 0.12, 0.15, 0.25, 0.20, 0.20, 0.25, 0.40,
105 0.30, 0.25, 0.25, 0.15, 0.20, 0.25
106 }; /*moisture content of extinction */
107
108 CELL *map_elev; /*full array for elevation map layer (for spotting) */
109 int nrows, ncols;
110 struct Cell_head window;
111
112
main(int argc,char * argv[])113 int main(int argc, char *argv[])
114 {
115
116 /***input of Rothermel equation (1972)***/
117 float h = 8000.0, /*heat of combustion, BTU/lb. */
118 rhop = 32, /*ovendry fuel density, lb./ft.^3 */
119 ST = 0.0555; /*fuel total mineral content, lb. minerals/lb. ovendry */
120
121 float sigma[14];
122
123 /***derived parameters of Rothermel equation (1972)***/
124 float R, /*rate of spread, ft./min.
125 R = IR*xi*(1+phiw+phis)/(rhob*epsilon*Qig) */
126 IR, /*reaction intensity, BTU/ft.^2/min.
127 IR = gamma*wn*h*etaM*etas */
128 gamma, /*optimum reation velosity, 1/min.
129 gamma = gammamax*(beta/betaop)^A*
130 exp(A(1-beta/betaop)) */
131 gammamax, /*maximum reation velosity, 1/min.
132 gammamax = sigma^1.5/(495+0.0594*sigma^1.5) */
133 betaop, /*optimum packing ratio
134 betaop = 3.348/(sigma^0.8189) */
135 A, /*A = 1/(4.77*sigma^0.1-7.27) */
136 etas = 0.174 / pow(0.01, 0.19), /*mineral damping coefficient */
137 xi, /*propagating flux ratio,
138 xi = exp((0.792+0.681*sigma^0.5)(beta+0.1))/
139 (192+0.2595*sigma) */
140 phiw, /*wind coefficient, phiw = C*U^B/(beta/betaop)^E */
141 C, /*C = 7.47/exp(0.133*sigma^0.55) */
142 B, /*B = 0.02526*sigma^.054 */
143 E, /*E = 0.715/exp(0.000359*sigma) */
144 phis, /*slope coefficient,phis = 5.275/beta^0.3*(tan(theta)^2) */
145 rhob, /*ovendry bulk density, lb./ft.^3, rohb = wo/delta */
146 epsilon[4][14], /*effective heating number, epsilon = 1/exp(138/sigma) */
147 Qig[14], /*heat of preignition, BTU/lb. Qig = 250+1116*Mf */
148 beta; /*packing ratio, beta = rhob/rhop */
149
150 /***intermediate variables***/
151 float R0, /*base ROS (w/out wind/slope) */
152 Rdir, sin_fac, cos_fac, Ffactor_all[4][14], /*in all fuel subclasses by sigma/WO */
153 Ffactor_in_dead[3][14], /*in dead fuel subclasses by sigma/WO */
154 Ffactor_dead[14], /*dead fuel weight by sigma/WO */
155 Ffactor_live[14], /*live fuel weight by sigma/WO */
156 Gfactor_in_dead[3][14], /*in dead fuel by the 6 classes */
157 G1, G2, G3, G4, G5, wo_dead[14], /*dead fuel total load */
158 wn_dead, /*net dead fuel total load */
159 wn_live, /*net live fuel (total) load */
160 class_sum, moisture[4], /*moistures of 1-h,10-h,100-h,live fuels */
161 Mf_dead, /*total moisture of dead fuels */
162 etaM_dead, /*dead fuel misture damping coefficient */
163 etaM_live, /*live fuel misture damping coefficient */
164 xmext, /*live fuel moisture of extinction */
165 phi_ws, /*wind and slope conbined coefficient */
166 wmfd, fdmois, fined, finel;
167
168 /*other local variables */
169 int col, row,
170 spotting,
171 model, class,
172 fuel_fd = 0,
173 mois_1h_fd = 0, mois_10h_fd = 0, mois_100h_fd = 0, mois_live_fd = 0,
174 vel_fd = 0, dir_fd = 0,
175 elev_fd = 0, slope_fd = 0, aspect_fd = 0,
176 base_fd = 0, max_fd = 0, maxdir_fd = 0, spotdist_fd = 0;
177
178 char *name_base, *name_max, *name_maxdir, *name_spotdist;
179
180 CELL *fuel, /*cell buffer for fuel model map layer */
181 *mois_1h, /*cell buffer for 1-hour fuel moisture map layer */
182 *mois_10h, /*cell buffer for 10-hour fuel moisture map layer */
183 *mois_100h, /*cell buffer for 100-hour fuel moisture map layer */
184 *mois_live, /*cell buffer for live fuel moisture map layer */
185 *vel, /*cell buffer for wind velocity map layer */
186 *dir, /*cell buffer for wind direction map layer */
187 *elev = NULL, /*cell buffer for elevation map layer (for spotting) */
188 *slope, /*cell buffer for slope map layer */
189 *aspect, /*cell buffer for aspect map layer */
190 *base, /*cell buffer for base ROS map layer */
191 *max, /*cell buffer for max ROS map layer */
192 *maxdir, /*cell buffer for max ROS direction map layer */
193 *spotdist = NULL; /*cell buffer for max spotting distance map layer */
194
195 extern struct Cell_head window;
196
197 struct
198 {
199 struct Option *model,
200 *mois_1h, *mois_10h, *mois_100h, *mois_live,
201 *vel, *dir, *elev, *slope, *aspect, *base, *max, *maxdir, *spotdist;
202 } parm;
203
204 /* please, remove before GRASS 7 released */
205 struct GModule *module;
206
207 /* initialize access to database and create temporary files */
208 G_gisinit(argv[0]);
209
210 /* Set description */
211 module = G_define_module();
212 G_add_keyword(_("raster"));
213 G_add_keyword(_("fire"));
214 G_add_keyword(_("spread"));
215 G_add_keyword(_("rate of spread"));
216 G_add_keyword(_("hazard"));
217 G_add_keyword(_("model"));
218 module->label = _("Generates rate of spread raster maps.");
219 module->description =
220 _("Generates three, or four raster map layers showing the base "
221 "(perpendicular) rate of spread (ROS), the maximum (forward) ROS, "
222 "the direction of the maximum ROS, and optionally the "
223 "maximum potential spotting distance for fire spread simulation.");
224
225 parm.model = G_define_standard_option(G_OPT_R_INPUT);
226 parm.model->key = "model";
227 parm.model->label = _("Raster map containing fuel models");
228 parm.model->description =
229 _("Name of an "
230 "existing raster map layer in the user's current mapset search path containing "
231 "the standard fuel models defined by the USDA Forest Service. Valid values "
232 "are 1-13; other numbers are recognized as barriers by r.ros.");
233
234 parm.mois_1h = G_define_standard_option(G_OPT_R_INPUT);
235 parm.mois_1h->key = "moisture_1h";
236 parm.mois_1h->required = NO;
237 parm.mois_1h->label =
238 _("Raster map containing the 1-hour fuel moisture (%)");
239 parm.mois_1h->description =
240 _("Name of an existing raster map layer in "
241 "the user's current mapset search path containing the 1-hour (<.25\") "
242 "fuel moisture (percentage content multiplied by 100).");
243
244 parm.mois_10h = G_define_standard_option(G_OPT_R_INPUT);
245 parm.mois_10h->key = "moisture_10h";
246 parm.mois_10h->required = NO;
247 parm.mois_10h->label =
248 _("Raster map containing the 10-hour fuel moisture (%)");
249 parm.mois_10h->description =
250 _("Name of an existing raster map layer in the "
251 "user's current mapset search path containing the 10-hour (.25-1\") fuel "
252 "moisture (percentage content multiplied by 100).");
253
254 parm.mois_100h = G_define_standard_option(G_OPT_R_INPUT);
255 parm.mois_100h->key = "moisture_100h";
256 parm.mois_100h->required = NO;
257 parm.mois_100h->label =
258 _("Raster map containing the 100-hour fuel moisture (%)");
259 parm.mois_100h->description =
260 _("Name of an existing raster map layer in the "
261 "user's current mapset search path containing the 100-hour (1-3\") fuel moisture "
262 "(percentage content multiplied by 100).");
263
264 parm.mois_live = G_define_standard_option(G_OPT_R_INPUT);
265 parm.mois_live->key = "moisture_live";
266 parm.mois_live->label =
267 _("Raster map containing live fuel moisture (%)");
268 parm.mois_live->description =
269 _("Name of an existing raster map layer in the "
270 "user's current mapset search path containing live (herbaceous) fuel "
271 "moisture (percentage content multiplied by 100).");
272
273 parm.vel = G_define_standard_option(G_OPT_R_INPUT);
274 parm.vel->key = "velocity";
275 parm.vel->required = NO;
276 parm.vel->label =
277 _("Raster map containing midflame wind velocities (ft/min)");
278 parm.vel->description =
279 _("Name of an existing raster map layer in the user's "
280 "current mapset search path containing wind velocities at half of the average "
281 "flame height (feet/minute).");
282
283 parm.dir = G_define_standard_option(G_OPT_R_INPUT);
284 parm.dir->key = "direction";
285 parm.dir->required = NO;
286 parm.dir->label =
287 _("Name of raster map containing wind directions (degree)");
288 parm.dir->description =
289 _("Name of an existing raster map "
290 "layer in the user's current mapset search path containing wind direction, "
291 "clockwise from north (degree).");
292
293 parm.slope = G_define_standard_option(G_OPT_R_INPUT);
294 parm.slope->key = "slope";
295 parm.slope->required = NO;
296 parm.slope->label =
297 _("Name of raster map containing slope (degree)");
298 parm.slope->description =
299 _("Name of an existing raster map layer "
300 "in the user's current mapset search path containing "
301 "topographic slope (degree).");
302
303 parm.aspect = G_define_standard_option(G_OPT_R_INPUT);
304 parm.aspect->key = "aspect";
305 parm.aspect->required = NO;
306 parm.aspect->label =
307 _("Raster map containing aspect (degree, CCW from E)");
308 parm.aspect->description =
309 _("Name of an existing "
310 "raster map layer in the user's current mapset search path containing "
311 "topographic aspect, counterclockwise from east (GRASS convention) "
312 "in degrees.");
313
314 parm.elev = G_define_standard_option(G_OPT_R_ELEV);
315 parm.elev->required = NO;
316 parm.elev->label =
317 _("Raster map containing elevation (m, required for spotting)");
318 parm.elev->description =
319 _("Name of an existing raster map "
320 "layer in the user's current mapset search path containing elevation (meters). "
321 "Option is required from spotting distance computation "
322 "(when spotting_distance option is provided)");
323
324 parm.base = G_define_standard_option(G_OPT_R_OUTPUT);
325 parm.base->key = "base_ros";
326 parm.base->required = YES;
327 parm.base->label =
328 _("Output raster map containing base ROS (cm/min)");
329 parm.base->description =
330 _("Base (perpendicular) rate of spread (ROS)");
331
332 parm.max = G_define_standard_option(G_OPT_R_OUTPUT);
333 parm.max->key = "max_ros";
334 parm.max->required = YES;
335 parm.max->label =
336 _("Output raster map containing maximal ROS (cm/min)");
337 parm.max->description =
338 _("The maximum (forward) rate of spread (ROS)");
339
340 parm.maxdir = G_define_standard_option(G_OPT_R_OUTPUT);
341 parm.maxdir->key = "direction_ros";
342 parm.maxdir->required = YES;
343 parm.maxdir->label =
344 _("Output raster map containing directions of maximal ROS (degree)");
345 parm.maxdir->description =
346 _("The direction of the maximal (forward) rate of spread (ROS)");
347
348 parm.spotdist = G_define_standard_option(G_OPT_R_OUTPUT);
349 parm.spotdist->key = "spotting_distance";
350 parm.spotdist->required = NO;
351 parm.spotdist->label =
352 _("Output raster map containing maximal spotting distance (m)");
353 parm.spotdist->description =
354 _("The maximal potential spotting distance"
355 " (requires elevation raster map to be provided).");
356
357 /* Parse command line */
358 if (G_parser(argc, argv))
359 exit(EXIT_FAILURE);
360
361 if (parm.spotdist->answer)
362 spotting = 1;
363 else
364 spotting = 0;
365
366 /* Check if input layers exists in data base */
367 if (G_find_raster2(parm.model->answer, "") == NULL)
368 G_fatal_error(_("Raster map <%s> not found"), parm.model->answer);
369
370 if (!
371 (parm.mois_1h->answer || parm.mois_10h->answer ||
372 parm.mois_100h->answer)) {
373 G_fatal_error(_("No dead fuel moisture is given. "
374 "At least one of the 1-h, 10-h, 100-h moisture layers is required."));
375 }
376
377 if (parm.mois_1h->answer) {
378 if (G_find_raster2(parm.mois_1h->answer, "") == NULL)
379 G_fatal_error(_("Raster map <%s> not found"),
380 parm.mois_1h->answer);
381 }
382 if (parm.mois_10h->answer) {
383 if (G_find_raster2(parm.mois_10h->answer, "") == NULL)
384 G_fatal_error(_("Raster map <%s> not found"),
385 parm.mois_10h->answer);
386 }
387 if (parm.mois_100h->answer) {
388 if (G_find_raster2(parm.mois_100h->answer, "") == NULL)
389 G_fatal_error(_("Raster map <%s> not found"),
390 parm.mois_100h->answer);
391 }
392
393 if (G_find_raster2(parm.mois_live->answer, "") == NULL)
394 G_fatal_error(_("Raster map <%s> not found"), parm.mois_live->answer);
395
396 if (parm.vel->answer && !(parm.dir->answer)) {
397 G_fatal_error(_("A wind direction layer should be "
398 "given if the wind velocity layer <%s> has been given"),
399 parm.vel->answer);
400 }
401 if (!(parm.vel->answer) && parm.dir->answer) {
402 G_fatal_error(_("A wind velocity layer should be given "
403 "if the wind direction layer <%s> has been given"),
404 parm.dir->answer);
405 }
406 if (parm.vel->answer) {
407 if (G_find_raster2(parm.vel->answer, "") == NULL)
408 G_fatal_error(_("Raster map <%s> not found"), parm.vel->answer);
409 }
410 if (parm.dir->answer) {
411 if (G_find_raster2(parm.dir->answer, "") == NULL)
412 G_fatal_error(_("Raster map <%s> not found"), parm.dir->answer);
413 }
414
415 if (parm.slope->answer && !(parm.aspect->answer)) {
416 G_fatal_error(_("An aspect layer should be given "
417 "if the slope layer <%s> has been given"),
418 parm.slope->answer);
419 }
420 if (!(parm.slope->answer) && parm.aspect->answer) {
421 G_fatal_error(_("A slope layer should be given "
422 "if the aspect layer <%s> has been given"),
423 parm.aspect->answer);
424 }
425 if (parm.slope->answer) {
426 if (G_find_raster2(parm.slope->answer, "") == NULL)
427 G_fatal_error(_("Raster map <%s> not found"), parm.slope->answer);
428 }
429 if (parm.aspect->answer) {
430 if (G_find_raster2(parm.aspect->answer, "") == NULL)
431 G_fatal_error(_("Raster map <%s> not found"),
432 parm.aspect->answer);
433 }
434
435 if (spotting) {
436 if (!(parm.elev->answer)) {
437 G_fatal_error(_("An elevation layer should be given "
438 "if considering spotting"));
439 }
440 else {
441 if (G_find_raster2(parm.elev->answer, "") == NULL)
442 G_fatal_error(_("Raster map <%s> not found"),
443 parm.elev->answer);
444 }
445 }
446
447 /*assign names of the three output ROS layers */
448 name_base = parm.base->answer;
449 name_max = parm.max->answer;
450 name_maxdir = parm.maxdir->answer;
451
452 /*assign a name to output SPOTTING distance layer */
453 if (spotting) {
454 name_spotdist = parm.spotdist->answer;
455 }
456
457 /* Get database window parameters */
458 G_get_window(&window);
459
460 /* find number of rows and columns in window */
461 nrows = Rast_window_rows();
462 ncols = Rast_window_cols();
463
464 fuel = Rast_allocate_c_buf();
465 mois_1h = Rast_allocate_c_buf();
466 mois_10h = Rast_allocate_c_buf();
467 mois_100h = Rast_allocate_c_buf();
468 mois_live = Rast_allocate_c_buf();
469 vel = Rast_allocate_c_buf();
470 dir = Rast_allocate_c_buf();
471 slope = Rast_allocate_c_buf();
472 aspect = Rast_allocate_c_buf();
473 base = Rast_allocate_c_buf();
474 max = Rast_allocate_c_buf();
475 maxdir = Rast_allocate_c_buf();
476 if (spotting) {
477 spotdist = Rast_allocate_c_buf();
478 elev = Rast_allocate_c_buf();
479 map_elev = (CELL *) G_calloc(nrows * ncols, sizeof(CELL));
480 }
481
482 /* Open input cell layers for reading */
483
484 fuel_fd = Rast_open_old(parm.model->answer, "");
485
486 if (parm.mois_1h->answer)
487 mois_1h_fd = Rast_open_old(parm.mois_1h->answer, "");
488
489 if (parm.mois_10h->answer)
490 mois_10h_fd = Rast_open_old(parm.mois_10h->answer,"");
491
492 if (parm.mois_100h->answer)
493 mois_100h_fd = Rast_open_old(parm.mois_100h->answer, "");
494
495 mois_live_fd = Rast_open_old(parm.mois_live->answer, "");
496
497 if (parm.vel->answer)
498 vel_fd = Rast_open_old(parm.vel->answer, "");
499
500 if (parm.dir->answer)
501 dir_fd = Rast_open_old(parm.dir->answer, "");
502
503 if (parm.slope->answer)
504 slope_fd = Rast_open_old(parm.slope->answer, "");
505
506 if (parm.aspect->answer)
507 aspect_fd = Rast_open_old(parm.aspect->answer, "");
508
509 if (spotting)
510 elev_fd = Rast_open_old(parm.elev->answer, "");
511
512 base_fd = Rast_open_c_new(name_base);
513 max_fd = Rast_open_c_new(name_max);
514 maxdir_fd = Rast_open_c_new(name_maxdir);
515 if (spotting)
516 spotdist_fd = Rast_open_c_new(name_spotdist);
517
518 /*compute weights, combined wo, and combined sigma */
519 /*wo[model] -- simple sum of WO[class][model] by all fuel subCLASS */
520 /*sigma[model] -- weighted sum of SIGMA[class][model] by all fuel subCLASS *epsilon[class][model] */
521 for (model = 1; model <= 13; model++) {
522 class_sum = 0.0;
523 wo_dead[model] = 0.0;
524 sigma[model] = 0.0;
525 for (class = 0; class <= 3; class++) {
526 class_sum = class_sum + WO[class][model] * SIGMA[class][model];
527 if (SIGMA[class][model] > 0.0) {
528 epsilon[class][model] = exp(-138.0 / SIGMA[class][model]);
529 }
530 else {
531 epsilon[class][model] = 0.0;
532 }
533 }
534 for (class = 0; class <= 3; class++) {
535 Ffactor_all[class][model] =
536 WO[class][model] * SIGMA[class][model] / class_sum;
537 sigma[model] =
538 sigma[model] +
539 SIGMA[class][model] * Ffactor_all[class][model];
540 }
541 class_sum = 0.0;
542 for (class = 0; class <= 2; class++) {
543 wo_dead[model] = wo_dead[model] + WO[class][model];
544 class_sum = class_sum + WO[class][model] * SIGMA[class][model];
545 }
546 for (class = 0; class <= 2; class++) {
547 Ffactor_in_dead[class][model] =
548 WO[class][model] * SIGMA[class][model] / class_sum;
549 }
550
551 /* compute G factor for each of the 6 subclasses */
552 G1 = 0.0;
553 G2 = 0.0;
554 G3 = 0.0;
555 G4 = 0.0;
556 G5 = 0.0;
557 for (class = 0; class <= 2; class++) {
558 if (SIGMA[class][model] >= 1200)
559 G1 = G1 + Ffactor_in_dead[class][model];
560 if (SIGMA[class][model] < 1200 && SIGMA[class][model] >= 192)
561 G2 = G2 + Ffactor_in_dead[class][model];
562 if (SIGMA[class][model] < 192 && SIGMA[class][model] >= 96)
563 G3 = G3 + Ffactor_in_dead[class][model];
564 if (SIGMA[class][model] < 96 && SIGMA[class][model] >= 48)
565 G4 = G4 + Ffactor_in_dead[class][model];
566 if (SIGMA[class][model] < 48 && SIGMA[class][model] >= 16)
567 G5 = G5 + Ffactor_in_dead[class][model];
568 }
569 for (class = 0; class <= 2; class++) {
570 if (SIGMA[class][model] >= 1200)
571 Gfactor_in_dead[class][model] = G1;
572 if (SIGMA[class][model] < 1200 && SIGMA[class][model] >= 192)
573 Gfactor_in_dead[class][model] = G2;
574 if (SIGMA[class][model] < 192 && SIGMA[class][model] >= 96)
575 Gfactor_in_dead[class][model] = G3;
576 if (SIGMA[class][model] < 96 && SIGMA[class][model] >= 48)
577 Gfactor_in_dead[class][model] = G4;
578 if (SIGMA[class][model] < 48 && SIGMA[class][model] >= 16)
579 Gfactor_in_dead[class][model] = G5;
580 if (SIGMA[class][model] < 16)
581 Gfactor_in_dead[class][model] = 0.0;
582 }
583
584 Ffactor_dead[model] =
585 class_sum / (class_sum + WO[3][model] * SIGMA[3][model]);
586 Ffactor_live[model] = 1 - Ffactor_dead[model];
587 }
588
589 /*if considering spotting, read elevation map into an array */
590 if (spotting)
591 for (row = 0; row < nrows; row++) {
592 Rast_get_c_row(elev_fd, elev, row);
593 for (col = 0; col < ncols; col++)
594 DATA(map_elev, row, col) = elev[col];
595 }
596
597 /*major computation: compute ROSs one cell a time */
598 for (row = 0; row < nrows; row++) {
599 G_percent(row, nrows, 2);
600 Rast_get_c_row(fuel_fd, fuel, row);
601 if (parm.mois_1h->answer)
602 Rast_get_c_row(mois_1h_fd, mois_1h, row);
603 if (parm.mois_10h->answer)
604 Rast_get_c_row(mois_10h_fd, mois_10h, row);
605 if (parm.mois_100h->answer)
606 Rast_get_c_row(mois_100h_fd, mois_100h, row);
607 Rast_get_c_row(mois_live_fd, mois_live, row);
608 if (parm.vel->answer)
609 Rast_get_c_row(vel_fd, vel, row);
610 if (parm.dir->answer)
611 Rast_get_c_row(dir_fd, dir, row);
612 if (parm.slope->answer)
613 Rast_get_c_row(slope_fd, slope, row);
614 if (parm.aspect->answer)
615 Rast_get_c_row(aspect_fd, aspect, row);
616
617 /*initialize cell buffers for output map layers */
618 for (col = 0; col < ncols; col++) {
619 base[col] = max[col] = maxdir[col] = 0;
620 if (spotting)
621 spotdist[col] = 0;
622 }
623
624 for (col = 0; col < ncols; col++) {
625 /*check if a fuel is within the 13 models,
626 *if not, no processing; useful when no data presents*/
627 if (fuel[col] < 1 || fuel[col] > 13)
628 continue;
629 if (parm.mois_1h->answer)
630 moisture[0] = 0.01 * mois_1h[col];
631 if (parm.mois_10h->answer)
632 moisture[1] = 0.01 * mois_10h[col];
633 if (parm.mois_100h->answer)
634 moisture[2] = 0.01 * mois_100h[col];
635 moisture[3] = 0.01 * mois_live[col];
636 if (parm.aspect->answer)
637 aspect[col] = (630 - aspect[col]) % 360;
638
639 /* assign some dead fuel moisture if not completely
640 * given based on empirical relationship*/
641 if (!(parm.mois_10h->answer || parm.mois_100h->answer)) {
642 moisture[1] = moisture[0] + 0.01;
643 moisture[2] = moisture[0] + 0.02;
644 }
645 if (!(parm.mois_1h->answer || parm.mois_100h->answer)) {
646 moisture[0] = moisture[1] - 0.01;
647 moisture[2] = moisture[1] + 0.01;
648 }
649 if (!(parm.mois_1h->answer || parm.mois_10h->answer)) {
650 moisture[0] = moisture[2] - 0.02;
651 moisture[1] = moisture[2] - 0.01;
652 }
653 if (!(parm.mois_1h->answer) && parm.mois_10h->answer &&
654 parm.mois_100h->answer)
655 moisture[0] = moisture[1] - 0.01;
656 if (!(parm.mois_10h->answer) && parm.mois_1h->answer &&
657 parm.mois_100h->answer)
658 moisture[1] = moisture[0] + 0.01;
659 if (!(parm.mois_100h->answer) && parm.mois_1h->answer &&
660 parm.mois_10h->answer)
661 moisture[2] = moisture[1] + 0.01;
662
663 /*compute xmext, moisture of extinction of live fuels */
664 wmfd = 0.0;
665 fined = 0.0;
666 if (SIGMA[3][fuel[col]] > 0.0) {
667 for (class = 0; class <= 2; class++) {
668 if (SIGMA[class][fuel[col]] == 0.0)
669 continue;
670 fined =
671 fined +
672 WO[class][fuel[col]] * exp(-138.0 /
673 SIGMA[class][fuel[col]]);
674 wmfd =
675 wmfd +
676 WO[class][fuel[col]] * exp(-138.0 /
677 SIGMA[class][fuel[col]]) *
678 moisture[class];
679 }
680 fdmois = wmfd / fined;
681 finel = WO[3][fuel[col]] * exp(-500.0 / SIGMA[3][fuel[col]]);
682 xmext =
683 2.9 * (fined / finel) * (1 - fdmois / MX[fuel[col]]) -
684 0.226;
685 }
686 else
687 xmext = MX[fuel[col]];
688 if (xmext < MX[fuel[col]])
689 xmext = MX[fuel[col]];
690
691 /*compute other intermediate values */
692 Mf_dead = 0.0;
693 wn_dead = 0.0;
694 class_sum = 0.0;
695 for (class = 0; class <= 2; class++) {
696 Mf_dead =
697 Mf_dead +
698 moisture[class] * Ffactor_in_dead[class][fuel[col]];
699 wn_dead =
700 wn_dead +
701 WO[class][fuel[col]] * Gfactor_in_dead[class][fuel[col]] *
702 (1 - ST);
703 Qig[class] = 250 + 1116 * moisture[class];
704 class_sum =
705 class_sum +
706 Ffactor_all[class][fuel[col]] *
707 epsilon[class][fuel[col]] * Qig[class];
708 }
709 etaM_dead =
710 1.0 - 2.59 * (Mf_dead / MX[fuel[col]]) +
711 5.11 * (Mf_dead / MX[fuel[col]]) * (Mf_dead / MX[fuel[col]]) -
712 3.52 * (Mf_dead / MX[fuel[col]]) * (Mf_dead / MX[fuel[col]]) *
713 (Mf_dead / MX[fuel[col]]);
714 if (Mf_dead >= MX[fuel[col]])
715 etaM_dead = 0.0;
716 etaM_live =
717 1.0 - 2.59 * (moisture[3] / xmext) +
718 5.11 * (moisture[3] / xmext) * (moisture[3] / xmext) -
719 3.52 * (moisture[3] / xmext) * (moisture[3] / xmext) *
720 (moisture[3] / xmext);
721 if (moisture[3] >= xmext)
722 etaM_live = 0.0;
723 wn_live = WO[3][fuel[col]] * (1 - ST);
724 Qig[3] = 250 + 1116 * moisture[3];
725 class_sum =
726 class_sum +
727 Ffactor_all[3][fuel[col]] * epsilon[3][fuel[col]] * Qig[3];
728
729 /*final computations */
730 rhob = (wo_dead[fuel[col]] + WO[3][fuel[col]]) / DELTA[fuel[col]];
731 beta = rhob / rhop;
732 betaop = 3.348 / pow(sigma[fuel[col]], 0.8189);
733 A = 133 / pow(sigma[fuel[col]], 0.7913);
734 gammamax =
735 pow(sigma[fuel[col]],
736 1.5) / (495 + 0.0594 * pow(sigma[fuel[col]], 1.5));
737 gamma =
738 gammamax * pow(beta / betaop,
739 A) * exp(A * (1 - beta / betaop));
740 xi = exp((0.792 + 0.681 * pow(sigma[fuel[col]], 0.5)) * (beta +
741 0.1)) /
742 (192 + 0.2595 * sigma[fuel[col]]);
743 IR = gamma * h * (wn_dead * etaM_dead +
744 wn_live * etaM_live) * etas;
745
746 R0 = IR * xi / (rhob * class_sum);
747
748 if (parm.vel->answer && parm.dir->answer) {
749 C = 7.47 * exp(-0.133 * pow(sigma[fuel[col]], 0.55));
750 B = 0.02526 * pow(sigma[fuel[col]], 0.54);
751 E = 0.715 * exp(-0.000359 * sigma[fuel[col]]);
752 phiw = C * pow((double)vel[col], B) / pow(beta / betaop, E);
753 }
754 else
755 phiw = 0.0;
756
757 if (parm.slope->answer && parm.aspect->answer) {
758 phis =
759 5.275 * pow(beta,
760 -0.3) * tan(M_D2R * slope[col]) *
761 tan(M_D2R * slope[col]);
762 }
763 else
764 phis = 0.0;
765
766 /*compute the maximum ROS R and its direction,
767 *vector adding for the effect of wind and slope*/
768 if (parm.dir->answer && parm.aspect->answer) {
769 sin_fac =
770 phiw * sin(M_D2R * dir[col]) +
771 phis * sin(M_D2R * aspect[col]);
772 cos_fac =
773 phiw * cos(M_D2R * dir[col]) +
774 phis * cos(M_D2R * aspect[col]);
775 phi_ws = sqrt(sin_fac * sin_fac + cos_fac * cos_fac);
776 Rdir = atan2(sin_fac, cos_fac) / M_D2R;
777 }
778 else if (parm.dir->answer && !(parm.aspect->answer)) {
779 phi_ws = phiw;
780 Rdir = dir[col];
781 }
782 else if (!(parm.dir->answer) && parm.aspect->answer) {
783 phi_ws = phis;
784 Rdir = (float)aspect[col];
785 }
786 else {
787 phi_ws = 0.0;
788 Rdir = 0.0;
789 }
790 R = R0 * (1 + phi_ws);
791 if (Rdir < 0.0)
792 Rdir = Rdir + 360;
793 /*printf("\nparm.dir->aanswer=%s, parm.aspect->aanswer=%s, phis=%f, phi_ws=%f, aspect[col]=%d,Rdir=%f",parm.dir->answer,parm.aspect->answer,phis,phi_ws,aspect[col],Rdir); */
794
795 /*maximum spotting distance */
796 if (spotting)
797 spotdist[col] =
798 spot_dist(fuel[col], R, vel[col], Rdir, row, col);
799
800 /*to avoid 0 R, convert ft./min to cm/min */
801 R0 = 30.5 * R0;
802 R = 30.5 * R;
803 /*4debugging R0 = R0/30.5/1.1; R = R/30.5/1.1; */
804
805 base[col] = (int)R0;
806 max[col] = (int)R;
807 maxdir[col] = (int)Rdir;
808 /*printf("(%d, %d)\nR0=%.2f, vel=%d, dir=%d, phiw=%.2f, s=%d, as=%d, phis=%.2f, R=%.1f, Rdir=%.0f\n", row, col, R0, vel[col], dir[col], phiw, slope[col], aspect[col], phis, R, Rdir); */
809 }
810 Rast_put_row(base_fd, base, CELL_TYPE);
811 Rast_put_row(max_fd, max, CELL_TYPE);
812 Rast_put_row(maxdir_fd, maxdir, CELL_TYPE);
813 if (spotting)
814 Rast_put_row(spotdist_fd, spotdist, CELL_TYPE);
815 }
816 G_percent(row, nrows, 2);
817
818 Rast_close(fuel_fd);
819 if (parm.mois_1h->answer)
820 Rast_close(mois_1h_fd);
821 if (parm.mois_10h->answer)
822 Rast_close(mois_10h_fd);
823 if (parm.mois_100h->answer)
824 Rast_close(mois_100h_fd);
825 Rast_close(mois_live_fd);
826 if (parm.vel->answer)
827 Rast_close(vel_fd);
828 if (parm.dir->answer)
829 Rast_close(dir_fd);
830 if (parm.slope->answer)
831 Rast_close(slope_fd);
832 if (parm.aspect->answer)
833 Rast_close(aspect_fd);
834 Rast_close(base_fd);
835 Rast_close(max_fd);
836 Rast_close(maxdir_fd);
837 if (spotting) {
838 Rast_close(spotdist_fd);
839 G_free(map_elev);
840 }
841
842 /*
843 for (model = 1; model <= 13; model++) {
844 if (model == 1)
845 printf("\n Grass and Grass-dominated\n");
846 else if (model == 4)
847 printf(" Chaparral and Shrubfields\n");
848 else if (model == 8)
849 printf(" Timber Litter\n");
850 else if (model == 11)
851 printf(" Logging Slash\n");
852 printf("Model %2d ", model);
853 for (class = 0; class <= 3; class++)
854 printf("%4.0f/%.3f ", SIGMA[class][model], WO[class][model]);
855 printf(" %.1f %.2f\n", DELTA[model], MX[model]);
856 }
857
858 printf("\nWeight in All Fuel Subclasses:\n");
859 for (model = 1; model <= 13; model++) {
860 printf("model %2d ", model);
861 for (class = 0; class <= 3; class++)
862 printf("%.2f ", Ffactor_all[class][model]);
863 printf("%4.0f/%.3f=%6.0f model %2d\n", sigma[model], wo_dead[model]+WO[3][model], sigma[model]/(wo_dead[model]+WO[3][model]), model);
864 }
865 printf("\nWeight in Dead Fuel Subclasses, Dead Weight/Live Weight:\n");for (model = 1; model <= 13; model++) {
866 printf("model %2d ", model);
867 for (class = 0; class <= 2; class++)
868 printf("%.2f ", Ffactor_in_dead[class][model]);
869 printf("%.2f/%.2f model %2d\n", Ffactor_dead[model], Ffactor_live[model], model);
870 }
871 */
872
873 G_done_msg(_("Raster maps <%s>, <%s> and <%s> created."), name_base, name_max, name_maxdir);
874
875 exit(EXIT_SUCCESS);
876 }
877