| 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 |
|
|---|
| 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 | }
|
|---|
