2
/****************************************************************************
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>
11
* This raster module creates three raster map layers:
12
* 1. Base (perpendicular) rate of spread (ROS);
13
* 2. Maximum (forward) ROS; and
14
* 3. Direction of the Maximum ROS.
15
* The calculation of the two ROS values for each raster
16
* cell is based on the Fortran code by Pat Andrews (1983)
17
* of the Northern Forest Fire Laboratory, USDA Forest
18
* Service. These three raster map layers are expected to
19
* be the inputs for a seperate GRASS raster module
22
* 'r.ros' can be run in two standard GRASS modes:
23
* interactive and command line. For an interactive run,
26
* and follow the prompts; for a command line mode,
28
* r.ros [-v] model=name [moisture_1h=name]
30
* [moisture_100h=name]
31
* moisture_live=name [velocity=name]
32
* [direction=name] [slope=name]
33
* [aspect=name] output=name
37
* model 1-13: the standard fuel models,
38
* all other numbers: same as barriers;
39
* moisture_1h 100*moisture_content;
40
* moisture_10h 100*moisture_content;
41
* moisture_100h 100*moisture_content;
42
* moisture_live 100*moisture_content;
46
* aspect degree starting from East, anti-clockwise;
48
* for Base ROS cm/min (technically not ft/min);
49
* for Max ROS cm/min (technically not ft/min);
50
* for Direction degree.
52
* Note that the name given as output will be used as
53
* PREFIX for several actual raster maps. For example, if
54
* my_ros is given as the output name, 'r.ros' will
55
* actually produce three raster maps named my_ros.base,
56
* my_ros.max, and my_ros.maxdir, or even my_ros.spotdist
59
* COPYRIGHT: (C) 2000-2006 by the GRASS Development Team
61
* This program is free software under the GNU General Public
62
* License (>=v2). Read the file COPYING that comes with GRASS
65
*****************************************************************************/
70
#include <grass/gis.h>
71
#include <grass/glocale.h>
72
#include "local_proto.h"
75
#define DATA(map, r, c) (map)[(r) * ncols + (c)]
79
/*measurements of the 13 fuel models, input of Rothermel equation (1972) */
81
{ {0, 0.034, 0.092, 0.138, 0.230, 0.046, 0.069, 0.052, 0.069, 0.134,
82
0.138, 0.069, 0.184, 0.322},
83
{0, 0, 0.046, 0, 0.184, 0.023, 0.115, 0.086, 0.046, 0.019, 0.092, 0.207,
85
{0, 0, 0.023, 0, 0.092, 0, 0.092, 0.069, 0.115, 0.007, 0.230, 0.253, 0.759,
87
{0, 0, 0.023, 0, 0.230, 0.092, 0, 0.017, 0, 0, 0.092, 0, 0}
90
/*ovendry fuel loading, lb./ft.^2 */
91
float DELTA[] = { 0, 1.0, 1.0, 2.5, 6.0, 2.0, 2.5, 2.5,
92
0.2, 0.2, 1.0, 1.0, 2.3, 3.0
93
}; /*fuel depth, ft. */
95
{ {0, 3500, 3000, 1500, 2000, 2000, 1750, 1750, 2000, 2500, 2000, 1500,
97
{0, 0, 109, 0, 109, 109, 109, 109, 109, 109, 109, 109, 109, 109},
98
{0, 0, 30, 0, 30, 0, 30, 30, 30, 30, 30, 30, 30, 30},
99
{0, 0, 1500, 0, 1500, 1500, 0, 1500, 0, 0, 1500, 0, 0, 0}
102
/*fuel particale surface-area-to-volume ratio, 1/ft. */
103
float MX[] = { 0, 0.12, 0.15, 0.25, 0.20, 0.20, 0.25, 0.40,
104
0.30, 0.25, 0.25, 0.15, 0.20, 0.25
105
}; /*moisture content of extinction */
107
CELL *map_elev; /*full array for elevation map layer (for spotting) */
109
struct Cell_head window;
112
int main(int argc, char *argv[])
115
/***input of Rothermel equation (1972)***/
116
float h = 8000.0, /*heat of combustion, BTU/lb. */
117
rhop = 32, /*ovendry fuel density, lb./ft.^3 */
118
ST = 0.0555; /*fuel total mineral content, lb. minerals/lb. ovendry */
122
/***derived parameters of Rothermel equation (1972)***/
123
float R, /*rate of spread, ft./min.
124
R = IR*xi*(1+phiw+phis)/(rhob*epsilon*Qig) */
125
IR, /*reaction intensity, BTU/ft.^2/min.
126
IR = gamma*wn*h*etaM*etas */
127
gamma, /*optimum reation velosity, 1/min.
128
gamma = gammamax*(beta/betaop)^A*
129
exp(A(1-beta/betaop)) */
130
gammamax, /*maximum reation velosity, 1/min.
131
gammamax = sigma^1.5/(495+0.0594*sigma^1.5) */
132
betaop, /*optimum packing ratio
133
betaop = 3.348/(sigma^0.8189) */
134
A, /*A = 1/(4.77*sigma^0.1-7.27) */
135
etas = 0.174 / pow(0.01, 0.19), /*mineral damping coefficient */
136
xi, /*propagating flux ratio,
137
xi = exp((0.792+0.681*sigma^0.5)(beta+0.1))/
138
(192+0.2595*sigma) */
139
phiw, /*wind coefficient, phiw = C*U^B/(beta/betaop)^E */
140
C, /*C = 7.47/exp(0.133*sigma^0.55) */
141
B, /*B = 0.02526*sigma^.054 */
142
E, /*E = 0.715/exp(0.000359*sigma) */
143
phis, /*slope coefficient,phis = 5.275/beta^0.3*(tan(theta)^2) */
144
rhob, /*ovendry bulk density, lb./ft.^3, rohb = wo/delta */
145
epsilon[4][14], /*effective heating number, epsilon = 1/exp(138/sigma) */
146
Qig[14], /*heat of preignition, BTU/lb. Qig = 250+1116*Mf */
147
beta; /*packing ratio, beta = rhob/rhop */
149
/***intermediate variables***/
150
float R0, /*base ROS (w/out wind/slope) */
151
Rdir, sin_fac, cos_fac, Ffactor_all[4][14], /*in all fuel subclasses by sigma/WO */
152
Ffactor_in_dead[3][14], /*in dead fuel subclasses by sigma/WO */
153
Ffactor_dead[14], /*dead fuel weight by sigma/WO */
154
Ffactor_live[14], /*live fuel weight by sigma/WO */
155
Gfactor_in_dead[3][14], /*in dead fuel by the 6 classes */
156
G1, G2, G3, G4, G5, wo_dead[14], /*dead fuel total load */
157
wn_dead, /*net dead fuel total load */
158
wn_live, /*net live fuel (total) load */
159
class_sum, moisture[4], /*moistures of 1-h,10-h,100-h,live fuels */
160
Mf_dead, /*total moisture of dead fuels */
161
etaM_dead, /*dead fuel misture damping coefficent */
162
etaM_live, /*live fuel misture damping coefficent */
163
xmext, /*live fuel moisture of extinction */
164
phi_ws, /*wind and slope conbined coefficient */
165
wmfd, fdmois, fined, finel;
167
/*other local variables */
172
mois_1h_fd = 0, mois_10h_fd = 0, mois_100h_fd = 0, mois_live_fd = 0,
173
vel_fd = 0, dir_fd = 0,
174
elev_fd = 0, slope_fd = 0, aspect_fd = 0,
175
base_fd = 0, max_fd = 0, maxdir_fd = 0, spotdist_fd = 0;
177
char name_base[60], name_max[60], name_maxdir[60], name_spotdist[60];
179
CELL *fuel, /*cell buffer for fuel model map layer */
180
*mois_1h, /*cell buffer for 1-hour fuel moisture map layer */
181
*mois_10h, /*cell buffer for 10-hour fuel moisture map layer */
182
*mois_100h, /*cell buffer for 100-hour fuel moisture map layer */
183
*mois_live, /*cell buffer for live fuel moisture map layer */
184
*vel, /*cell buffer for wind velocity map layer */
185
*dir, /*cell buffer for wind direction map layer */
186
*elev = NULL, /*cell buffer for elevation map layer (for spotting) */
187
*slope, /*cell buffer for slope map layer */
188
*aspect, /*cell buffer for aspect map layer */
189
*base, /*cell buffer for base ROS map layer */
190
*max, /*cell buffer for max ROS map layer */
191
*maxdir, /*cell buffer for max ROS direction map layer */
192
*spotdist = NULL; /*cell buffer for max spotting distance map layer */
194
extern struct Cell_head window;
198
struct Option *model,
199
*mois_1h, *mois_10h, *mois_100h, *mois_live,
200
*vel, *dir, *elev, *slope, *aspect, *output;
203
/* please, remove before GRASS 7 released */
204
struct Flag *flag1, *flag2;
205
struct GModule *module;
207
/* initialize access to database and create temporary files */
210
/* Set description */
211
module = G_define_module();
212
module->keywords = _("raster, fire");
213
module->description =
214
_("Generates three, or four raster map layers showing 1) the base "
215
"(perpendicular) rate of spread (ROS), 2) the maximum (forward) ROS, "
216
"3) the direction of the maximum ROS, and optionally 4) the "
217
"maximum potential spotting distance.");
219
parm.model = G_define_option();
220
parm.model->key = "model";
221
parm.model->type = TYPE_STRING;
222
parm.model->required = YES;
223
parm.model->gisprompt = "old,cell,raster";
224
parm.model->description = _("Name of raster map containing fuel MODELs");
226
parm.mois_1h = G_define_option();
227
parm.mois_1h->key = "moisture_1h";
228
parm.mois_1h->type = TYPE_STRING;
229
parm.mois_1h->gisprompt = "old,cell,raster";
230
parm.mois_1h->description =
231
_("Name of raster map containing the 1-HOUR fuel MOISTURE (%)");
233
parm.mois_10h = G_define_option();
234
parm.mois_10h->key = "moisture_10h";
235
parm.mois_10h->type = TYPE_STRING;
236
parm.mois_10h->gisprompt = "old,cell,raster";
237
parm.mois_10h->description =
238
_("Name of raster map containing the 10-HOUR fuel MOISTURE (%)");
240
parm.mois_100h = G_define_option();
241
parm.mois_100h->key = "moisture_100h";
242
parm.mois_100h->type = TYPE_STRING;
243
parm.mois_100h->gisprompt = "old,cell,raster";
244
parm.mois_100h->description =
245
_("Name of raster map containing the 100-HOUR fuel MOISTURE (%)");
247
parm.mois_live = G_define_option();
248
parm.mois_live->key = "moisture_live";
249
parm.mois_live->type = TYPE_STRING;
250
parm.mois_live->required = YES;
251
parm.mois_live->gisprompt = "old,cell,raster";
252
parm.mois_live->description =
253
_("Name of raster map containing LIVE fuel MOISTURE (%)");
255
parm.vel = G_define_option();
256
parm.vel->key = "velocity";
257
parm.vel->type = TYPE_STRING;
258
parm.vel->gisprompt = "old,cell,raster";
259
parm.vel->description =
260
_("Name of raster map containing midflame wind VELOCITYs (ft/min)");
262
parm.dir = G_define_option();
263
parm.dir->key = "direction";
264
parm.dir->type = TYPE_STRING;
265
parm.dir->gisprompt = "old,cell,raster";
266
parm.dir->description =
267
_("Name of raster map containing wind DIRECTIONs (degree)");
269
parm.slope = G_define_option();
270
parm.slope->key = "slope";
271
parm.slope->type = TYPE_STRING;
272
parm.slope->gisprompt = "old,cell,raster";
273
parm.slope->description =
274
_("Name of raster map containing SLOPE (degree)");
276
parm.aspect = G_define_option();
277
parm.aspect->key = "aspect";
278
parm.aspect->type = TYPE_STRING;
279
parm.aspect->gisprompt = "old,cell,raster";
280
parm.aspect->description =
281
_("Name of raster map containing ASPECT (degree, anti-clockwise from E)");
283
parm.elev = G_define_option();
284
parm.elev->key = "elevation";
285
parm.elev->type = TYPE_STRING;
286
parm.elev->gisprompt = "old,cell,raster";
287
parm.elev->description =
288
_("Name of raster map containing ELEVATION (m) (required w/ -s)");
290
parm.output = G_define_option();
291
parm.output->key = "output";
292
parm.output->type = TYPE_STRING;
293
parm.output->required = YES;
294
parm.output->gisprompt = "new,cell,raster";
295
parm.output->description =
296
_("Name of raster map to contain results (several new layers)");
298
/* please, remove before GRASS 7 released */
299
flag1 = G_define_flag();
301
flag1->description = _("Run verbosely");
303
flag2 = G_define_flag();
305
flag2->description = _("Also produce maximum SPOTTING distance");
307
/* Parse command line */
308
if (G_parser(argc, argv))
311
/* please, remove before GRASS 7 released */
313
putenv("GRASS_VERBOSE=3");
314
G_warning(_("The '-v' flag is superseded and will be removed "
315
"in future. Please use '--verbose' instead."));
318
spotting = flag2->answer;
320
/* Check if input layers exists in data base */
321
if (G_find_cell2(parm.model->answer, "") == NULL)
322
G_fatal_error(_("Raster map <%s> not found"), parm.model->answer);
325
(parm.mois_1h->answer || parm.mois_10h->answer ||
326
parm.mois_100h->answer)) {
328
("no dead fuel moisture is given. At least one of the 1-h, 10-h, 100-h moisture layers is required.");
333
if (parm.mois_1h->answer) {
334
if (G_find_cell2(parm.mois_1h->answer, "") == NULL)
335
G_fatal_error(_("Raster map <%s> not found"),
336
parm.mois_1h->answer);
338
if (parm.mois_10h->answer) {
339
if (G_find_cell2(parm.mois_10h->answer, "") == NULL)
340
G_fatal_error(_("Raster map <%s> not found"),
341
parm.mois_10h->answer);
343
if (parm.mois_100h->answer) {
344
if (G_find_cell2(parm.mois_100h->answer, "") == NULL)
345
G_fatal_error(_("Raster map <%s> not found"),
346
parm.mois_100h->answer);
349
if (G_find_cell2(parm.mois_live->answer, "") == NULL)
350
G_fatal_error(_("Raster map <%s> not found"), parm.mois_live->answer);
352
if (parm.vel->answer && !(parm.dir->answer)) {
354
("a wind direction layer should be given if the wind velocity layer--%s-- has been given\n",
359
if (!(parm.vel->answer) && parm.dir->answer) {
361
("a wind velocity layer should be given if the wind direction layer--%s-- has been given\n",
366
if (parm.vel->answer) {
367
if (G_find_cell2(parm.vel->answer, "") == NULL)
368
G_fatal_error(_("Raster map <%s> not found"), parm.vel->answer);
370
if (parm.dir->answer) {
371
if (G_find_cell2(parm.dir->answer, "") == NULL)
372
G_fatal_error(_("Raster map <%s> not found"), parm.dir->answer);
375
if (parm.slope->answer && !(parm.aspect->answer)) {
377
("an aspect layer should be given if the slope layer--%s-- has been given\n",
382
if (!(parm.slope->answer) && parm.aspect->answer) {
384
("a slope layer should be given if the aspect layer--%s-- has been given\n",
385
parm.aspect->answer);
389
if (parm.slope->answer) {
390
if (G_find_cell2(parm.slope->answer, "") == NULL)
391
G_fatal_error(_("Raster map <%s> not found"), parm.slope->answer);
393
if (parm.aspect->answer) {
394
if (G_find_cell2(parm.aspect->answer, "") == NULL)
395
G_fatal_error(_("Raster map <%s> not found"),
396
parm.aspect->answer);
400
if (!(parm.elev->answer)) {
402
("an elevation layer should be given if considering spotting\n");
407
if (G_find_cell2(parm.elev->answer, "") == NULL)
408
G_fatal_error(_("Raster map <%s> not found"),
413
/* Check if specified output layer name IS LEGAL */
414
if (G_legal_filename(parm.output->answer) < 0)
415
G_fatal_error("%s - illegal name", parm.output->answer);
417
/*assign names of the three output ROS layers */
418
sprintf(name_base, "%s.base", parm.output->answer);
419
sprintf(name_max, "%s.max", parm.output->answer);
420
sprintf(name_maxdir, "%s.maxdir", parm.output->answer);
422
/*check if the output layer names EXIST */
423
if (G_find_cell2(name_base, G_mapset()))
424
G_fatal_error(_("Raster map <%s> already exists in mapset <%s>, select another name"),
425
name_base, G_mapset());
427
if (G_find_cell2(name_max, G_mapset()))
428
G_fatal_error(_("Raster map <%s> already exists in mapset <%s>, select another name"),
429
name_max, G_mapset());
431
if (G_find_cell2(name_maxdir, G_mapset()))
432
G_fatal_error(_("Raster map <%s> already exists in mapset <%s>, select another name"),
433
name_maxdir, G_mapset());
435
/*assign a name to output SPOTTING distance layer */
437
sprintf(name_spotdist, "%s.spotdist", parm.output->answer);
438
if (G_find_cell2(name_spotdist, G_mapset()))
439
G_fatal_error(_("Raster map <%s> already exists in mapset <%s>, select another name"),
440
name_spotdist, G_mapset());
443
/* Get database window parameters */
444
if (G_get_window(&window) < 0)
445
G_fatal_error("can't read current window parameters");
447
/* find number of rows and columns in window */
448
nrows = G_window_rows();
449
ncols = G_window_cols();
451
fuel = G_allocate_cell_buf();
452
mois_1h = G_allocate_cell_buf();
453
mois_10h = G_allocate_cell_buf();
454
mois_100h = G_allocate_cell_buf();
455
mois_live = G_allocate_cell_buf();
456
vel = G_allocate_cell_buf();
457
dir = G_allocate_cell_buf();
458
slope = G_allocate_cell_buf();
459
aspect = G_allocate_cell_buf();
460
base = G_allocate_cell_buf();
461
max = G_allocate_cell_buf();
462
maxdir = G_allocate_cell_buf();
464
spotdist = G_allocate_cell_buf();
465
elev = G_allocate_cell_buf();
466
map_elev = (CELL *) G_calloc(nrows * ncols, sizeof(CELL));
469
/* Open input cell layers for reading */
472
G_open_cell_old(parm.model->answer,
473
G_find_cell2(parm.model->answer, ""));
475
G_fatal_error(_("Unable to open raster map <%s>"),
478
if (parm.mois_1h->answer) {
480
G_open_cell_old(parm.mois_1h->answer,
481
G_find_cell2(parm.mois_1h->answer, ""));
483
G_fatal_error(_("Unable to open raster map <%s>"),
484
parm.mois_1h->answer);
486
if (parm.mois_10h->answer) {
488
G_open_cell_old(parm.mois_10h->answer,
489
G_find_cell2(parm.mois_10h->answer, ""));
491
G_fatal_error(_("Unable to open raster map <%s>"),
492
parm.mois_10h->answer);
494
if (parm.mois_100h->answer) {
496
G_open_cell_old(parm.mois_100h->answer,
497
G_find_cell2(parm.mois_100h->answer, ""));
498
if (mois_100h_fd < 0)
499
G_fatal_error(_("Unable to open raster map <%s>"),
500
parm.mois_100h->answer);
504
G_open_cell_old(parm.mois_live->answer,
505
G_find_cell2(parm.mois_live->answer, ""));
506
if (mois_live_fd < 0)
507
G_fatal_error(_("Unable to open raster map <%s>"),
508
parm.mois_live->answer);
510
if (parm.vel->answer) {
512
G_open_cell_old(parm.vel->answer,
513
G_find_cell2(parm.vel->answer, ""));
515
G_fatal_error(_("Unable to open raster map <%s>"),
518
if (parm.dir->answer) {
520
G_open_cell_old(parm.dir->answer,
521
G_find_cell2(parm.dir->answer, ""));
523
G_fatal_error(_("Unable to open raster map <%s>"),
527
if (parm.slope->answer) {
529
G_open_cell_old(parm.slope->answer,
530
G_find_cell2(parm.slope->answer, ""));
532
G_fatal_error(_("Unable to open raster map <%s>"),
535
if (parm.aspect->answer) {
537
G_open_cell_old(parm.aspect->answer,
538
G_find_cell2(parm.aspect->answer, ""));
540
G_fatal_error(_("Unable to open raster map <%s>"),
541
parm.aspect->answer);
546
G_open_cell_old(parm.elev->answer,
547
G_find_cell2(parm.elev->answer, ""));
549
G_fatal_error(_("Unable to open raster map <%s>"),
553
base_fd = G_open_cell_new(name_base);
554
max_fd = G_open_cell_new(name_max);
555
maxdir_fd = G_open_cell_new(name_maxdir);
557
spotdist_fd = G_open_cell_new(name_spotdist);
559
/*compute weights, combined wo, and combined sigma */
560
/*wo[model] -- simple sum of WO[class][model] by all fuel subCLASS */
561
/*sigma[model] -- weighted sum of SIGMA[class][model] by all fuel subCLASS *epsilon[class][model] */
562
for (model = 1; model <= 13; model++) {
564
wo_dead[model] = 0.0;
566
for (class = 0; class <= 3; class++) {
567
class_sum = class_sum + WO[class][model] * SIGMA[class][model];
568
if (SIGMA[class][model] > 0.0) {
569
epsilon[class][model] = exp(-138.0 / SIGMA[class][model]);
572
epsilon[class][model] = 0.0;
575
for (class = 0; class <= 3; class++) {
576
Ffactor_all[class][model] =
577
WO[class][model] * SIGMA[class][model] / class_sum;
580
SIGMA[class][model] * Ffactor_all[class][model];
583
for (class = 0; class <= 2; class++) {
584
wo_dead[model] = wo_dead[model] + WO[class][model];
585
class_sum = class_sum + WO[class][model] * SIGMA[class][model];
587
for (class = 0; class <= 2; class++) {
588
Ffactor_in_dead[class][model] =
589
WO[class][model] * SIGMA[class][model] / class_sum;
592
/* compute G factor for each of the 6 subclasses */
598
for (class = 0; class <= 2; class++) {
599
if (SIGMA[class][model] >= 1200)
600
G1 = G1 + Ffactor_in_dead[class][model];
601
if (SIGMA[class][model] < 1200 && SIGMA[class][model] >= 192)
602
G2 = G2 + Ffactor_in_dead[class][model];
603
if (SIGMA[class][model] < 192 && SIGMA[class][model] >= 96)
604
G3 = G3 + Ffactor_in_dead[class][model];
605
if (SIGMA[class][model] < 96 && SIGMA[class][model] >= 48)
606
G4 = G4 + Ffactor_in_dead[class][model];
607
if (SIGMA[class][model] < 48 && SIGMA[class][model] >= 16)
608
G5 = G5 + Ffactor_in_dead[class][model];
610
for (class = 0; class <= 2; class++) {
611
if (SIGMA[class][model] >= 1200)
612
Gfactor_in_dead[class][model] = G1;
613
if (SIGMA[class][model] < 1200 && SIGMA[class][model] >= 192)
614
Gfactor_in_dead[class][model] = G2;
615
if (SIGMA[class][model] < 192 && SIGMA[class][model] >= 96)
616
Gfactor_in_dead[class][model] = G3;
617
if (SIGMA[class][model] < 96 && SIGMA[class][model] >= 48)
618
Gfactor_in_dead[class][model] = G4;
619
if (SIGMA[class][model] < 48 && SIGMA[class][model] >= 16)
620
Gfactor_in_dead[class][model] = G5;
621
if (SIGMA[class][model] < 16)
622
Gfactor_in_dead[class][model] = 0.0;
625
Ffactor_dead[model] =
626
class_sum / (class_sum + WO[3][model] * SIGMA[3][model]);
627
Ffactor_live[model] = 1 - Ffactor_dead[model];
630
/*if considering spotting, read elevation map into an array */
632
for (row = 0; row < nrows; row++) {
633
if (G_get_map_row(elev_fd, elev, row) < 0)
634
G_fatal_error("cannot get map row!");
635
for (col = 0; col < ncols; col++)
636
DATA(map_elev, row, col) = elev[col];
639
/*major computation: compute ROSs one cell a time */
640
G_message(_("Percent Completed ... "));
642
for (row = 0; row < nrows; row++) {
643
G_percent(row, nrows, 2);
644
if (G_get_map_row(fuel_fd, fuel, row) < 0)
645
G_fatal_error("cannot get map row: %d!", row);
646
if (parm.mois_1h->answer)
647
if (G_get_map_row(mois_1h_fd, mois_1h, row) < 0)
648
G_fatal_error("cannot get map row: %d!", row);
649
if (parm.mois_10h->answer)
650
if (G_get_map_row(mois_10h_fd, mois_10h, row) < 0)
651
G_fatal_error("cannot get map row: %d!", row);
652
if (parm.mois_100h->answer)
653
if (G_get_map_row(mois_100h_fd, mois_100h, row) < 0)
654
G_fatal_error("cannot get map row: %d!", row);
655
if (G_get_map_row(mois_live_fd, mois_live, row) < 0)
656
G_fatal_error("cannot get map row: %d!", row);
657
if (parm.vel->answer)
658
if (G_get_map_row(vel_fd, vel, row) < 0)
659
G_fatal_error("cannot get map row: %d!", row);
660
if (parm.dir->answer)
661
if (G_get_map_row(dir_fd, dir, row) < 0)
662
G_fatal_error("cannot get map row: %d!", row);
663
if (parm.slope->answer)
664
if (G_get_map_row(slope_fd, slope, row) < 0)
665
G_fatal_error("cannot get map row: %d!", row);
666
if (parm.aspect->answer)
667
if (G_get_map_row(aspect_fd, aspect, row) < 0)
668
G_fatal_error("cannot get map row: %d!", row);
670
/*initialize cell buffers for output map layers */
671
for (col = 0; col < ncols; col++) {
672
base[col] = max[col] = maxdir[col] = 0;
677
for (col = 0; col < ncols; col++) {
678
/*check if a fuel is within the 13 models,
679
*if not, no processing; useful when no data presents*/
680
if (fuel[col] < 1 || fuel[col] > 13)
682
if (parm.mois_1h->answer)
683
moisture[0] = 0.01 * mois_1h[col];
684
if (parm.mois_10h->answer)
685
moisture[1] = 0.01 * mois_10h[col];
686
if (parm.mois_100h->answer)
687
moisture[2] = 0.01 * mois_100h[col];
688
moisture[3] = 0.01 * mois_live[col];
689
if (parm.aspect->answer)
690
aspect[col] = (630 - aspect[col]) % 360;
692
/* assign some dead fuel moisture if not completely
693
* given based on emperical relationship*/
694
if (!(parm.mois_10h->answer || parm.mois_100h->answer)) {
695
moisture[1] = moisture[0] + 0.01;
696
moisture[2] = moisture[0] + 0.02;
698
if (!(parm.mois_1h->answer || parm.mois_100h->answer)) {
699
moisture[0] = moisture[1] - 0.01;
700
moisture[2] = moisture[1] + 0.01;
702
if (!(parm.mois_1h->answer || parm.mois_10h->answer)) {
703
moisture[0] = moisture[2] - 0.02;
704
moisture[1] = moisture[2] - 0.01;
706
if (!(parm.mois_1h->answer) && parm.mois_10h->answer &&
707
parm.mois_100h->answer)
708
moisture[0] = moisture[1] - 0.01;
709
if (!(parm.mois_10h->answer) && parm.mois_1h->answer &&
710
parm.mois_100h->answer)
711
moisture[1] = moisture[0] + 0.01;
712
if (!(parm.mois_100h->answer) && parm.mois_1h->answer &&
713
parm.mois_10h->answer)
714
moisture[2] = moisture[1] + 0.01;
716
/*compute xmext, moisture of extinction of live fuels */
719
if (SIGMA[3][fuel[col]] > 0.0) {
720
for (class = 0; class <= 2; class++) {
721
if (SIGMA[class][fuel[col]] == 0.0)
725
WO[class][fuel[col]] * exp(-138.0 /
726
SIGMA[class][fuel[col]]);
729
WO[class][fuel[col]] * exp(-138.0 /
730
SIGMA[class][fuel[col]]) *
733
fdmois = wmfd / fined;
734
finel = WO[3][fuel[col]] * exp(-500.0 / SIGMA[3][fuel[col]]);
736
2.9 * (fined / finel) * (1 - fdmois / MX[fuel[col]]) -
740
xmext = MX[fuel[col]];
741
if (xmext < MX[fuel[col]])
742
xmext = MX[fuel[col]];
744
/*compute other intermediate values */
748
for (class = 0; class <= 2; class++) {
751
moisture[class] * Ffactor_in_dead[class][fuel[col]];
754
WO[class][fuel[col]] * Gfactor_in_dead[class][fuel[col]] *
756
Qig[class] = 250 + 1116 * moisture[class];
759
Ffactor_all[class][fuel[col]] *
760
epsilon[class][fuel[col]] * Qig[class];
763
1.0 - 2.59 * (Mf_dead / MX[fuel[col]]) +
764
5.11 * (Mf_dead / MX[fuel[col]]) * (Mf_dead / MX[fuel[col]]) -
765
3.52 * (Mf_dead / MX[fuel[col]]) * (Mf_dead / MX[fuel[col]]) *
766
(Mf_dead / MX[fuel[col]]);
767
if (Mf_dead >= MX[fuel[col]])
770
1.0 - 2.59 * (moisture[3] / xmext) +
771
5.11 * (moisture[3] / xmext) * (moisture[3] / xmext) -
772
3.52 * (moisture[3] / xmext) * (moisture[3] / xmext) *
773
(moisture[3] / xmext);
774
if (moisture[3] >= xmext)
776
wn_live = WO[3][fuel[col]] * (1 - ST);
777
Qig[3] = 250 + 1116 * moisture[3];
780
Ffactor_all[3][fuel[col]] * epsilon[3][fuel[col]] * Qig[3];
782
/*final computations */
783
rhob = (wo_dead[fuel[col]] + WO[3][fuel[col]]) / DELTA[fuel[col]];
785
betaop = 3.348 / pow(sigma[fuel[col]], 0.8189);
786
A = 133 / pow(sigma[fuel[col]], 0.7913);
788
pow(sigma[fuel[col]],
789
1.5) / (495 + 0.0594 * pow(sigma[fuel[col]], 1.5));
791
gammamax * pow(beta / betaop,
792
A) * exp(A * (1 - beta / betaop));
793
xi = exp((0.792 + 0.681 * pow(sigma[fuel[col]], 0.5)) * (beta +
795
(192 + 0.2595 * sigma[fuel[col]]);
796
IR = gamma * h * (wn_dead * etaM_dead +
797
wn_live * etaM_live) * etas;
799
R0 = IR * xi / (rhob * class_sum);
801
if (parm.vel->answer && parm.dir->answer) {
802
C = 7.47 * exp(-0.133 * pow(sigma[fuel[col]], 0.55));
803
B = 0.02526 * pow(sigma[fuel[col]], 0.54);
804
E = 0.715 * exp(-0.000359 * sigma[fuel[col]]);
805
phiw = C * pow((double)vel[col], B) / pow(beta / betaop, E);
810
if (parm.slope->answer && parm.aspect->answer) {
813
-0.3) * tan(0.01745 * slope[col]) *
814
tan(0.01745 * slope[col]);
819
/*compute the maximum ROS R and its direction,
820
*vector adding for the effect of wind and slope*/
821
if (parm.dir->answer && parm.aspect->answer) {
823
phiw * sin(0.01745 * dir[col]) +
824
phis * sin(0.01745 * aspect[col]);
826
phiw * cos(0.01745 * dir[col]) +
827
phis * cos(0.01745 * aspect[col]);
828
phi_ws = sqrt(sin_fac * sin_fac + cos_fac * cos_fac);
829
Rdir = atan2(sin_fac, cos_fac) / 0.01745;
831
else if (parm.dir->answer && !(parm.aspect->answer)) {
835
else if (!(parm.dir->answer) && parm.aspect->answer) {
837
Rdir = (float)aspect[col];
843
R = R0 * (1 + phi_ws);
846
/*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); */
848
/*maximum spotting distance */
851
spot_dist(fuel[col], R, vel[col], Rdir, row, col);
853
/*to avoid 0 R, convert ft./min to cm/min */
856
/*4debugging R0 = R0/30.5/1.1; R = R/30.5/1.1; */
860
maxdir[col] = (int)Rdir;
861
/*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); */
863
G_put_raster_row(base_fd, base, CELL_TYPE);
864
G_put_raster_row(max_fd, max, CELL_TYPE);
865
G_put_raster_row(maxdir_fd, maxdir, CELL_TYPE);
867
G_put_raster_row(spotdist_fd, spotdist, CELL_TYPE);
869
G_percent(row, nrows, 2);
871
G_close_cell(fuel_fd);
872
if (parm.mois_1h->answer)
873
G_close_cell(mois_1h_fd);
874
if (parm.mois_10h->answer)
875
G_close_cell(mois_10h_fd);
876
if (parm.mois_100h->answer)
877
G_close_cell(mois_100h_fd);
878
G_close_cell(mois_live_fd);
879
if (parm.vel->answer)
880
G_close_cell(vel_fd);
881
if (parm.dir->answer)
882
G_close_cell(dir_fd);
883
if (parm.slope->answer)
884
G_close_cell(slope_fd);
885
if (parm.aspect->answer)
886
G_close_cell(aspect_fd);
887
G_close_cell(base_fd);
888
G_close_cell(max_fd);
889
G_close_cell(maxdir_fd);
891
G_close_cell(spotdist_fd);
896
for (model = 1; model <= 13; model++) {
898
printf("\n Grass and Grass-dominated\n");
900
printf(" Chaparral and Shrubfields\n");
902
printf(" Timber Litter\n");
903
else if (model == 11)
904
printf(" Logging Slash\n");
905
printf("Model %2d ", model);
906
for (class = 0; class <= 3; class++)
907
printf("%4.0f/%.3f ", SIGMA[class][model], WO[class][model]);
908
printf(" %.1f %.2f\n", DELTA[model], MX[model]);
911
printf("\nWeight in All Fuel Subclasses:\n");
912
for (model = 1; model <= 13; model++) {
913
printf("model %2d ", model);
914
for (class = 0; class <= 3; class++)
915
printf("%.2f ", Ffactor_all[class][model]);
916
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);
918
printf("\nWeight in Dead Fuel Subclasses, Dead Weight/Live Weight:\n");for (model = 1; model <= 13; model++) {
919
printf("model %2d ", model);
920
for (class = 0; class <= 2; class++)
921
printf("%.2f ", Ffactor_in_dead[class][model]);
922
printf("%.2f/%.2f model %2d\n", Ffactor_dead[model], Ffactor_live[model], model);
added
removed
raster/wildfire/r.ros/main.c
