diff --git a/src/fwi/C6calc.ts b/src/fwi/C6calc.ts new file mode 100644 index 0000000..c867f28 --- /dev/null +++ b/src/fwi/C6calc.ts @@ -0,0 +1,103 @@ +/** + * C-6 Conifer Plantation Fire Spread Calculator + * + * Calculate C6 (Conifer plantation) Fire Spread. C6 is a special + * case, and thus has its own function. To calculate C6 fire spread, this + * function also calculates and can return ROS, CFB, RSC, or RSI by specifying + * in the option parameter. + * All variables names are laid out in the same manner as Forestry Canada Fire + * Danger Group (FCFDG) (1992). Development and Structure of the Canadian + * Forest Fire Behavior Prediction System." Technical Report ST-X-3, Forestry + * Canada, Ottawa, Ontario. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param ISI Initial Spread Index + * @param BUI Buildup Index + * @param FMC Foliar Moisture Content + * @param SFC Surface Fuel Consumption + * @param CBH Crown Base Height + * @param ROS Rate of Spread + * @param CFB Crown Fraction Burned + * @param RSC Crown Fire Spread Rate (m/min) + * @param option Which variable to calculate ("ROS", "RSC", "RSI", or "CFB" (the default).) + * + * @return ROS, CFB, RSC or RSI depending on which option was selected + */ + +import { criticalSurfaceIntensity, surfaceFireRateOfSpread, crownFractionBurned } from "./CFBcalc" +import { buildupEffect } from "./buildup_effect"; + +function intermediateSurfaceRateOfSpreadC6(ISI: number): number { + // Eq. 62 (FCFDG 1992) Intermediate surface fire spread rate + return 30 * Math.pow(1 - Math.exp(-0.08 * ISI), 3.0); +} + +function surfaceRateOfSpreadC6(RSI: number, BUI: number): number { + // Eq. 63 (FCFDG 1992) Surface fire spread rate (m/min) + return RSI * buildupEffect("C6", BUI); +} + +function crownRateOfSpreadC6(ISI: number, FMC: number): number { + // Average foliar moisture effect + const FMEavg = 0.778; + // Eq. 59 (FCFDG 1992) Crown flame temperature (degrees K) + const tt = 1500 - 2.75 * FMC; + // Eq. 60 (FCFDG 1992) Head of ignition (kJ/kg) + const H = 460 + 25.9 * FMC; + // Eq. 61 (FCFDG 1992) Average foliar moisture effect + const FME = Math.pow(1.5 - 0.00275 * FMC, 4.0) / (460 + 25.9 * FMC) * 1000; + // Eq. 64 (FCFDG 1992) Crown fire spread rate (m/min) + return 60 * (1 - Math.exp(-0.0497 * ISI)) * FME / FMEavg; +} + +function crownFractionBurnedC6(RSC: number, RSS: number, RSO: number): number { + return (RSC > RSS) && (RSS > RSO) ? crownFractionBurned(RSS, RSO) : 0; +} + +function rateOfSpreadC6(RSC: number, RSS: number, CFB: number): number { + // Eq. 65 (FCFDG 1992) Calculate Rate of spread (m/min) + return RSC > RSS ? RSS + CFB * (RSC - RSS) : RSS; +} + +export function C6calc( + FUELTYPE: string, + ISI: number, + BUI: number, + FMC: number, + SFC: number, + CBH: number, + option: string = "CFB" +): number { + // This should be validated for the FUELTYPE + // if (FUELTYPE !== "C6") { + // throw new Error("Invalid FUELTYPE for C6 calculation"); + // } + + const RSI = intermediateSurfaceRateOfSpreadC6(ISI); + if (option === "RSI") { + console.warn("Deprecated: intermediateSurfaceRateOfSpreadC6"); + return RSI; + } + + const RSC = crownRateOfSpreadC6(ISI, FMC); + if (option === "RSC") { + console.warn("Deprecated: crownRateOfSpreadC6"); + return RSC; + } + + const RSS = surfaceRateOfSpreadC6(RSI, BUI); + const CSI = criticalSurfaceIntensity(FMC, CBH); + const RSO = surfaceFireRateOfSpread(CSI, SFC); + const CFB = crownFractionBurnedC6(RSC, RSS, RSO); + if (option === "CFB") { + console.warn("Deprecated: crownFractionBurnedC6"); + return CFB; + } + + console.warn("Deprecated: rateOfSpreadC6"); + const ROS = rateOfSpreadC6(RSC, RSS, CFB); + return ROS; +} + +// Exporting functions for external usage +export { intermediateSurfaceRateOfSpreadC6, surfaceRateOfSpreadC6, crownRateOfSpreadC6, crownFractionBurnedC6, rateOfSpreadC6 }; diff --git a/src/fwi/CFBcalc.ts b/src/fwi/CFBcalc.ts new file mode 100644 index 0000000..7801821 --- /dev/null +++ b/src/fwi/CFBcalc.ts @@ -0,0 +1,95 @@ +/** + * Crown Fraction Burned Calculator + * + * Calculate Calculate Crown Fraction Burned. To calculate CFB, we + * also need to calculate Critical surface intensity (CSI), and Surface fire + * rate of spread (RSO). The value of each of these equations can be returned + * to the calling function without unnecessary additional calculations. + * + * All variables names are laid out in the same manner as Forestry Canada Fire + * Danger Group (FCFDG) (1992). Development and Structure of the Canadian + * Forest Fire Behavior Prediction System." Technical Report ST-X-3, Forestry + * Canada, Ottawa, Ontario. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param FMC Foliar Moisture Content + * @param SFC Surface Fuel Consumption + * @param CBH Crown Base Height + * @param ROS Rate of Spread + * @param option Which variable to calculate ("ROS", "RSC", "RSI", or "CFB" (the default).) + * @returns CFB, CSI, RSO depending on which option was selected. + */ + +/** + * Calculate Critical Surface Intensity (CSI). + * + * @param FMC Foliar Moisture Content + * @param CBH Crown Base Height + * @returns Critical Surface Intensity (CSI) + */ +function criticalSurfaceIntensity(FMC: number, CBH: number): number { + return 0.001 * (Math.pow(CBH, 1.5)) * Math.pow((460 + 25.9 * FMC), 1.5); +} + +/** + * Calculate Surface Fire Rate of Spread (RSO). + * + * @param CSI Critical Surface Intensity + * @param SFC Surface Fuel Consumption + * @returns Surface Fire Rate of Spread (RSO) + */ +function surfaceFireRateOfSpread(CSI: number, SFC: number): number { + return CSI / (300 * SFC); +} + +/** + * Calculate Crown Fraction Burned (CFB). + * + * @param ROS Rate of Spread + * @param RSO Surface Fire Rate of Spread + * @returns Crown Fraction Burned (CFB) + */ +function crownFractionBurned(ROS: number, RSO: number): number { + return ROS > RSO ? 1 - Math.exp(-0.23 * (ROS - RSO)) : 0; +} + +/** + * Main function to calculate various fire behavior metrics. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param FMC Foliar Moisture Content + * @param SFC Surface Fuel Consumption + * @param ROS Rate of Spread + * @param CBH Crown Base Height + * @param option Which variable to calculate ("ROS", "RSC", "RSI", or "CFB" (the default).) + * @returns CFB, CSI, RSO depending on which option was selected. + */ +function CFBcalc( + FUELTYPE: string, + FMC: number, + SFC: number, + ROS: number, + CBH: number, + option: string = "CFB" +): number { + const CSI = criticalSurfaceIntensity(FMC, CBH); + + if (option === "CSI") { + console.warn("Deprecated: criticalSurfaceIntensity"); + return CSI; + } + + const RSO = surfaceFireRateOfSpread(CSI, SFC); + + if (option === "RSO") { + console.warn("Deprecated: surfaceFireRateOfSpread"); + return RSO; + } + + const CFB = crownFractionBurned(ROS, RSO); + console.warn("Deprecated: crownFractionBurned"); + return CFB; +} + +// Export the function for external usage +export { CFBcalc, criticalSurfaceIntensity, surfaceFireRateOfSpread, crownFractionBurned }; diff --git a/src/fwi/Slopecalc.ts b/src/fwi/Slopecalc.ts new file mode 100644 index 0000000..fe31d4b --- /dev/null +++ b/src/fwi/Slopecalc.ts @@ -0,0 +1,329 @@ +/** + * Slope Adjusted wind speed or slope direction of spread calculation + * + * Calculate the net effective windspeed (WSV), the net effective + * wind direction (RAZ) or the wind azimuth (WAZ). + * + * All variables names are laid out in the same manner as FCFDG (1992) and + * Wotton (2009). + * + * Forestry Canada Fire Danger Group (FCFDG) (1992). "Development and + * Structure of the Canadian Forest Fire Behavior Prediction System." + * Technical Report ST-X-3, Forestry Canada, Ottawa, Ontario. + * + * Wotton, B.M., Alexander, M.E., Taylor, S.W. 2009. Updates and revisions to + * the 1992 Canadian forest fire behavior prediction system. Nat. Resour. + * Can., Can. For. Serv., Great Lakes For. Cent., Sault Ste. Marie, Ontario, + * Canada. Information Report GLC-X-10, 45p. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param FFMC Fine Fuel Moisture Code + * @param BUI The Buildup Index value + * @param WS Windspeed (km/h) + * @param WAZ Wind Azimuth + * @param GS Ground Slope (%) + * @param SAZ Slope Azimuth + * @param FMC Foliar Moisture Content + * @param SFC Surface Fuel Consumption (kg/m^2) + * @param PC Percent Conifer (%) + * @param PDF Percent Dead Balsam Fir (%) + * @param CC Constant + * @param CBH Crown Base Height (m) + * @param ISI Initial Spread Index + * + * @returns RAZ and WSV + * - Rate of spread azimuth (degrees) and Wind Slope speed (km/hr) + */ + +import { initialSpreadIndex } from "./initial_spread_index"; +import { rateOfSpread } from "./rate_of_spread"; + +interface SlopeAdjustmentResult { + WSV: number | null; + RAZ: number | null; +} + +interface ROSResult { + ROS: number[]; + CFB: number[]; + CSI: number[]; + RSO: number[]; +} + +/** + * Helper function to check if an element is in an array + */ +function inArray(element: T, array: T[]): boolean { + return array.includes(element); +} + +function slopeAdjustment( + FUELTYPE: string[], + FFMC: number[], + BUI: number[], + WS: number[], + WAZ: number[], + GS: number[], + SAZ: number[], + FMC: number[], + SFC: number[], + PC: number[], + PDF: number[], + CC: number[], + CBH: number[], + ISI: number[] +): SlopeAdjustmentResult[] { + const NoBUI = FFMC.map(() => -1); + + const SF = GS.map(gs => gs >= 70 ? 10 : Math.exp(3.533 * Math.pow(gs / 100, 1.2))); + const ISZ = FFMC.map(ffmc => initialSpreadIndex(ffmc, 0)); + const RSZ = rateOfSpread({ FUELTYPE, ISI: ISZ, BUI: NoBUI, FMC, SFC, PC, PDF, CC, CBH }); + const RSF = RSZ.map((rsz, i) => rsz * SF[i]); + + type FuelType = "C1" | "C2" | "C3" | "C4" | "C5" | "C6" | "C7" | + "D1" | "M1" | "M2" | "M3" | "M4" | "S1" | "S2" | "S3" | "O1A" | "O1B"; + + const d: FuelType[] = [ + "C1", "C2", "C3", "C4", "C5", "C6", "C7", + "D1", "M1", "M2", "M3", "M4", "S1", "S2", "S3", "O1A", "O1B" + ]; + + const a: Record = { + "C1": 90, "C2": 110, "C3": 110, "C4": 110, "C5": 30, "C6": 30, "C7": 45, + "D1": 30, "M1": 0, "M2": 0, "M3": 120, "M4": 100, "S1": 75, "S2": 40, "S3": 55, + "O1A": 190, "O1B": 250 + }; + + const b: Record = { + "C1": 0.0649, "C2": 0.0282, "C3": 0.0444, "C4": 0.0293, "C5": 0.0697, "C6": 0.0800, + "C7": 0.0305, "D1": 0.0232, "M1": 0, "M2": 0, "M3": 0.0572, "M4": 0.0404, + "S1": 0.0297, "S2": 0.0438, "S3": 0.0829, "O1A": 0.0310, "O1B": 0.0350 + }; + + const c0: Record = { + "C1": 4.5, "C2": 1.5, "C3": 3.0, "C4": 1.5, "C5": 4.0, "C6": 3.0, "C7": 2.0, + "D1": 1.6, "M1": 0, "M2": 0, "M3": 1.4, "M4": 1.48, "S1": 1.3, "S2": 1.7, + "S3": 3.2, "O1A": 1.4, "O1B": 1.7 + }; + + const isBasic = FUELTYPE.map(fuel => ["C1", "C2", "C3", "C4", "C5", "C6", "C7", "D1", "S1", "S2", "S3"].includes(fuel)); + + let ISF = new Array(FFMC.length).fill(-99); + + ISF = ISF.map((isf, i) => { + if (isBasic[i]) { + const rsf = RSF[i]; + const fuel = FUELTYPE[i] as FuelType; + const expValue = 1 - Math.pow(rsf / a[fuel], 1 / c0[fuel]); + return expValue >= 0.01 ? Math.log(expValue) / -b[fuel] : Math.log(0.01) / -b[fuel]; + } + return isf; + }); + + // When calculating the M1/M2 types, we are going to calculate for both C2 and D1 types, and combine + if (FUELTYPE.includes("M1") || FUELTYPE.includes("M2")) { + const RSZ_C2 = rateOfSpread({ + FUELTYPE: FUELTYPE.map(fuel => (fuel === "M1" || fuel === "M2" ? "C2" : fuel)), + ISI: ISZ, + BUI: NoBUI, + FMC, + SFC, + PC, + PDF, + CC, + CBH + }); + const RSF_C2 = RSZ_C2.map((rsz, i) => rsz * SF[i]); + + const RSZ_D1 = rateOfSpread({ + FUELTYPE: FUELTYPE.map(fuel => (fuel === "M1" || fuel === "M2" ? "D1" : fuel)), + ISI: ISZ, + BUI: NoBUI, + FMC, + SFC, + PC, + PDF, + CC, + CBH + }); + const RSF_D1 = RSZ_D1.map((rsz, i) => rsz * SF[i]); + + ISF = ISF.map((isf, i) => { + if (FUELTYPE[i] === "M1" || FUELTYPE[i] === "M2") { + const RSF0_C2 = 1 - Math.pow(RSF_C2[i] / a["C2"], 1 / c0["C2"]); + const ISF_C2 = RSF0_C2 >= 0.01 ? Math.log(RSF0_C2) / -b["C2"] : Math.log(0.01) / -b["C2"]; + + const RSF0_D1 = 1 - Math.pow(RSF_D1[i] / a["D1"], 1 / c0["D1"]); + const ISF_D1 = RSF0_D1 >= 0.01 ? Math.log(RSF0_D1) / -b["D1"] : Math.log(0.01) / -b["D1"]; + + return PC[i] / 100 * ISF_C2 + (1 - PC[i] / 100) * ISF_D1; + } + return isf; + }); + } + + // Set % Dead Balsam Fir to 100% for M3/M4 types + const PDF100 = ISI.map(() => 100); + + if (FUELTYPE.includes("M3")) { + const RSZ_M3 = rateOfSpread({ + FUELTYPE: FUELTYPE.map(fuel => (fuel === "M3" ? "M3" : fuel)), + ISI: ISZ, + BUI: NoBUI, + FMC, + SFC, + PC, + PDF: PDF100, + CC, + CBH + }); + const RSF_M3 = RSZ_M3.map((rsz, i) => rsz * SF[i]); + + const RSZ_D1 = rateOfSpread({ + FUELTYPE: FUELTYPE.map(fuel => (fuel === "M3" ? "D1" : fuel)), + ISI: ISZ, + BUI: NoBUI, + FMC, + SFC, + PC, + PDF: PDF100, + CC, + CBH + }); + const RSF_D1 = RSZ_D1.map((rsz, i) => rsz * SF[i]); + + ISF = ISF.map((isf, i) => { + if (FUELTYPE[i] === "M3") { + const RSF0_M3 = 1 - Math.pow(RSF_M3[i] / a["M3"], 1 / c0["M3"]); + const ISF_M3 = RSF0_M3 >= 0.01 ? Math.log(RSF0_M3) / -b["M3"] : Math.log(0.01) / -b["M3"]; + + const RSF0_D1 = 1 - Math.pow(RSF_D1[i] / a["D1"], 1 / c0["D1"]); + const ISF_D1 = RSF0_D1 >= 0.01 ? Math.log(RSF0_D1) / -b["D1"] : Math.log(0.01) / -b["D1"]; + + return PDF[i] / 100 * ISF_M3 + (1 - PDF[i] / 100) * ISF_D1; + } + return isf; + }); + } + + if (FUELTYPE.includes("M4")) { + const RSZ_M4 = rateOfSpread({ + FUELTYPE: FUELTYPE.map(fuel => (fuel === "M4" ? "M4" : fuel)), + ISI: ISZ, + BUI: NoBUI, + FMC, + SFC, + PC, + PDF: PDF100, + CC, + CBH + }); + const RSF_M4 = RSZ_M4.map((rsz, i) => rsz * SF[i]); + + const RSZ_D1 = rateOfSpread({ + FUELTYPE: FUELTYPE.map(fuel => (fuel === "M4" ? "D1" : fuel)), + ISI: ISZ, + BUI: NoBUI, + FMC, + SFC, + PC, + PDF: PDF100, + CC, + CBH + }); + const RSF_D1 = RSZ_D1.map((rsz, i) => rsz * SF[i]); + + ISF = ISF.map((isf, i) => { + if (FUELTYPE[i] === "M4") { + const RSF0_M4 = 1 - Math.pow(RSF_M4[i] / a["M4"], 1 / c0["M4"]); + const ISF_M4 = RSF0_M4 >= 0.01 ? Math.log(RSF0_M4) / -b["M4"] : Math.log(0.01) / -b["M4"]; + + const RSF0_D1 = 1 - Math.pow(RSF_D1[i] / a["D1"], 1 / c0["D1"]); + const ISF_D1 = RSF0_D1 >= 0.01 ? Math.log(RSF0_D1) / -b["D1"] : Math.log(0.01) / -b["D1"]; + + return PDF[i] / 100 * ISF_M4 + (1 - PDF[i] / 100) * ISF_D1; + } + return isf; + }); + } + + const CF = FUELTYPE.map((fuel, i) => { + if (fuel === "O1A" || fuel === "O1B") { + return CC[i] < 58.8 ? 0.005 * (Math.exp(0.061 * CC[i]) - 1) : 0.176 + 0.02 * (CC[i] - 58.8); + } + return -99; + }); + + ISF = ISF.map((isf, i) => { + if (FUELTYPE[i] === "O1A" || FUELTYPE[i] === "O1B") { + const rsf = RSF[i]; + const fuel = FUELTYPE[i]; + const expValue = 1 - Math.pow(rsf / (CF[i] * a[fuel]), 1 / c0[fuel]); + return expValue >= 0.01 ? Math.log(expValue) / -b[fuel] : Math.log(0.01) / -b[fuel]; + } + return isf; + }); + + if (FUELTYPE.includes("NF") || FUELTYPE.includes("WA")) { + return FFMC.map(() => ({ + WSV: null, + RAZ: null + })); + } + + const m = FFMC.map(ffmc => 101 - ffmc); + const fF = m.map(mVal => 91.9 * Math.exp(-0.1386 * mVal) * (1 + Math.pow(mVal, 5.31) / 4.93e7)); + const WSE = fF.map((ff, i) => 1 / 0.05039 * Math.log(ISF[i] / (0.208 * ff))); + + const adjustedWSE = WSE.map((wse, i) => { + if (wse > 40 && ISF[i] < (0.999 * 2.496 * fF[i])) { + return 28 - (1 / 0.0818 * Math.log(1 - ISF[i] / (2.496 * fF[i]))); + } else if (wse > 40 && ISF[i] >= (0.999 * 2.496 * fF[i])) { + return 112.45; + } + return wse; + }); + + const WSX = WS.map((ws, i) => ws * Math.sin(WAZ[i]) + adjustedWSE[i] * Math.sin(SAZ[i])); + const WSY = WS.map((ws, i) => ws * Math.cos(WAZ[i]) + adjustedWSE[i] * Math.cos(SAZ[i])); + const WSV = WSX.map((wsx, i) => Math.sqrt(wsx * wsx + WSY[i] * WSY[i])); + const RAZ = WSV.map((wsv, i) => Math.acos(WSY[i] / wsv)); + const adjustedRAZ = WSX.map((wsx, i) => wsx < 0 ? 2 * Math.PI - RAZ[i] : RAZ[i]); + + return FFMC.map((_, i) => ({ + WSV: WSV[i], + RAZ: adjustedRAZ[i] + })); +} + +function Slopecalc( + FUELTYPE: string[], + FFMC: number[], + BUI: number[], + WS: number[], + WAZ: number[], + GS: number[], + SAZ: number[], + FMC: number[], + SFC: number[], + PC: number[], + PDF: number[], + CC: number[], + CBH: number[], + ISI: number[], + output: string = "RAZ" +): SlopeAdjustmentResult[] | null { + + const validOutTypes = ["RAZ", "WAZ", "WSV"]; + if (!validOutTypes.includes(output)) { + console.warn("Deprecated: slopeAdjustment"); + throw new Error(`In 'Slopecalc()', '${output}' is an invalid 'output' type.`); + } + + const values = slopeAdjustment(FUELTYPE, FFMC, BUI, WS, WAZ, GS, SAZ, FMC, SFC, PC, PDF, CC, CBH, ISI); + values[0][output as keyof SlopeAdjustmentResult] || null; + return values; +} + +// Exporting functions for external usage +export { slopeAdjustment, Slopecalc }; \ No newline at end of file diff --git a/src/fwi/back_rate_of_spread.ts b/src/fwi/back_rate_of_spread.ts new file mode 100644 index 0000000..effe833 --- /dev/null +++ b/src/fwi/back_rate_of_spread.ts @@ -0,0 +1,77 @@ +/** + * Back Fire Rate of Spread Calculator + * + * Calculate the Back Fire Spread Rate. All variables names are + * laid out in the same manner as Forestry Canada Fire Danger Group (FCFDG) + * (1992). + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param FFMC Fine Fuel Moisture Code + * @param BUI Buildup Index + * @param WSV Wind Speed Vector + * @param FMC Foliar Moisture Content + * @param SFC Surface Fuel Consumption + * @param PC Percent Conifer + * @param PDF Percent Dead Balsam Fir + * @param CC Degree of Curing + * @param CBH Crown Base Height + * + * @return BROS: Back Fire Rate of Spread + */ + +import { rateOfSpread } from "./rate_of_spread"; + + +function backRateOfSpread( + FUELTYPE: string[], + FFMC: number[], + BUI: number[], + WSV: number[], + FMC: number[], + SFC: number[], + PC: number[], + PDF: number[], + CC: number[], + CBH: number[] +): number[] { + // Eq. 46 (FCFDG 1992) + // Calculate the FFMC function from the ISI equation + const FFMC_COEFFICIENT = 1; // Define the appropriate coefficient value + const m = FFMC_COEFFICIENT * (101 - FFMC[0]) / (59.5 + FFMC[0]); + + // Eq. 45 (FCFDG 1992) + const fF = 91.9 * Math.exp(-0.1386 * m) * (1.0 + Math.pow(m, 5.31) / 4.93e7); + + // Eq. 75 (FCFDG 1992) + // Calculate the Back fire wind function + const BfW = Math.exp(-0.05039 * WSV[0]); + + // Calculate the ISI associated with the back fire spread rate + // Eq. 76 (FCFDG 1992) + const BISI = [0.208 * BfW * fF]; + + // Eq. 77 (FCFDG 1992) + // Calculate final Back fire spread rate FUELTYPE, ISI: ISZ, BUI: NoBUI, FMC, SFC, PC, PDF, CC, CBH + const BROS = rateOfSpread({FUELTYPE, ISI: BISI, BUI, FMC, SFC, PC, PDF, CC, CBH}); + return BROS; +} + +function BROScalc( + FUELTYPE: string[], + FFMC: number[], + BUI: number[], + WSV: number[], + FMC: number[], + SFC: number[], + PC: number[], + PDF: number[], + CC: number[], + CBH: number[] +): number[] { + console.warn("Deprecated: backRateOfSpread"); + return backRateOfSpread(FUELTYPE, FFMC, BUI, WSV, FMC, SFC, PC, PDF, CC, CBH); +} + + +// Exporting functions for external usage +export { backRateOfSpread, BROScalc }; diff --git a/src/fwi/buildup_effect.ts b/src/fwi/buildup_effect.ts new file mode 100644 index 0000000..12cc15e --- /dev/null +++ b/src/fwi/buildup_effect.ts @@ -0,0 +1,52 @@ +/** + * Build Up Effect Calculator + * + * Computes the Buildup Effect on Fire Spread Rate. All variables + * names are laid out in the same manner as Forestry Canada Fire Danger Group + * (FCFDG)(1992). + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param BUI The Buildup Index value + * + * @return BE: Build up effect + */ + +function buildupEffect(FUELTYPE: string, BUI: number): number { + // Fuel Type String representations + const d = [ + "C1", "C2", "C3", "C4", "C5", "C6", "C7", + "D1", "M1", "M2", "M3", "M4", "S1", "S2", "S3", "O1A", "O1B" + ]; + + // The average BUI for the fuel type - as referenced by the "d" list above + const BUIo: { [key: string]: number } = { + "C1": 72, "C2": 64, "C3": 62, "C4": 66, "C5": 56, "C6": 62, "C7": 106, + "D1": 32, "M1": 50, "M2": 50, "M3": 50, "M4": 50, "S1": 38, "S2": 63, + "S3": 31, "O1A": 1, "O1B": 1 + }; + + // Proportion of maximum possible spread rate that is reached at a standard BUI + const Q: { [key: string]: number } = { + "C1": 0.9, "C2": 0.7, "C3": 0.75, "C4": 0.8, "C5": 0.8, "C6": 0.8, "C7": 0.85, + "D1": 0.9, "M1": 0.8, "M2": 0.8, "M3": 0.8, "M4": 0.8, "S1": 0.75, "S2": 0.75, + "S3": 0.75, "O1A": 1.0, "O1B": 1.0 + }; + + // Eq. 54 (FCFDG 1992) The Buildup Effect + let BE: number; + if (BUI > 0 && BUIo[FUELTYPE] > 0) { + BE = Math.exp(50 * Math.log(Q[FUELTYPE]) * (1 / BUI - 1 / BUIo[FUELTYPE])); + } else { + BE = 1; + } + + return BE; +} + +function BEcalc(FUELTYPE: string, BUI: number): number { + console.warn("Deprecated: buildupEffect"); + return buildupEffect(FUELTYPE, BUI); +} + +// Exporting functions for external usage +export { buildupEffect, BEcalc }; diff --git a/src/fwi/buildup_index.ts b/src/fwi/buildup_index.ts new file mode 100644 index 0000000..5aaa851 --- /dev/null +++ b/src/fwi/buildup_index.ts @@ -0,0 +1,35 @@ +/** + * Buildup Index Calculator + * + * @description Buildup Index Calculation. Based on the equations and FORTRAN program + * for the Canadian Forest Fire Weather Index System (Van Wagner, C.E.; Pickett, T.L. 1985). + * + * @param dmc Duff Moisture Code + * @param dc Drought Code + * + * @return A single Buildup Index (BUI) value + */ + +function buildupIndex(dmc: number, dc: number): number { + // Eq. 27a + let bui1 = (dmc === 0 && dc === 0) ? 0 : 0.8 * dc * dmc / (dmc + 0.4 * dc); + + // Eq. 27b - next 3 lines + const p = (dmc === 0) ? 0 : (dmc - bui1) / dmc; + const cc = 0.92 + Math.pow(0.0114 * dmc, 1.7); + let bui0 = dmc - cc * p; + + // Constraints + bui0 = Math.max(0, bui0); + bui1 = Math.min(dmc, bui1); + + return bui1; +} + +// Deprecated wrapper function +function buiCalc(...args: any[]): number { + console.warn(".Deprecated: 'buildupIndex' function is deprecated"); + return buildupIndex(args[0], args[1]); +} + +export { buildupIndex, buiCalc }; diff --git a/src/fwi/crown_base_height.ts b/src/fwi/crown_base_height.ts new file mode 100644 index 0000000..77cb428 --- /dev/null +++ b/src/fwi/crown_base_height.ts @@ -0,0 +1,36 @@ +/** + * Crown Base Height Calculator + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param CBH Crown Base Height + * @param SD Stand Density + * @param SH Stand Height + * + * @return CBH: Adjusted Crown Base Height + */ + +function crownBaseHeight(FUELTYPE: string, CBH: number, SD: number, SH: number): number { + // Default Crown Base Heights for different fuel types + const CBHs: { [key: string]: number } = { + "C1": 2, "C2": 3, "C3": 8, "C4": 4, "C5": 18, "C6": 7, "C7": 10, + "D1": 0, "M1": 6, "M2": 6, "M3": 6, "M4": 6, "S1": 0, "S2": 0, + "S3": 0, "O1A": 0, "O1B": 0 + }; + + // Adjust Crown Base Height based on conditions + if (CBH <= 0 || CBH > 50 || isNaN(CBH)) { + if (FUELTYPE === "C6" && SD > 0 && SH > 0) { + CBH = -11.2 + 1.06 * SH + 0.0017 * SD; + } else { + CBH = CBHs[FUELTYPE] !== undefined ? CBHs[FUELTYPE] : CBH; + } + } + + // Ensure CBH is not less than a very small positive number + CBH = CBH < 0 ? 1e-07 : CBH; + + return CBH; +} + +// Exporting function for external usage +export { crownBaseHeight }; diff --git a/src/fwi/crown_fuel_load.ts b/src/fwi/crown_fuel_load.ts new file mode 100644 index 0000000..e66249b --- /dev/null +++ b/src/fwi/crown_fuel_load.ts @@ -0,0 +1,27 @@ +/** + * Crown Fuel Load function + * + * @param FUELTYPE Character fuel type indicator + * @param CFL Crown Fuel Load + * + * @return Adjusted Crown Fuel Load + */ + +function crownFuelLoad(FUELTYPE: string, CFL: number): number { + // Default Crown Fuel Loads for different fuel types + const CFLs: { [key: string]: number } = { + "C1": 0.75, "C2": 0.8, "C3": 1.15, "C4": 1.2, "C5": 1.2, "C6": 1.8, "C7": 0.5, + "D1": 0, "M1": 0.8, "M2": 0.8, "M3": 0.8, "M4": 0.8, "S1": 0, "S2": 0, + "S3": 0, "O1A": 0, "O1B": 0 + }; + + // Adjust Crown Fuel Load based on conditions + if (CFL <= 0 || CFL > 2 || isNaN(CFL)) { + CFL = CFLs[FUELTYPE] !== undefined ? CFLs[FUELTYPE] : CFL; + } + + return CFL; +} + +// Exporting function for external usage +export { crownFuelLoad }; diff --git a/src/fwi/direction.ts b/src/fwi/direction.ts new file mode 100644 index 0000000..216c154 --- /dev/null +++ b/src/fwi/direction.ts @@ -0,0 +1,55 @@ +/** + * Direction definer + * + * @description New DIRECTION function to determine clockwise or + * counter-clockwise "interpretation" + * + * @param bearingT1T2 Bearing between T1 and T2 + * @param bearingT1T3 Bearing between T1 and T3 + * @param ThetaAdeg Direction + * + * @return DIR - a direction in degrees + */ + +function direction(bearingT1T2: number[], bearingT1T3: number[], ThetaAdeg: number): number[] { + const DIR: number[] = new Array(bearingT1T2.length).fill(NaN); + + for (let i = 0; i < bearingT1T2.length; i++) { + let T1T2 = bearingT1T2[i]; + let T1T3 = bearingT1T3[i]; + + if (T1T2 > 0 && T1T3 > 0 && T1T2 > T1T3) { + DIR[i] = T1T2 - ThetaAdeg; + } else if (T1T2 > 0 && T1T3 > 0 && T1T2 < T1T3) { + DIR[i] = T1T2 + ThetaAdeg; + } else if (T1T2 < 0 && T1T3 < 0 && T1T2 > T1T3) { + DIR[i] = T1T2 - ThetaAdeg; + } else if (T1T2 < 0 && T1T3 < 0 && T1T2 < T1T3) { + DIR[i] = T1T2 + ThetaAdeg; + } else if (T1T2 > 0 && T1T2 < 90 && T1T3 < 0 && T1T3 > -90) { + DIR[i] = T1T2 - ThetaAdeg; + } else if (T1T2 < 0 && T1T2 > -90 && T1T3 > 0 && T1T3 < 90) { + DIR[i] = T1T2 + ThetaAdeg; + } else if (T1T2 > 90 && T1T3 < -90 && T1T2 + ThetaAdeg > 180) { + DIR[i] = T1T2 + ThetaAdeg - 360; + } else if (T1T2 > 90 && T1T3 < -90 && T1T2 + ThetaAdeg < 180) { + DIR[i] = T1T2 + ThetaAdeg; + } else if (T1T2 < -90 && T1T3 > 90 && T1T2 - ThetaAdeg < -180) { + DIR[i] = T1T2 - ThetaAdeg + 360; + } else if (T1T2 < -90 && T1T3 > 90 && T1T2 - ThetaAdeg > -180) { + DIR[i] = T1T2 - ThetaAdeg; + } + + if (DIR[i] < -180) { + DIR[i] = DIR[i] + 10; + } + if (DIR[i] > 180) { + DIR[i] = DIR[i] - 10; + } + } + + return DIR; +} + +// Exporting function for external usage +export { direction }; diff --git a/src/fwi/distance_at_time.ts b/src/fwi/distance_at_time.ts new file mode 100644 index 0000000..08ac79e --- /dev/null +++ b/src/fwi/distance_at_time.ts @@ -0,0 +1,41 @@ +/** + * Distance at time t calculator + * + * @description Calculate the Head fire spread distance at time t. In the + * documentation this variable is just "D". + * + * All variables names are laid out in the same manner as Forestry Canada Fire + * Danger Group (FCFDG) (1992). Development and Structure of the Canadian + * Forest Fire Behavior Prediction System." Technical Report ST-X-3, + * Forestry Canada, Ottawa, Ontario. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param ROSeq The predicted equilibrium rate of spread (m/min) + * @param HR The elapsed time (min) + * @param CFB Crown Fraction Burned + * + * @return DISTt - Head fire spread distance at time t + */ + +function distanceAtTime(FUELTYPE: string, ROSeq: number, HR: number, CFB: number): number { + // Eq. 72 (FCFDG 1992) + // Calculate the alpha constant for the DISTt calculation + let alpha: number; + if (["C1", "O1A", "O1B", "S1", "S2", "S3", "D1"].includes(FUELTYPE)) { + alpha = 0.115; + } else { + alpha = 0.115 - 18.8 * Math.pow(CFB, 2.5) * Math.exp(-8 * CFB); + } + + // Eq. 71 (FCFDG 1992) Calculate Head fire spread distance + const DISTt = ROSeq * (HR + Math.exp(-alpha * HR) / alpha - 1 / alpha); + return DISTt; +} + +function DISTtcalc(...args: [string, number, number, number]): number { + console.warn("Deprecated: distanceAtTime"); + return distanceAtTime(...args); +} + +// Exporting functions for external usage +export { distanceAtTime, DISTtcalc }; diff --git a/src/fwi/drought_code.ts b/src/fwi/drought_code.ts new file mode 100644 index 0000000..420d3ab --- /dev/null +++ b/src/fwi/drought_code.ts @@ -0,0 +1,73 @@ +/** + * Drought Code Calculator + * + * @description Drought Code Calculation. All code is based on a C code library + * that was written by Canadian Forest Service Employees, which was originally + * based on the Fortran code listed in the reference below. + * All equations in this code refer to that document. Equations and FORTRAN + * program for the Canadian Forest Fire Weather Index System. 1985. Van Wagner, + * C.E.; Pickett, T.L. Canadian Forestry Service, Petawawa National Forestry + * Institute, Chalk River, Ontario. Forestry Technical Report 33. 18 p. + * Additional reference on FWI system Development and structure of the Canadian + * Forest Fire Weather Index System. 1987. Van Wagner, C.E. Canadian Forestry + * Service, Headquarters, Ottawa. Forestry Technical Report 35. 35 p. + * + * @param dc_yda The Drought Code from previous iteration + * @param temp Temperature (centigrade) + * @param rh Relative Humidity (%) + * @param prec Precipitation(mm) + * @param lat Latitude (decimal degrees) + * @param mon Month (1-12) + * @param latAdjust Latitude adjustment (TRUE, FALSE, default=TRUE) + * + * @return A single drought code value + */ + +function droughtCode(dc_yda: number, temp: number, rh: number, prec: number, lat: number, mon: number, latAdjust: boolean = true): number { + // Day length factor for DC Calculations + // 20N: North of 20 degrees N + const fl01: number[] = [-1.6, -1.6, -1.6, 0.9, 3.8, 5.8, 6.4, 5, 2.4, 0.4, -1.6, -1.6]; + // 20S: South of 20 degrees S + const fl02: number[] = [6.4, 5, 2.4, 0.4, -1.6, -1.6, -1.6, -1.6, -1.6, 0.9, 3.8, 5.8]; + // Constrain temperature + temp = temp < -2.8 ? -2.8 : temp; + + // Eq. 22 - Potential Evapotranspiration + let pe = (0.36 * (temp + 2.8) + fl01[mon - 1]) / 2; + + // Daylength factor adjustment by latitude for Potential Evapotranspiration + if (latAdjust) { + if (lat <= -20) { + pe = (0.36 * (temp + 2.8) + fl02[mon - 1]) / 2; + } else if (lat > -20 && lat <= 20) { + pe = (0.36 * (temp + 2.8) + 1.4) / 2; + } + } + + // Cap potential evapotranspiration at 0 for negative winter DC values + pe = pe < 0 ? 0 : pe; + + const ra = prec; + // Eq. 18 - Effective Rainfall + const rw = 0.83 * ra - 1.27; + // Eq. 19 + const smi = 800 * Math.exp(-1 * dc_yda / 400); + // Alteration to Eq. 21 + let dr0 = dc_yda - 400 * Math.log(1 + 3.937 * rw / smi); + dr0 = dr0 < 0 ? 0 : dr0; + // if precip is less than 2.8 then use yesterday's DC + const dr = prec <= 2.8 ? dc_yda : dr0; + // Alteration to Eq. 23 + let dc1 = dr + pe; + dc1 = dc1 < 0 ? 0 : dc1; + + return dc1; +} + +function dcCalc(...args: [number, number, number, number, number, number, boolean?]): number { + console.warn("Deprecated: droughtCode"); + return droughtCode(...args); +} + +// Exporting functions for external usage +export { droughtCode, dcCalc }; diff --git a/src/fwi/duff_moisture_code.ts b/src/fwi/duff_moisture_code.ts new file mode 100644 index 0000000..9b9aedc --- /dev/null +++ b/src/fwi/duff_moisture_code.ts @@ -0,0 +1,94 @@ +/** + * Duff Moisture Code Calculator + * + * @description Duff Moisture Code Calculation. All code is based on a C code + * library that was written by Canadian Forest Service Employees, which was + * originally based on the Fortran code listed in the reference below. All + * equations in this code refer to that document. + * + * Equations and FORTRAN program for the Canadian Forest Fire Weather Index + * System. 1985. Van Wagner, C.E.; Pickett, T.L. Canadian Forestry Service, + * Petawawa National Forestry Institute, Chalk River, Ontario. Forestry + * Technical Report 33. 18 p. + * + * Additional reference on FWI system + * + * Development and structure of the Canadian Forest Fire Weather Index System. + * 1987. Van Wagner, C.E. Canadian Forestry Service, Headquarters, Ottawa. + * Forestry Technical Report 35. 35 p. + * + * @param dmc_yda The Duff Moisture Code from previous iteration + * @param temp Temperature (centigrade) + * @param rh Relative Humidity (%) + * @param prec Precipitation(mm) + * @param lat Latitude (decimal degrees) + * @param mon Month (1-12) + * @param latAdjust Latitude adjustment (TRUE, FALSE, default=TRUE) + * + * @return A single drought moisture code value + */ + +function duffMoistureCode(dmc_yda: number, temp: number, rh: number, prec: number, lat: number, mon: number, latAdjust: boolean = true): number { + // Reference latitude for DMC day length adjustment + // 46N: Canadian standard, latitude >= 30N (Van Wagner 1987) + const ell01: number[] = [6.5, 7.5, 9, 12.8, 13.9, 13.9, 12.4, 10.9, 9.4, 8, 7, 6]; + // 20N: For 30 > latitude >= 10 + const ell02: number[] = [7.9, 8.4, 8.9, 9.5, 9.9, 10.2, 10.1, 9.7, 9.1, 8.6, 8.1, 7.8]; + // 20S: For -10 > latitude >= -30 + const ell03: number[] = [10.1, 9.6, 9.1, 8.5, 8.1, 7.8, 7.9, 8.3, 8.9, 9.4, 9.9, 10.2]; + // 40S: For -30 > latitude + const ell04: number[] = [11.5, 10.5, 9.2, 7.9, 6.8, 6.2, 6.5, 7.4, 8.7, 10, 11.2, 11.8]; + + // Constrain low end of temperature + temp = temp < -1.1 ? -1.1 : temp; + + // Eq. 16 - The log drying rate + let rk = 1.894 * (temp + 1.1) * (100 - rh) * ell01[mon - 1] * 1e-04; + + // Adjust the day length and thus the drying rate based on latitude and month + if (latAdjust) { + if (lat <= 30 && lat > 10) { + rk = 1.894 * (temp + 1.1) * (100 - rh) * ell02[mon - 1] * 1e-04; + } else if (lat <= -10 && lat > -30) { + rk = 1.894 * (temp + 1.1) * (100 - rh) * ell03[mon - 1] * 1e-04; + } else if (lat <= -30 && lat >= -90) { + rk = 1.894 * (temp + 1.1) * (100 - rh) * ell04[mon - 1] * 1e-04; + } else if (lat <= 10 && lat > -10) { + rk = 1.894 * (temp + 1.1) * (100 - rh) * 9 * 1e-04; + } + } + + // Constrain P + let pr = prec <= 1.5 ? dmc_yda : (() => { + const ra = prec; + // Eq. 11 - Net rain amount + const rw = 0.92 * ra - 1.27; + // Alteration to Eq. 12 to calculate more accurately + const wmi = 20 + 280 / Math.exp(0.023 * dmc_yda); + // Eqs. 13a, 13b, 13c + const b = dmc_yda <= 33 + ? 100 / (0.5 + 0.3 * dmc_yda) + : dmc_yda <= 65 + ? 14 - 1.3 * Math.log(dmc_yda) + : 6.2 * Math.log(dmc_yda) - 17.2; + // Eq. 14 - Moisture content after rain + const wmr = wmi + 1000 * rw / (48.77 + b * rw); + // Alteration to Eq. 15 to calculate more accurately + return 43.43 * (5.6348 - Math.log(wmr - 20)); + })(); + + pr = pr < 0 ? 0 : pr; + // Calculate final P (DMC) + let dmc1 = pr + rk; + dmc1 = dmc1 < 0 ? 0 : dmc1; + + return dmc1; +} + +function dmcCalc(...args: [number, number, number, number, number, number, boolean?]): number { + console.warn("Deprecated: duffMoistureCode"); + return duffMoistureCode(...args); +} + +// Exporting functions for external usage +export { duffMoistureCode, dmcCalc }; diff --git a/src/fwi/fbp.ts b/src/fwi/fbp.ts new file mode 100644 index 0000000..8c35f62 --- /dev/null +++ b/src/fwi/fbp.ts @@ -0,0 +1,409 @@ +/* #' Fire Behavior Prediction System function +#' +#' @description \code{fbp} calculates the outputs from the Canadian Forest Fire +#' Behavior Prediction (FBP) System (Forestry Canada Fire Danger Group 1992) +#' based on given fire weather and fuel moisture conditions (from the Canadian +#' Forest Fire Weather Index (FWI) System (Van Wagner 1987)), fuel type, date, +#' and slope. Fire weather, for the purpose of FBP System calculation, +#' comprises observations of 10 m wind speed and direction at the time of the +#' fire, and two associated outputs from the Fire Weather Index System, the +#' Fine Fuel Moisture Content (FFMC) and Buildup Index (BUI). FWI System +#' components can be calculated with the sister function \code{\link{fwi}}. +#' +#' The Canadian Forest Fire Behavior Prediction (FBP) System (Forestry Canada +#' Fire Danger Group 1992) is a subsystem of the Canadian Forest Fire Danger +#' Rating System, which also includes the Canadian Forest Fire Weather Index +#' (FWI) System. The FBP System provides quantitative estimates of head fire +#' spread rate, fuel consumption, fire intensity, and a basic fire description +#' (e.g., surface, crown) for 16 different important forest and rangeland types +#' across Canada. Using a simple conceptual model of the growth of a point +#' ignition as an ellipse through uniform fuels and under uniform weather +#' conditions, the system gives, as a set of secondary outputs, estimates of +#' flank and back fire behavior and consequently fire area perimeter length and +#' growth rate. +#' +#' The FBP System evolved since the mid-1970s from a series of regionally +#' developed burning indexes to an interim edition of the nationally develop +#' FBP system issued in 1984. Fire behavior models for spread rate and fuel +#' consumption were derived from a database of over 400 experimental, wild and +#' prescribed fire observations. The FBP System, while providing quantitative +#' predictions of expected fire behavior is intended to supplement the +#' experience and judgment of operational fire managers (Hirsch 1996). +#' +#' The FBP System was updated with some minor corrections and revisions in 2009 +#' (Wotton et al. 2009) with several additional equations that were initially +#' not included in the system. This fbp function included these updates and +#' corrections to the original equations and provides a complete suite of fire +#' behavior prediction variables. Default values of optional input variables +#' provide a reasonable mid-range setting. Latitude, longitude, elevation, and +#' the date are used to calculate foliar moisture content, using a set of +#' models defined in the FBP System; note that this latitude/longitude-based +#' function is only valid for Canada. If the Foliar Moisture Content (FMC) is +#' specified directly as an input, the fbp function will use this value +#' directly rather than calculate it. This is also true of other input +#' variables. +#' +#' Note that Wind Direction (WD) is the compass direction from which wind is +#' coming. Wind azimuth (not an input) is the direction the wind is blowing to +#' and is 180 degrees from wind direction; in the absence of slope, the wind +#' azimuth is coincident with the direction the head fire will travel (the +#' spread direction azimuth, RAZ). Slope aspect is the main compass direction +#' the slope is facing. Slope azimuth (not an input) is the direction a head +#' fire will spread up slope (in the absence of wind effects) and is 180 +#' degrees from slope aspect (Aspect). Wind direction and slope aspect are the +#' commonly used directional identifiers when specifying wind and slope +#' orientation respectively. The input theta specifies an angle (given as a +#' compass bearing) at which a user is interested in fire behavior predictions; +#' it is typically some angle off of the final spread rate direction since if +#' for instance theta=RAZ (the final spread azimuth of the fire) then the rate +#' of spread at angle theta (TROS) will be equivalent to ROS. +#' +#' @param input The input data, a data.frame containing fuel types, fire +#' weather component, and slope (see below). Each vector of inputs defines a +#' single FBP System prediction for a single fuel type and set of weather +#' conditions. The data.frame can be used to evaluate the FBP System for a +#' single fuel type and instant in time, or multiple records for a single point +#' (e.g., one weather station, either hourly or daily for instance) or multiple +#' points (multiple weather stations or a gridded surface). All input variables +#' have to be named as listed below, but they are case insensitive, and do not +#' have to be in any particular order. Fuel type is of type character; other +#' arguments are numeric. Missing values in numeric variables could either be +#' assigned as NA or leave as blank.\cr\cr +#' +#' +#' \tabular{lll}{ +#' \bold{Required Inputs:}\tab\tab\cr +#' \bold{Input} \tab \bold{Description/Full name} \tab \bold{Defaults}\cr +#' +#' \var{FuelType} +#' \tab FBP System Fuel Type including "C-1",\cr +#' \tab"C-2", "C-3", "C-4","C-5", "C-6", "C-7",\cr +#' \tab "D-1", "M-1", "M-2", "M-3", "M-4", "NF",\cr +#' \tab "D-1", "S-2", "S-3", "O-1a", "O-1b", and\cr +#' \tab "WA", where "WA" and "NF" stand for \cr +#' \tab "water" and "non-fuel", respectively.\cr\cr +#' +#' \var{LAT} \tab Latitude [decimal degrees] \tab 55\cr +#' \var{LONG} \tab Longitude [decimal degrees] \tab -120\cr +#' \var{FFMC} \tab Fine fuel moisture code [FWI System component] \tab 90\cr +#' \var{BUI} \tab Buildup index [FWI System component] \tab 60\cr +#' \var{WS} \tab Wind speed [km/h] \tab 10\cr +#' \var{GS} \tab Ground Slope [percent] \tab 0\cr +#' \var{Dj} \tab Julian day \tab 180\cr +#' \var{Aspect} \tab Aspect of the slope [decimal degrees] \tab 0\cr\cr +#' +#' \bold{Optional Inputs (1):} +#' \tab Variables associated with certain fuel \cr +#' \tab types. These could be skipped if relevant \cr +#' \tab fuel types do not appear in the input data.\cr\cr +#' +#' \bold{Input} \tab \bold{Full names of inputs} \tab \bold{Defaults}\cr +#' +#' \var{PC} \tab Percent Conifer for M1/M2 [percent] \tab 50\cr +#' \var{PDF} \tab Percent Dead Fir for M3/M4 [percent] \tab 35\cr +#' \var{cc} \tab Percent Cured for O1a/O1b [percent] \tab 80\cr +#' \var{GFL} \tab Grass Fuel Load [kg/m^2] \tab 0.35\cr\cr +#' +#' \bold{Optional Inputs (2):} +#' \tab Variables that could be ignored without \cr +#' \tab causing major impacts to the primary outputs\cr\cr +#' +#' \bold{Input} \tab \bold{Full names of inputs} \tab \bold{Defaults}\cr +#' \var{CBH} \tab Crown to Base Height [m] \tab 3\cr +#' \var{WD} \tab Wind direction [decimal degrees] \tab 0\cr +#' \var{Accel} \tab Acceleration: 1 = point, 0 = line \tab 0\cr +#' \var{ELV*} \tab Elevation [meters above sea level] \tab NA\cr +#' \var{BUIEff}\tab Buildup Index effect: 1=yes, 0=no \tab 1\cr +#' \var{D0} \tab Julian day of minimum Foliar Moisture Content \tab 0\cr +#' \var{hr} \tab Hours since ignition \tab 1\cr +#' \var{ISI} \tab Initial spread index \tab 0\cr +#' \var{CFL} \tab Crown Fuel Load [kg/m^2]\tab 1.0\cr +#' \var{FMC} \tab Foliar Moisture Content if known [percent] \tab 0\cr +#' \var{SH} \tab C-6 Fuel Type Stand Height [m] \tab 0\cr +#' \var{SD} \tab C-6 Fuel Type Stand Density [stems/ha] \tab 0\cr +#' \var{theta} \tab Elliptical direction of calculation [degrees] \tab 0\cr\cr } +#' +#' @param output FBP output offers 3 options (see details in \bold{Values} +#' section): +#' +#' \tabular{lc}{ \bold{Outputs} \tab \bold{Number of outputs}\cr +#' \var{Primary(\bold{default})} \tab 8\cr +#' \var{Secondary} \tab 34\cr +#' \var{All} \tab 42\cr\cr} +#' +#' @param m Optimal number of pixels at each iteration of computation when +#' \code{nrow(input) >= 1000}. Default \code{m = NULL}, where the function will +#' assign \code{m = 1000} when \code{nrow(input)} is between 1000 and 500,000, +#' and \code{m = 3000} otherwise. By including this option, the function is +#' able to process large dataset more efficiently. The optimal value may vary +#' with different computers. +#' +#' @param cores Number of CPU cores (integer) used in the computation, default +#' is 1. By signing \code{cores > 1}, the function will apply parallel +#' computation technique provided by the \code{foreach} package, which +#' significantly reduces the computation time for large input data (over a +#' million records). For small dataset, \code{cores=1} is actually faster. +#' +#' @return \code{fbp} returns a dataframe with primary, secondary, or all +#' output variables, a combination of the primary and secondary outputs. +#' +#' \bold{Primary} FBP output includes the following 8 variables: +#' +#' \item{CFB}{Crown Fraction Burned by the head fire} +#' \item{CFC}{Crown Fuel Consumption [kg/m^2]} +#' \item{FD}{Fire description (1=Surface, 2=Intermittent, 3=Crown)} +#' \item{HFI}{Head Fire Intensity [kW/m]} +#' \item{RAZ}{Spread direction azimuth [degrees]} +#' \item{ROS}{Equilibrium Head Fire Rate of Spread [m/min]} +#' \item{SFC}{Surface Fuel Consumption [kg/m^2]} +#' \item{TFC}{Total Fuel Consumption [kg/m^2]} +#' +#' \bold{Secondary} FBP System outputs include the following 34 raster layers. +#' In order to calculate the reliable secondary outputs, depending on the +#' outputs, optional inputs may have to be provided. +#' +#' \item{BE}{BUI effect on spread rate} +#' \item{SF}{Slope Factor (multiplier for ROS increase upslope)} +#' \item{ISI}{Initial Spread Index} +#' \item{FFMC}{Fine fuel moisture code [FWI System component]} +#' \item{FMC}{Foliar Moisture Content [\%]} +#' \item{Do}{Julian Date of minimum FMC} +#' \item{RSO}{Critical spread rate for crowning [m/min]} +#' \item{CSI}{Critical Surface Intensity for crowning [kW/m]} +#' \item{FROS}{Equilibrium Flank Fire Rate of Spread [m/min]} +#' \item{BROS}{Equilibrium Back Fire Rate of Spread [m/min]} +#' \item{HROSt}{Head Fire Rate of Spread at time hr [m/min]} +#' \item{FROSt}{Flank Fire Rate of Spread at time hr [m/min]} +#' \item{BROSt}{Back Fire Rate of Spread at time hr [m/min]} +#' \item{FCFB}{Flank Fire Crown Fraction Burned} +#' \item{BCFB}{Back Fire Crown Fraction Burned} +#' \item{FFI}{Equilibrium Spread Flank Fire Intensity [kW/m]} +#' \item{BFI}{Equilibrium Spread Back Fire Intensity [kW/m]} +#' \item{FTFC}{Flank Fire Total Fuel Consumption [kg/m^2] } +#' \item{BTFC}{Back Fire Total Fuel Consumption [kg/m^2] } +#' \item{DH}{Head Fire Spread Distance after time hr [m] } +#' \item{DB}{Back Fire Spread Distance after time hr [m] } +#' \item{DF}{Flank Fire Spread Distance after time hr [m] } +#' \item{TI}{Time to Crown Fire Initiation [hrs since ignition] } +#' \item{FTI}{Time to Flank Fire Crown initiation [hrs since ignition]} +#' \item{BTI}{Time to Back Fire Crown initiation [hrs since ignition]} +#' \item{LB}{Length to Breadth ratio} +#' \item{LBt}{Length to Breadth ratio after elapsed time hr } +#' \item{WSV}{Net vectored wind speed [km/hr]} +#' \item{TROS*}{Equilibrium Rate of Spread at bearing theta [m/min] } +#' \item{TROSt*}{Rate of Spread at bearing theta at time t [m/min] } +#' \item{TCFB*}{Crown Fraction Burned at bearing theta } +#' \item{TFI*}{Fire Intensity at bearing theta [kW/m] } +#' \item{TTFC*}{Total Fuel Consumption at bearing theta [kg/m^2] } +#' \item{TTI*}{Time to Crown Fire initiation at bearing theta [hrs since +#' ignition] } +#' +#' *These outputs represent fire behaviour at a point on the perimeter of an +#' elliptical fire defined by a user input angle theta. theta represents the +#' bearing of a line running between the fire ignition point and a point on the +#' perimeter of the fire. It is important to note that in this formulation the +#' theta is a bearing and does not represent the angle from the semi-major axis +#' (spread direction) of the ellipse. This formulation is similar but not +#' identical to methods presented in Wotton et al (2009) and Tymstra et al +#' (2009). +#' @author Xianli Wang, Alan Cantin, Marc-André Parisien, Mike Wotton, Kerry +#' Anderson, and Mike Flannigan +#' @seealso \code{\link{fwi}, \link{fbpRaster}} +#' @references 1. Hirsch K.G. 1996. Canadian Forest Fire Behavior Prediction +#' (FBP) System: user's guide. Nat. Resour. Can., Can. For. Serv., Northwest +#' Reg., North. For. Cent., Edmonton, Alberta. Spec. Rep. 7. 122p. +#' +#' 2. Forestry Canada Fire Danger Group. 1992. Development and structure of +#' the Canadian Forest Fire Behavior Prediction System. Forestry Canada, +#' Ottawa, Ontario Information Report ST-X-3. 63 p. +#' \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/10068.pdf} +#' +#' 3. Wotton, B.M., Alexander, M.E., Taylor, S.W. 2009. Updates and revisions +#' to the 1992 Canadian forest fire behavior prediction system. Nat. Resour. +#' Can., Can. For. Serv., Great Lakes For. Cent., Sault Ste. Marie, Ontario, +#' Canada. Information Report GLC-X-10, 45p. +#' \url{ +#' https://publications.gc.ca/collections/collection_2010/nrcan/ +#' Fo123-2-10-2009-eng.pdf} +#' +#' 4. Tymstra, C., Bryce, R.W., Wotton, B.M., Armitage, O.B. 2009. Development +#' and structure of Prometheus: the Canadian wildland fire growth simulation +#' Model. Nat. Resour. Can., Can. For. Serv., North. For. Cent., Edmonton, AB. +#' Inf. Rep. NOR-X-417. +#' \url{https://d1ied5g1xfgpx8.cloudfront.net/pdfs/31775.pdf} +#' +#' @importFrom parallel makeCluster stopCluster +#' @importFrom doParallel registerDoParallel +#' @importFrom foreach foreach registerDoSEQ %dopar% +#' @importFrom data.table rbindlist setkey +#' +#' @keywords methods +#' @examples +#' +#' library(cffdrs) +#' # The dataset is the standard test data for FPB system +#' # provided by Wotton et al (2009) +#' data("test_fbp") +#' head(test_fbp) +#' # id FuelType LAT LONG ELV FFMC BUI WS WD GS Dj D0 hr +#' # 1 1 C-1 55 110 NA 90 130 20.0 0 15 182 NA 0.33333333 +#' # 2 2 C2 50 90 NA 97 119 20.4 0 75 121 NA 0.33333333 +#' # 3 3 C-3 55 110 NA 95 30 50.0 0 0 182 NA 0.08333333 +#' # 4 4 C-4 55 105 200 85 82 0.0 NA 75 182 NA 0.50000000 +#' # 5 5 c5 55 105 NA 88 56 3.4 0 23 152 145 0.50000000 +#' # id PC PDF GFL cc theta Accel Aspect BUIEff CBH CFL ISI +#' # 1 NA NA NA NA 0 1 270 1 NA NA 0 +#' # 2 NA NA NA NA 0 1 315 1 NA NA 0 +#' # 3 NA NA NA NA 0 1 180 1 NA NA 0 +#' # 4 NA NA NA NA 0 1 315 1 NA NA 0 +#' # 5 NA NA NA NA 0 1 180 1 NA NA 0 +#' +#' # Primary output (default) +#' fbp(test_fbp) +#' # or +#' fbp(test_fbp, output = "Primary") +#' # or +#' fbp(test_fbp, "P") +#' # Secondary output +#' fbp(test_fbp, "Secondary") +#' # or +#' fbp(test_fbp, "S") +#' # All output +#' fbp(test_fbp, "All") +#' # or +#' fbp(test_fbp, "A") +#' # For a single record: +#' fbp(test_fbp[7, ]) +#' # For a section of the records: +#' fbp(test_fbp[8:13, ]) +#' # fbp function produces the default values if no data is fed to +#' # the function: +#' fbp() +#' +#' @export fbp */ + +/** + * Fire Behavior Prediction System function + * + * @param input The input data, an array of objects containing fuel types, fire weather component, and slope. + * @param output FBP output options: "Primary", "Secondary", or "All". + * @param m Optimal number of pixels at each iteration of computation. + * @param cores Number of CPU cores used in the computation. + * + * @return An array of objects with primary, secondary, or all output variables. + */ + +async function fbp( + input: Record[] | null = null, + output: string = "Primary", + m: number | null = null, + cores: number = 1 + ): Promise[]> { + let ID: any = null; + + // If input is not provided, then calculate FBP with default values + if (input === null) { + input = [{ + FUELTYPE: "C2", + ACCEL: 0, + DJ: 180, + D0: 0, + ELV: 0, + BUIEFF: 1, + HR: 1, + FFMC: 90, + ISI: 0, + BUI: 60, + WS: 10, + WD: 0, + GS: 0, + ASPECT: 0, + PC: 50, + PDF: 35, + CC: 80, + GFL: 0.35, + CBH: 3, + CFL: 1, + LAT: 55, + LONG: -120, + FMC: 0, + THETA: 0 + }]; + input = input.map(record => ({ ...record, FUELTYPE: record.FUELTYPE.toString() })); + return fireBehaviourPrediction(input); + } else { + input = input.map(record => { + let newRecord: Record = {}; + for (let key in record) { + newRecord[key.toUpperCase()] = record[key]; + } + return newRecord; + }); + + // Determine optimal number of pixels to process at each iteration + if (m === null) { + m = input.length > 500000 ? 3000 : 1000; + } + m = input.length >= m ? m : input.length; + const n0 = Math.round(input.length / m); + const n = m * n0 >= input.length ? n0 : n0 + 1; + + // Set up parallel processing, if # of cores is entered + if (cores > 1) { + // Create and register a set of parallel instances for foreach + const ca = await Promise.all(Array.from({ length: n }, async (_, i) => { + if (i + 1 === n) { + return fireBehaviourPrediction(input!.slice((i * m!), input!.length), output); + } else { + return fireBehaviourPrediction(input!.slice((i * m!), (i + 1) * m!), output); + } + })); + + // Flatten the array of arrays + const fullList = ca.flat(); + + // Create a single keyed data table + const keyedList: { [key: string]: any } = fullList.reduce((acc, item) => { + acc[item.ID] = item; + return acc; + }, {}); + + // Convert to array to keep backwards compatibility + return Object.values(keyedList); + } else { + const ca: Record[][] = []; + for (let i = 0; i < n; i++) { + let foo; + if (i + 1 === n) { + foo = input.slice((i * m!), input.length); + } else { + foo = input.slice((i * m!), (i + 1) * m!); + } + ca.push(await fireBehaviourPrediction(foo, output)); + } + + // Flatten the array of arrays + const fullList = ca.flat(); + + // Create a single keyed data table + const keyedList: { [key: string]: any } = fullList.reduce((acc, item) => { + acc[item.ID] = item; + return acc; + }, {}); + + // Convert to array to keep backwards compatibility + return Object.values(keyedList); + } + } + } + + // Dummy function to simulate fire_behaviour_prediction in TypeScript + async function fireBehaviourPrediction(input: Record[], output: string = "Primary"): Promise[]> { + // This function should be implemented based on actual logic of fire_behaviour_prediction + // For now, it just returns the input array for demonstration purposes + return input; + } + + // Exporting the function for external usage + export { fbp }; + \ No newline at end of file diff --git a/src/fwi/fbpRaster.ts b/src/fwi/fbpRaster.ts new file mode 100644 index 0000000..bb6b34a --- /dev/null +++ b/src/fwi/fbpRaster.ts @@ -0,0 +1,459 @@ +/* #' Raster-based Fire Behavior Prediction System Calculations +#' +#' @description \code{fbpRaster} calculates the outputs from the Canadian Forest +#' Fire Behavior Prediction (FBP) System (Forestry Canada Fire Danger Group +#' 1992) based on raster format fire weather and fuel moisture conditions (from +#' the Canadian Forest Fire Weather Index (FWI) System (Van Wagner 1987)), fuel +#' type, date, and slope. Fire weather, for the purpose of FBP System +#' calculation, comprises observations of 10 m wind speed and direction at the +#' time of the fire, and two associated outputs from the Fire Weather Index +#' System, the Fine Fuel Moisture Content (FFMC) and Buildup Index (BUI). +#' Raster-based FWI System components can be calculated with the sister +#' function \code{\link{fwiRaster}}. +#' +#' The Canadian Forest Fire Behavior Prediction (FBP) System (Forestry Canada +#' Fire Danger Group 1992) is a subsystem of the Canadian Forest Fire Danger +#' Rating System, which also includes the Canadian Forest Fire Weather Index +#' (FWI) System. The FBP System provides quantitative estimates of head fire +#' spread rate, fuel consumption, fire intensity, and a basic fire description +#' (e.g., surface, crown) for 16 different important forest and rangeland types +#' across Canada. Using a simple conceptual model of the growth of a point +#' ignition as an ellipse through uniform fuels and under uniform weather +#' conditions, the system gives, as a set of secondary outputs, estimates of +#' flank and back fire behavior and consequently fire area perimeter length and +#' growth rate. +#' +#' The FBP System evolved since the mid-1970s from a series of regionally +#' developed burning indexes to an interim edition of the nationally develop +#' FBP system issued in 1984. Fire behavior models for spread rate and fuel +#' consumption were derived from a database of over 400 experimental, wild and +#' prescribed fire observations. The FBP System, while providing quantitative +#' predictions of expected fire behavior is intended to supplement the +#' experience and judgment of operational fire managers (Hirsch 1996). +#' +#' The FBP System was updated with some minor corrections and revisions in 2009 +#' (Wotton et al. 2009) with several additional equations that were initially +#' not included in the system. This fbp function included these updates and +#' corrections to the original equations and provides a complete suite of fire +#' behavior prediction variables. Default values of optional input variables +#' provide a reasonable mid-range setting. Latitude, longitude, elevation, and +#' the date are used to calculate foliar moisture content, using a set of +#' models defined in the FBP System; note that this latitude/longitude-based +#' function is only valid for Canada. If the Foliar Moisture Content (FMC) is +#' specified directly as an input, the fbp function will use this value +#' directly rather than calculate it. This is also true of other input +#' variables. +#' +#' Note that Wind Direction (WD) is the compass direction from which wind is +#' coming. Wind azimuth (not an input) is the direction the wind is blowing to +#' and is 180 degrees from wind direction; in the absence of slope, the wind +#' azimuth is coincident with the direction the head fire will travel (the +#' spread direction azimuth, RAZ). Slope aspect is the main compass direction +#' the slope is facing. Slope azimuth (not an input) is the direction a head +#' fire will spread up slope (in the absence of wind effects) and is 180 +#' degrees from slope aspect (Aspect). Wind direction and slope aspect are the +#' commonly used directional identifiers when specifying wind and slope +#' orientation respectively. The input theta specifies an angle (given as a +#' compass bearing) at which a user is interested in fire behavior predictions; +#' it is typically some angle off of the final spread rate direction since if +#' for instance theta=RAZ (the final spread azimuth of the fire) then the rate +#' of spread at angle theta (TROS) will be equivalent to ROS. +#' +#' Because raster format data cannot hold characters, we have to code these fuel +#' types in numeric codes. In sequence, the codes are c(1:19). FuelType could +#' also be converted as factor and assigned to the raster layer, the function +#' will still work. +#' +#' \tabular{ll}{ +#' \bold{Fuel Type} \tab \bold{code} \cr +#' \verb{C-1} \tab 1 \cr +#' \verb{C-2} \tab 2 \cr +#' \verb{C-3} \tab 3 \cr +#' \verb{C-4} \tab 4 \cr +#' \verb{C-5} \tab 5 \cr +#' \verb{C-6} \tab 6 \cr +#' \verb{C-7} \tab 7 \cr +#' \verb{D-1} \tab 8 \cr +#' \verb{M-1} \tab 9 \cr +#' \verb{M-2} \tab 10 \cr +#' \verb{M-3} \tab 11 \cr +#' \verb{M-4} \tab 12 \cr +#' \verb{NF} \tab 13 \cr +#' \verb{O-1a} \tab 14 \cr +#' \verb{O-1b} \tab 15 \cr +#' \verb{S-1} \tab 16 \cr +#' \verb{S-2} \tab 17 \cr +#' \verb{S-3} \tab 18 \cr +#' \verb{WA} \tab 19 \cr\cr} +#' +#' @param input The input data, a RasterStack containing fuel types, fire +#' weather component, and slope layers (see below). Each vector of inputs +#' defines a single FBP System prediction for a single fuel type and set of +#' weather conditions. The RasterStack can be used to evaluate the FBP System +#' for a single fuel type and instant in time, or multiple records for a single +#' point (e.g., one weather station, either hourly or daily for instance) or +#' multiple points (multiple weather stations or a gridded surface). All input +#' variables have to be named as listed below, but they are case insensitive, +#' and do not have to be in any particular order. Fuel type is of type +#' character; other arguments are numeric. Missing values in numeric variables +#' could either be assigned as NA or leave as blank. +#' +#' \tabular{lll}{ +#' \bold{Required Inputs:}\tab\tab\cr +#' \bold{Input} \tab \bold{Description/Full name} \tab \bold{Defaults}\cr +#' +#' \var{FuelType} +#' \tab FBP System Fuel Type including "C-1",\cr +#' \tab"C-2", "C-3", "C-4","C-5", "C-6", "C-7",\cr +#' \tab "D-1", "M-1", "M-2", "M-3", "M-4", "NF",\cr +#' \tab "D-1", "S-2", "S-3", "O-1a", "O-1b", and\cr +#' \tab "WA", where "WA" and "NF" stand for \cr +#' \tab "water" and "non-fuel", respectively.\cr\cr +#' +#' \var{LAT} \tab Latitude [decimal degrees] \tab 55\cr +#' \var{LONG} \tab Longitude [decimal degrees] \tab -120\cr +#' \var{FFMC} \tab Fine fuel moisture code [FWI System component] \tab 90\cr +#' \var{BUI} \tab Buildup index [FWI System component] \tab 60\cr +#' \var{WS} \tab Wind speed [km/h] \tab 10\cr +#' \var{GS} \tab Ground Slope [percent] \tab 0\cr +#' \var{Dj} \tab Julian day \tab 180\cr +#' \var{Aspect} \tab Aspect of the slope [decimal degrees] \tab 0\cr\cr +#' +#' \bold{Optional Inputs (1):} +#' \tab Variables associated with certain fuel \cr +#' \tab types. These could be skipped if relevant \cr +#' \tab fuel types do not appear in the input data.\cr\cr +#' +#' \bold{Input} \tab \bold{Full names of inputs} \tab \bold{Defaults}\cr +#' +#' \var{PC} \tab Percent Conifer for M1/M2 [percent] \tab 50\cr +#' \var{PDF} \tab Percent Dead Fir for M3/M4 [percent] \tab 35\cr +#' \var{cc} \tab Percent Cured for O1a/O1b [percent] \tab 80\cr +#' \var{GFL} \tab Grass Fuel Load [kg/m^2] \tab 0.35\cr\cr +#' +#' \bold{Optional Inputs (2):} +#' \tab Variables that could be ignored without \cr +#' \tab causing major impacts to the primary outputs\cr\cr +#' +#' \bold{Input} \tab \bold{Full names of inputs} \tab \bold{Defaults}\cr +#' \var{CBH} \tab Crown to Base Height [m] \tab 3\cr +#' \var{WD} \tab Wind direction [decimal degrees] \tab 0\cr +#' \var{Accel} \tab Acceleration: 1 = point, 0 = line \tab 0\cr +#' \var{ELV*} \tab Elevation [meters above sea level] \tab NA\cr +#' \var{BUIEff}\tab Buildup Index effect: 1=yes, 0=no \tab 1\cr +#' \var{D0} \tab Julian day of minimum Foliar Moisture Content \tab 0\cr +#' \var{hr} \tab Hours since ignition \tab 1\cr +#' \var{ISI} \tab Initial spread index \tab 0\cr +#' \var{CFL} \tab Crown Fuel Load [kg/m^2]\tab 1.0\cr +#' \var{FMC} \tab Foliar Moisture Content if known [percent] \tab 0\cr +#' \var{SH} \tab C-6 Fuel Type Stand Height [m] \tab 0\cr +#' \var{SD} \tab C-6 Fuel Type Stand Density [stems/ha] \tab 0\cr +#' \var{theta} \tab Elliptical direction of calculation [degrees] \tab 0\cr\cr } +#' @param output FBP output offers 3 options (see details in \bold{Values} +#' section): +#' +#' \tabular{lc}{ \bold{Outputs} \tab \bold{Number of outputs}\cr \var{Primary +#' (\bold{default})} \tab 8\cr \var{Secondary} \tab 34\cr \var{All} \tab 42\cr +#' } +#' @param select Selected outputs +#' +#' @param m Optimal number of pixels at each iteration of computation when +#' \code{ncell(input) >= 1000}. Default m = NULL, where the function will +#' assign m = 1000 when \code{ncell(input)} is between 1000 and 500,000, and +#' m=3000 otherwise. By including this option, the function is able to process +#' large dataset more efficiently. The optimal value may vary with different +#' computers. +#' +#' @param cores Number of CPU cores (integer) used in the computation, default +#' is 1. By signing \code{cores > 1}, the function will apply parallel +#' computation technique provided by the \code{foreach} package, which +#' significantly reduces the computation time for large input data (over a +#' million grid points). For small dataset, \code{cores=1} is actually faster. +#' +#' @return \code{fbpRaster} returns a RasterStack with primary, secondary, or +#' all output variables, a combination of the primary and secondary outputs. +#' Primary FBP output includes the following 8 raster layers: +#' +#' \item{CFB}{Crown Fraction Burned by the head fire} +#' \item{CFC}{Crown Fuel Consumption [kg/m^2]} +#' \item{FD}{Fire description (1=Surface, 2=Intermittent, 3=Crown)} +#' \item{HFI}{Head Fire Intensity [kW/m]} +#' \item{RAZ}{Spread direction azimuth [degrees]} +#' \item{ROS}{Equilibrium Head Fire Rate of Spread [m/min]} +#' \item{SFC}{Surface Fuel Consumption [kg/m^2]} +#' \item{TFC}{Total Fuel Consumption [kg/m^2]} +#' +#' Secondary FBP System outputs include the following 34 raster layers. In order +#' to calculate the reliable secondary outputs, depending on the outputs, +#' optional inputs may have to be provided. +#' +#' \item{BE}{BUI effect on spread rate} +#' \item{SF}{Slope Factor (multiplier for ROS increase upslope)} +#' \item{ISI}{Initial Spread Index} +#' \item{FFMC}{Fine fuel moisture code [FWI System component]} +#' \item{FMC}{Foliar Moisture Content [\%]} +#' \item{Do}{Julian Date of minimum FMC} +#' \item{RSO}{Critical spread rate for crowning [m/min]} +#' \item{CSI}{Critical Surface Intensity for crowning [kW/m]} +#' \item{FROS}{Equilibrium Flank Fire Rate of Spread [m/min]} +#' \item{BROS}{Equilibrium Back Fire Rate of Spread [m/min]} +#' \item{HROSt}{Head Fire Rate of Spread at time hr [m/min]} +#' \item{FROSt}{Flank Fire Rate of Spread at time hr [m/min]} +#' \item{BROSt}{Back Fire Rate of Spread at time hr [m/min]} +#' \item{FCFB}{Flank Fire Crown Fraction Burned} +#' \item{BCFB}{Back Fire Crown Fraction Burned} +#' \item{FFI}{Equilibrium Spread Flank Fire Intensity [kW/m]} +#' \item{BFI}{Equilibrium Spread Back Fire Intensity [kW/m]} +#' \item{FTFC}{Flank Fire Total Fuel Consumption [kg/m^2] } +#' \item{BTFC}{Back Fire Total Fuel Consumption [kg/m^2] } +#' \item{DH}{Head Fire Spread Distance after time hr [m] } +#' \item{DB}{Back Fire Spread Distance after time hr [m] } +#' \item{DF}{Flank Fire Spread Distance after time hr [m] } +#' \item{TI}{Time to Crown Fire Initiation [hrs since ignition] } +#' \item{FTI}{Time to Flank Fire Crown initiation [hrs since ignition]} +#' \item{BTI}{Time to Back Fire Crown initiation [hrs since ignition]} +#' \item{LB}{Length to Breadth ratio} +#' \item{LBt}{Length to Breadth ratio after elapsed time hr } +#' \item{WSV}{Net vectored wind speed [km/hr]} +#' \item{TROS*}{Equilibrium Rate of Spread at bearing theta [m/min] } +#' \item{TROSt*}{Rate of Spread at bearing theta at time t [m/min] } +#' \item{TCFB*}{Crown Fraction Burned at bearing theta } +#' \item{TFI*}{Fire Intensity at bearing theta [kW/m] } +#' \item{TTFC*}{Total Fuel Consumption at bearing theta [kg/m^2] } +#' \item{TTI*}{Time to Crown Fire initiation at bearing theta [hrs since +#' ignition] } +#' +#' *These outputs represent fire behaviour at a point on the perimeter of an +#' elliptical fire defined by a user input angle theta. theta represents the +#' bearing of a line running between the fire ignition point and a point on the +#' perimeter of the fire. It is important to note that in this formulation the +#' theta is a bearing and does not represent the angle from the semi-major axis +#' (spread direction) of the ellipse. This formulation is similar but not +#' identical to methods presented in Wotton et al (2009) and Tymstra et al +#' (2009). +#' +#' @author Xianli Wang, Alan Cantin, Marc-André Parisien, Mike Wotton, Kerry +#' Anderson, and Mike Flannigan +#' @seealso \code{\link{fbp}, \link{fwiRaster}, \link{hffmcRaster}} +#' @references 1. Hirsch K.G. 1996. Canadian Forest Fire Behavior Prediction +#' (FBP) System: user's guide. Nat. Resour. Can., Can. For. Serv., Northwest +#' Reg., North. For. Cent., Edmonton, Alberta. Spec. Rep. 7. 122p. +#' +#' 2. Forestry Canada Fire Danger Group. 1992. Development and structure of +#' the Canadian Forest Fire Behavior Prediction System. Forestry Canada, +#' Ottawa, Ontario Information Report ST-X-3. 63 p. +#' \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/10068.pdf} +#' +#' 3. Wotton, B.M., Alexander, M.E., Taylor, S.W. 2009. Updates and revisions +#' to the 1992 Canadian forest fire behavior prediction system. Nat. Resour. +#' Can., Can. For. Serv., Great Lakes For. Cent., Sault Ste. Marie, Ontario, +#' Canada. Information Report GLC-X-10, 45p. +#' \url{ +#' https://publications.gc.ca/collections/collection_2010/nrcan/ +#' Fo123-2-10-2009-eng.pdf} +#' +#' 4. Tymstra, C., Bryce, R.W., Wotton, B.M., Armitage, O.B. 2009. Development +#' and structure of Prometheus: the Canadian wildland fire growth simulation +#' Model. Nat. Resour. Can., Can. For. Serv., North. For. Cent., Edmonton, AB. +#' Inf. Rep. NOR-X-417. +#' \url{https://d1ied5g1xfgpx8.cloudfront.net/pdfs/31775.pdf} +#' +#' @importFrom foreach registerDoSEQ +#' @importFrom terra rast ncell values setValues writeRaster +#' @importFrom methods is +#' @import sf +#' @importFrom data.table as.data.table data.table := +#' +#' @keywords methods +#' @examples +#' +#' # The dataset is the standard test data for FBP system +#' # provided by Wotton et al (2009), and randomly assigned +#' # to a stack of raster layers +#' test_fbpRaster <- rast( +#' system.file("extdata", "test_fbpRaster.tif", package = "cffdrs") +#' ) +#' input <- test_fbpRaster +#' # Stack doesn't hold the raster layer names, we have to assign +#' # them: +#' names(input) <- c( +#' "FuelType", "LAT", "LONG", "ELV", "FFMC", "BUI", "WS", "WD", "GS", +#' "Dj", "D0", "hr", "PC", "PDF", "GFL", "cc", "theta", "Accel", "Aspect", +#' "BUIEff", "CBH", "CFL", "ISI" +#' ) +#' # Primary outputs: +#' system.time(foo <- fbpRaster(input = input)) +#' # Using the "select" option: +#' system.time(foo <- fbpRaster(input = input, select = c("HFI", "TFC", "ROS"))) +#' # Secondary outputs: +#' system.time(foo <- fbpRaster(input = input, output = "S")) +#' # All outputs: +#' system.time(foo <- fbpRaster(input = input, output = "A")) +#' +#' ### Additional, longer running examples ### +#' # Keep only the required input layers, the other layers would be +#' # assigned with default values: +#' # keep only the required inputs: +#' dat0 <- input[[c( +#' "FuelType", "LAT", "LONG", "FFMC", "BUI", "WS", "GS", "Dj", "Aspect" +#' )]] +#' system.time(foo <- fbpRaster(input = dat0, output = "A")) +#' +#' @export fbpRaster +#' */ + +type InputType = { [key: string]: any }; + +interface FuelCrossType { + FUELTYPE0: string; + code: string; // corrected to string to match the FUEL code type in the input +} + +async function fbpRaster( + input: InputType, + output: string = "Primary", + select: string[] | null = null, + m: number | null = null, + cores: number = 1 +): Promise { + let output_orig = output; + output = output.toUpperCase(); + + // Ensure 'input' is not already attached + if (typeof input !== "undefined") { + console.warn("Attached dataset 'input' is being detached to use fbp() function."); + } + + // Determine optimal number of pixels to process at each iteration + if (m === null) { + m = input.cellCount > 500000 ? 3000 : 1000; + } + + // Setup correct output names + const allNames = [ + "CFB", "CFC", "FD", "HFI", "RAZ", "ROS", "SFC", "TFC", "BE", "SF", "ISI", + "FFMC", "FMC", "D0", "RSO", "CSI", "FROS", "BROS", "HROSt", "FROSt", + "BROSt", "FCFB", "BCFB", "FFI", "BFI", "FTFC", "BTFC", "TI", "FTI", + "BTI", "LB", "LBt", "WSV", "DH", "DB", "DF", "TROS", "TROSt", "TCFB", + "TFI", "TTFC", "TTI" + ]; + const primaryNames = allNames.slice(0, 8); + const secondaryNames = allNames.slice(8); + + // If outputs are specified, then check if they exist + if (select !== null) { + select = select.map(name => name.toUpperCase()); + select = [...new Set(select)]; + + const invalidSelection = (selection: string[], validNames: string[]) => { + return !selection.every(item => validNames.includes(item)); + }; + + if (["SECONDARY", "S"].includes(output)) { + if (invalidSelection(select, secondaryNames)) { + throw new Error("Selected variables are not in the outputs"); + } + } + + if (["PRIMARY", "P"].includes(output)) { + if (invalidSelection(select, primaryNames)) { + throw new Error("Selected variables are not in the outputs"); + } + } + + if (["ALL", "A"].includes(output)) { + if (invalidSelection(select, allNames)) { + throw new Error("Selected variables are not in the outputs"); + } + } + } + + // If caller specifies select outputs, then create a raster stack that contains only those outputs + if (select === null) { + if (["PRIMARY", "P"].includes(output)) { + select = primaryNames; + } else if (["SECONDARY", "S"].includes(output)) { + select = secondaryNames; + } else if (["ALL", "A"].includes(output)) { + select = allNames; + } else { + throw new Error(`Invalid output selected. ${output_orig} cannot be returned.`); + } + } + + input.names = input.names.map((name: string) => name.toUpperCase()); + + // Create or update the input with LAT and LON + let r: InputType; + if (!input.names.includes("LAT")) { + r = toDataTable(input); + const coords = transformCoordinates(r.x, r.y, input.crs, 4326); + r.LON = coords.x; + r.LAT = coords.y; + + if (Math.max(...r.LAT) > 90 || Math.min(...r.LAT) < -90) { + console.warn("Input projection is not in lat/long, consider re-projection or include LAT as input"); + } + } else { + r = toDataTable(input); + } + + // Unique IDs + r.ID = r.x.map((_: any, index: number) => index + 1); + + // Merge fuel codes with integer values + const fuelCross: FuelCrossType[] = [ + { FUELTYPE0: "C-1", code: "1" }, { FUELTYPE0: "C-2", code: "2" }, { FUELTYPE0: "C-3", code: "3" }, { FUELTYPE0: "C-4", code: "4" }, { FUELTYPE0: "C-5", code: "5" }, + { FUELTYPE0: "C-6", code: "6" }, { FUELTYPE0: "C-7", code: "7" }, { FUELTYPE0: "D-1", code: "8" }, { FUELTYPE0: "M-1", code: "9" }, { FUELTYPE0: "M-2", code: "10" }, + { FUELTYPE0: "M-3", code: "11" }, { FUELTYPE0: "M-4", code: "12" }, { FUELTYPE0: "S-1", code: "13" }, { FUELTYPE0: "S-2", code: "14" }, { FUELTYPE0: "S-3", code: "15" }, + { FUELTYPE0: "O-1a", code: "16" }, { FUELTYPE0: "O-1b", code: "17" }, { FUELTYPE0: "WA", code: "18" }, { FUELTYPE0: "NF", code: "19" } + ]; + + r.FUELTYPE = r.FUEL.map((fuel: string) => fuelCross.find(fc => fc.code === fuel)?.FUELTYPE0); + + // Calculate FBP through the fbp() function + const FBP = await fbp(r, output, m, cores); + + // Reassign character representation of Fire Type S/I/C to a numeric value 1/2/3 + if (!["SECONDARY", "S"].includes(output)) { + FBP.FD = FBP.FD.map((fd: string) => parseInt(fd.replace(/S/g, "1").replace(/I/g, "2").replace(/C/g, "3"))); + } + + if (select.includes("FD")) { + console.log("FD = 1,2,3 representing Surface (S), Intermittent (I), and Crown (C) fire"); + } + + let out: InputType[] = []; + select.forEach(name => { + out.push({ [name]: FBP[name] }); + }); + + return out; +} + +/*************************************************************************************** */ + +// Dummy function to simulate toDataTable and transformCoordinates in TypeScript +function toDataTable(input: any): any { + // This function should convert input to data table format + // Dummy implementation for demonstration + return input; +} + +function transformCoordinates(x: number[], y: number[], fromCRS: string, toCRS: number): { x: number[], y: number[] } { + // This function should transform coordinates from one CRS to another + // Dummy implementation for demonstration + return { x, y }; +} + +// Dummy function to simulate fbp in TypeScript +async function fbp(input: InputType, output: string, m: number, cores: number): Promise { + // This function should implement the fire behaviour prediction logic + // Dummy implementation for demonstration + return input; +} + +/*************************************************************************************** */ + +// Exporting the function for external usage +export { fbpRaster }; diff --git a/src/fwi/fine_fuel_moisture_code.ts b/src/fwi/fine_fuel_moisture_code.ts new file mode 100644 index 0000000..402e61a --- /dev/null +++ b/src/fwi/fine_fuel_moisture_code.ts @@ -0,0 +1,66 @@ +const FFMC_COEFFICIENT = 250.0 * 59.5 / 101.0; + +/** + * Fine Fuel Moisture Code Calculation + * + * @param ffmc_yda The Fine Fuel Moisture Code from the previous iteration + * @param temp Temperature (centigrade) + * @param rh Relative Humidity (%) + * @param ws Wind speed (km/h) + * @param prec Precipitation (mm) + * + * @return A single fine fuel moisture code value + */ +function fineFuelMoistureCode(ffmc_yda: number, temp: number, rh: number, ws: number, prec: number): number { + let wmo = FFMC_COEFFICIENT * (101 - ffmc_yda) / (59.5 + ffmc_yda); + + const ra = prec > 0.5 ? prec - 0.5 : prec; + + if (prec > 0.5) { + if (wmo > 150) { + wmo = wmo + 0.0015 * (wmo - 150) * (wmo - 150) * Math.sqrt(ra) + + 42.5 * ra * Math.exp(-100 / (251 - wmo)) * (1 - Math.exp(-6.93 / ra)); + } else { + wmo = wmo + 42.5 * ra * Math.exp(-100 / (251 - wmo)) * (1 - Math.exp(-6.93 / ra)); + } + } + + wmo = wmo > 250 ? 250 : wmo; + + const ed = 0.942 * Math.pow(rh, 0.679) + 11 * Math.exp((rh - 100) / 10) + + 0.18 * (21.1 - temp) * (1 - 1 / Math.exp(rh * 0.115)); + + const ew = 0.618 * Math.pow(rh, 0.753) + 10 * Math.exp((rh - 100) / 10) + + 0.18 * (21.1 - temp) * (1 - 1 / Math.exp(rh * 0.115)); + + let z = (wmo < ed && wmo < ew) ? + (0.424 * (1 - Math.pow((100 - rh) / 100, 1.7)) + 0.0694 * Math.sqrt(ws) * (1 - Math.pow((100 - rh) / 100, 8))) : + 0; + + let x = z * 0.581 * Math.exp(0.0365 * temp); + + let wm = (wmo < ed && wmo < ew) ? ew - (ew - wmo) / Math.pow(10, x) : wmo; + + z = (wmo > ed) ? + (0.424 * (1 - Math.pow(rh / 100, 1.7)) + 0.0694 * Math.sqrt(ws) * (1 - Math.pow(rh / 100, 8))) : + z; + + x = z * 0.581 * Math.exp(0.0365 * temp); + + wm = (wmo > ed) ? ed + (wmo - ed) / Math.pow(10, x) : wm; + + let ffmc1 = (59.5 * (250 - wm)) / (FFMC_COEFFICIENT + wm); + + ffmc1 = ffmc1 > 101 ? 101 : (ffmc1 < 0 ? 0 : ffmc1); + + return ffmc1; +} + +// Dummy function to simulate .ffmcCalc in TypeScript +function ffmcCalc(...args: [number, number, number, number, number]): number { + console.warn("Deprecated: fineFuelMoistureCode"); + return fineFuelMoistureCode(...args); +} + +// Exporting the functions for external usage +export { fineFuelMoistureCode, ffmcCalc }; diff --git a/src/fwi/fire_intensity.ts b/src/fwi/fire_intensity.ts new file mode 100644 index 0000000..fcac56b --- /dev/null +++ b/src/fwi/fire_intensity.ts @@ -0,0 +1,30 @@ +/** + * Fire Intensity Calculator + * + * @description Calculate the Predicted Fire Intensity + * + * All variables names are laid out in the same manner as Forestry Canada Fire + * Danger Group (FCFDG) (1992). Development and Structure of the Canadian Forest + * Fire Behavior Prediction System." Technical Report ST-X-3, Forestry Canada, + * Ottawa, Ontario. + * + * @param FC Fuel Consumption (kg/m^2) + * @param ROS Rate of Spread (m/min) + * + * @return FI: Fire Intensity (kW/m) + */ +function fireIntensity(FC: number, ROS: number): number { + // Eq. 69 (FCFDG 1992) Fire Intensity (kW/m) + const FI = 300 * FC * ROS; + return FI; + } + + // Dummy function to simulate .FIcalc in TypeScript + function FIcalc(...args: [number, number]): number { + console.warn("Deprecated: fireIntensity"); + return fireIntensity(...args); + } + + // Exporting the functions for external usage + export { fireIntensity, FIcalc }; + \ No newline at end of file diff --git a/src/fwi/fire_season.ts b/src/fwi/fire_season.ts new file mode 100644 index 0000000..28e3b2d --- /dev/null +++ b/src/fwi/fire_season.ts @@ -0,0 +1,291 @@ +/* #' Fire Season Start and End +#' +#' @description \code{\link{fire_season}} calculates the start and end fire +#' season dates for a given weather station. The current method used in the +#' function is based on three consecutive daily maximum temperature thresholds +#' (Wotton and Flannigan 1993, Lawson and Armitage 2008). This function process +#' input from a single weather station. +#' +#' An important aspect to consider when calculating Fire Weather Index (FWI) +#' System variables is a definition of the fire season start and end dates +#' (Lawson and Armitage 2008). If a user starts calculations on a fire season +#' too late in the year, the FWI System variables may take too long to reach +#' equilibrium, thus throwing off the resulting indices. This function presents +#' two method of calculating these start and end dates, adapted from Wotton and +#' Flannigan (1993), and Lawson and Armitage (2008). The approach taken in this +#' function starts the fire season after three days of maximum temperature +#' greater than 12 degrees Celsius. The end of the fire season is determined +#' after three consecutive days of maximum temperature less than 5 degrees +#' Celsius. The two temperature thresholds can be adjusted as parameters in +#' the function call. In regions where temperature thresholds will not end a +#' fire season, it is possible for the fire season to span multiple years, in +#' this case setting the multi.year parameter to TRUE will allow these +#' calculations to proceed. +#' +#' This fire season length definition can also feed in to the overwinter DC +#' calculations (\link{overwinter_drought_code}). View the cffdrs package help +#' files for an example of using the \code{fire_season}, +#' \link{overwinter_drought_code}, and \link{fwi} functions in conjunction. +#' +#' @param input A data.frame containing input variables of including the +#' date/time and daily maximum temperature. Variable names have to be the same +#' as in the following list, but they are case insensitive. The order in which +#' the input variables are entered is not important either. +#' +#' \tabular{lll}{ +#' \var{yr} \tab (required) \tab Year of the observations\cr +#' \var{mon} \tab (required) \tab Month of the observations\cr +#' \var{day} \tab (required) \tab Day of the observations\cr +#' \var{tmax} \tab (required) \tab Maximum Daily Temperature (degrees C)\cr +#' \var{snow_depth} +#' \tab (optional) \tab Is consistent snow data in the input?\cr +#' }. +#' +#' @param fs.start Temperature threshold (degrees C) to start the fire season +#' (default=12) +#' @param fs.end Temperature threshold (degrees C) to end the fire season +#' (default=5) +#' @param method Method of fire season calculation. Options are "wf93"" or +#' "la08" (default=WF93) +#' @param consistent.snow Is consistent snow data in the input? (default=FALSE) +#' @param multi.year Should the fire season span multiple years? +#' (default=FALSE) +#' @return \link{fire_season} returns a data frame of season and start and end +#' dates. Columns in data frame are described below. +#' +#' Primary FBP output includes the following 8 variables: +#' \item{yr }{Year of the fire season start/end date} +#' \item{mon }{Month of the fire season start/end date} +#' \item{day }{Day of the fire season start/end date} +#' \item{fsdatetype }{ +#' Fire season date type (values are either "start" or "end")} +#' \item{date}{Full date value} +#' +#' @author Alan Cantin, Xianli Wang, Mike Wotton, and Mike Flannigan +#' +#' @seealso \code{\link{fwi}, \link{overwinter_drought_code}} +#' +#' @references Wotton, B.M. and Flannigan, M.D. (1993). Length of the fire +#' season in a changing climate. Forestry Chronicle, 69, 187-192. +#' +#' \url{https://www.ualberta.ca/~flanniga/publications/1993_Wotton_Flannigan.pdf} +#' +#' Lawson, B.D. and O.B. Armitage. 2008. Weather guide for the Canadian Forest +#' Fire Danger Rating System. Nat. Resour. Can., Can. For. Serv., North. For. +#' Cent., Edmonton, AB \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/29152.pdf} +#' @keywords methods +#' @examples +#' +#' library(cffdrs) +#' # The standard test data: +#' data("test_wDC") +#' print(head(test_wDC)) +#' ## Sort the data: +#' input <- with(test_wDC, test_wDC[order(id, yr, mon, day), ]) +#' +#' # Using the default fire season start and end temperature +#' # thresholds: +#' a_fs <- fire_season(input[input$id == 1, ]) +#' +#' # Check the result: +#' a_fs +#' +#' # yr mon day fsdatetype +#' # 1 1999 5 4 start +#' # 2 1999 5 12 end +#' # 3 1999 5 18 start +#' # 4 1999 5 25 end +#' # 5 1999 5 30 start +#' # 6 1999 10 6 end +#' # 7 2000 6 27 start +#' # 8 2000 10 7 end +#' +#' # In the resulting data frame, the fire season starts +#' # and ends multiple times in the first year. It is up to the user for how +#' # to interpret this. +#' +#' # modified fire season start and end temperature thresholds +#' a_fs <- fire_season(input[input$id == 1, ], fs.start = 10, fs.end = 3) +#' a_fs +#' # yr mon day fsdatetype +#' # 1 1999 5 2 start +#' # 2 1999 10 20 end +#' # 3 2000 6 16 start +#' # 4 2000 10 7 end +#' # select another id value, specify method explicitly +#' b_fs <- fire_season(input[input$id == 2, ], method = "WF93") +#' # print the calculated fire_season +#' b_fs +#' # yr mon day fsdatetype +#' # 1 1980 4 21 start +#' # 2 1980 9 19 end +#' # 3 1980 10 6 start +#' # 4 1980 10 16 end +#' # 5 1981 5 21 start +#' # 6 1981 10 13 end +#' +#' @export fire_season + +fire_season <- function( + input, + fs.start = 12, + fs.end = 5, + method = "WF93", + consistent.snow = FALSE, + multi.year = FALSE) { + ############################################################################# + # Description: + # Calculation of fire season. The intent of this function is to allow + # for calculation of the start and end dates of a fire season. Knowledge + # of the start and end dates allows for more accurate calculation of Fire + # Weather Indices. Only 1 method is presented here, however future + # modifications will allow the inclusion of other methods. + # + # The current method contained within are contained in the following + # article: + # + # Wotton B.M. and Flannigan, M.D. 1993. Length of the fire season in a + # changing climate. Forestry Chronicle, 69: 187-192. + + # **Note that snow depth/cover is not currently included due to a lack of + # widespread snow coverage data, however this will be a future addition. + # + # Args: + # input: Single station weather data input stream - view fire_season.Rd + # documentation for full description. + # fs.start: Temperature threshold to start the fire season (deg celcius) + # fs.end: Temperature threshold to end the fire season (deg celcius) + # method: Fire Season calculation method. The default and only current + # method being used is "WF93". + # consistent.snow TRUE/FALSE value if consistent snow data is available + # if it is not, then we will default to WF93. + # + # + # Returns: + # seasonStartEnd: data.frame containing start and end dates of fire + # seasons for a particular weather station. + # + ############################################################################# */ + +type InputType = { + yr: number[]; + mon: number[]; + day: number[]; + tmax: number[]; + snow_depth?: number[]; + }; + + type SeasonStartEndType = { + yr: number; + mon: number; + day: number; + fsdatetype: string; + date?: Date; + }; + + function fireSeason( + input: InputType, + fsStart = 12, + fsEnd = 5, + method = "WF93", + consistentSnow = false, + multiYear = false + ): SeasonStartEndType[] { + const names = Object.keys(input).map(key => key.toLowerCase()); + const yr = input.yr; + const mon = input.mon; + const day = input.day; + const tmax = input.tmax; + const snowDepth = input.snow_depth || []; + + if (!yr) throw new Error("Year was not provided, year is required for this function."); + if (!mon) throw new Error("Month was not provided, month is required for this function."); + if (!day) throw new Error("Day was not provided, day is required for this function."); + if (!tmax) throw new Error("Maximum Daily Temperature (tmax) was not provided, daily tmax is required for this function."); + + method = method.toLowerCase(); + if (method !== "wf93" && method !== "la08") { + throw new Error(`Selected method "${method}" is unavailable, read documentation for available methods.`); + } + + const n0 = tmax.length; + let seasonActive = false; + let seasonStartEnd: SeasonStartEndType[] = []; + + if (method === "wf93") { + for (let k = 3; k < n0; k++) { + if (!seasonActive && tmax.slice(k - 3, k).every(temp => temp > fsStart)) { + seasonActive = true; + let theday = day[k]; + if (!multiYear && mon[k] === 1 && day[k] === 4) { + theday = day[k - 3]; + } + seasonStartEnd.push({ + yr: yr[k], + mon: mon[k], + day: theday, + fsdatetype: "start" + }); + } + + if (seasonActive && tmax.slice(k - 3, k).every(temp => temp < fsEnd)) { + seasonActive = false; + seasonStartEnd.push({ + yr: yr[k], + mon: mon[k], + day: day[k], + fsdatetype: "end" + }); + } + } + } else if (method === "la08") { + if (consistentSnow) { + if (!snowDepth.length) { + throw new Error(`Snow depth is required for the selected method "${method}", read documentation for appropriate use.`); + } + + for (let k = 3; k < n0; k++) { + if (!seasonActive && snowDepth.slice(k - 2, k + 1).every(depth => depth <= 0)) { + seasonActive = true; + let theday = day[k]; + if (!multiYear && mon[k] === 1 && day[k] === 4) { + theday = day[k - 3]; + } + seasonStartEnd.push({ + yr: yr[k], + mon: mon[k], + day: theday, + fsdatetype: "start" + }); + } + + if (seasonActive && (snowDepth[k] > 0 || (mon[k] === 12 && tmax.slice(k - 2, k + 1).every(temp => temp < fsEnd)))) { + seasonActive = false; + seasonStartEnd.push({ + yr: yr[k], + mon: mon[k], + day: day[k], + fsdatetype: "end" + }); + } + } + } else { + return fireSeason(input, fsStart, fsEnd, "WF93", consistentSnow, multiYear); + } + } + + seasonStartEnd = seasonStartEnd.filter((value, index, self) => + index === self.findIndex((t) => ( + t.yr === value.yr && t.mon === value.mon && t.day === value.day + )) + ); + + seasonStartEnd.forEach(record => { + record.date = new Date(`${record.yr}-${record.mon}-${record.day}`); + }); + + return seasonStartEnd; + } + + export { fireSeason }; + \ No newline at end of file diff --git a/src/fwi/fire_weather_index.ts b/src/fwi/fire_weather_index.ts new file mode 100644 index 0000000..8923449 --- /dev/null +++ b/src/fwi/fire_weather_index.ts @@ -0,0 +1,45 @@ +/** + * Fire Weather Index Calculation. + * + * @description All code is based on a C code library that was written by + * Canadian Forest Service Employees, which was originally based on the Fortran + * code listed in the reference below. All equations in this code refer to that + * document. + * + * Equations and FORTRAN program for the Canadian Forest Fire + * Weather Index System. 1985. Van Wagner, C.E.; Pickett, T.L. + * Canadian Forestry Service, Petawawa National Forestry + * Institute, Chalk River, Ontario. Forestry Technical Report 33. + * 18 p. + * + * Additional reference on FWI system + * + * Development and structure of the Canadian Forest Fire Weather + * Index System. 1987. Van Wagner, C.E. Canadian Forestry Service, + * Headquarters, Ottawa. Forestry Technical Report 35. 35 p. + * + * @param isi Initial Spread Index + * @param bui Buildup Index + * + * @return A single fwi value + */ + +function fireWeatherIndex(isi: number, bui: number): number { + // Eqs. 28b, 28a, 29 + const bb = bui > 80 + ? 0.1 * isi * (1000 / (25 + 108.64 / Math.exp(0.023 * bui))) + : 0.1 * isi * (0.626 * Math.pow(bui, 0.809) + 2); + + // Eqs. 30b, 30a + const fwi = bb <= 1 ? bb : Math.exp(2.72 * Math.pow(0.434 * Math.log(bb), 0.647)); + return fwi; + } + + // Dummy function to simulate .fwiCalc in TypeScript + function fwiCalc(...args: [number, number]): number { + console.warn("Deprecated: fireWeatherIndex"); + return fireWeatherIndex(...args); + } + + // Exporting the functions for external usage + export { fireWeatherIndex, fwiCalc }; diff --git a/src/fwi/flank_rate_of_spread.ts b/src/fwi/flank_rate_of_spread.ts new file mode 100644 index 0000000..5f8613d --- /dev/null +++ b/src/fwi/flank_rate_of_spread.ts @@ -0,0 +1,32 @@ +/** + * Flank Fire Rate of Spread Calculator + * + * @description Calculate the Flank Fire Spread Rate. + * + * All variables names are laid out in the same manner as Forestry Canada + * Fire Danger Group (FCFDG) (1992). Development and Structure of the + * Canadian Forest Fire Behavior Prediction System." Technical Report + * ST-X-3, Forestry Canada, Ottawa, Ontario. + * + * @param ROS Fire Rate of Spread (m/min) + * @param BROS Back Fire Rate of Spread (m/min) + * @param LB Length to breadth ratio + * + * @return FROS Flank Fire Spread Rate (m/min) value + */ + +function flankRateOfSpread(ROS: number, BROS: number, LB: number): number { + // Eq. 89 (FCFDG 1992) + const FROS = (ROS + BROS) / LB / 2; + return FROS; + } + + // Dummy function to simulate .FROScalc in TypeScript + function FROScalc(...args: [number, number, number]): number { + console.warn("Deprecated: flankRateOfSpread"); + return flankRateOfSpread(...args); + } + + // Exporting the functions for external usage + export { flankRateOfSpread, FROScalc }; + \ No newline at end of file diff --git a/src/fwi/foliar_moisture_content.ts b/src/fwi/foliar_moisture_content.ts new file mode 100644 index 0000000..7f99d68 --- /dev/null +++ b/src/fwi/foliar_moisture_content.ts @@ -0,0 +1,69 @@ +/** + * Foliar Moisture Content Calculator + * + * @description Calculate Foliar Moisture Content on a specified day. + * All variables names are laid out in the same manner as Forestry Canada + * Fire Danger Group (FCFDG) (1992). Development and Structure of the + * Canadian Forest Fire Behavior Prediction System." Technical Report + * ST-X-3, Forestry Canada, Ottawa, Ontario. + * + * @param LAT Latitude (decimal degrees) + * @param LONG Longitude (decimal degrees) + * @param ELV Elevation (metres) + * @param DJ Day of year (often referred to as Julian date) + * @param D0 Date of minimum foliar moisture content. If D0, date of min + * FMC, is not known then D0 = null. + * + * @return FMC: Foliar Moisture Content value + */ + +function foliarMoistureContent(LAT: number[], LONG: number[], ELV: number[], DJ: number[], D0: (number | null)[]): number[] { + const FMC: number[] = new Array(LAT.length).fill(-1); + let LATN: number[] = new Array(LAT.length).fill(0); + + // Calculate Normalized Latitude + // Eqs. 1 & 3 (FCFDG 1992) + LATN = LAT.map((lat, i) => { + if (D0[i] !== null && D0[i]! > 0) return LATN[i]; + return (ELV[i] <= 0) + ? 46 + 23.4 * Math.exp(-0.0360 * (150 - LONG[i])) + : 43 + 33.7 * Math.exp(-0.0351 * (150 - LONG[i])); + }); + + // Calculate Date of minimum foliar moisture content + // Eqs. 2 & 4 (FCFDG 1992) + D0 = D0.map((d0, i) => { + if (d0 !== null && d0 > 0) return d0; + return (ELV[i] <= 0) + ? 151 * (LAT[i] / LATN[i]) + : 142.1 * (LAT[i] / LATN[i]) + 0.0172 * ELV[i]; + }); + + // Round D0 to the nearest integer because it is a date + D0 = D0.map(d0 => Math.round(d0!)); + + // Number of days between day of year and date of min FMC + // Eq. 5 (FCFDG 1992) + const ND = DJ.map((dj, i) => Math.abs(dj - D0[i]!)); + + // Calculate final FMC + // Eqs. 6, 7, & 8 (FCFDG 1992) + ND.forEach((nd, i) => { + FMC[i] = (nd < 30) + ? 85 + 0.0189 * Math.pow(nd, 2) + : (nd >= 30 && nd < 50) + ? 32.9 + 3.17 * nd - 0.0288 * Math.pow(nd, 2) + : 120; + }); + + return FMC; + } + + // Dummy function to simulate .FMCcalc in TypeScript + function FMCcalc(...args: [number[], number[], number[], number[], (number | null)[]]): number[] { + console.warn("Deprecated: foliarMoistureContent"); + return foliarMoistureContent(...args); + } + + // Exporting the functions for external usage + export { foliarMoistureContent, FMCcalc }; diff --git a/src/fwi/fwi.ts b/src/fwi/fwi.ts index 5a19dcb..f20a021 100644 --- a/src/fwi/fwi.ts +++ b/src/fwi/fwi.ts @@ -1,38 +1,458 @@ -/** - * Fire Weather Index Calculation. All code - * is based on a C code library that was written by Canadian - * Forest Service Employees, which was originally based on - * the Fortran code listed in the reference below. All equations - * in this code refer to that document. - * - * Equations and FORTRAN program for the Canadian Forest Fire - * Weather Index System. 1985. Van Wagner, C.E.; Pickett, T.L. - * Canadian Forestry Service, Petawawa National Forestry - * Institute, Chalk River, Ontario. Forestry Technical Report 33. - * 18 p. - * - * Additional reference on FWI system - * - * Development and structure of the Canadian Forest Fire Weather - * Index System. 1987. Van Wagner, C.E. Canadian Forestry Service, - * Headquarters, Ottawa. Forestry Technical Report 35. 35 p. - * - * @param {number} isi - Initial Spread Index - * @param {number} bui - Buildup Index - * @returns {number} A single fwi value - */ -export function fwi(isi: number, bui: number) { - // FIX: works for now - const exp = Math.exp; - const log = Math.log; - //Eqs. 28b, 28a, 29 - const bb = - bui > 80 - ? 0.1 * isi * (1000 / (25 + 108.64 / exp(0.023 * bui))) - : 0.1 * isi * (0.626 * (bui ** 0.809) + 2); - //Eqs. 30b, 30a - const fwi1 = bb <= 1 ? bb : exp(2.72 * ((0.434 * log(bb)) ** 0.647)); - return fwi1; -} +/* #' Fire Weather Index System +#' +#' @description \code{fwi} is used to calculate the outputs of the Canadian +#' Forest Fire Weather Index (FWI) System for one day or one fire season based +#' on noon local standard time (LST) weather observations of temperature, +#' relative humidity, wind speed, and 24-hour rainfall, as well as the previous +#' day's fuel moisture conditions. This function could be used for either one +#' weather station or for multiple weather stations. +#' +#' The Canadian Forest Fire Weather Index (FWI) System is a major subsystem of +#' the Canadian Forest Fire Danger Rating System, which also includes Canadian +#' Forest Fire Behavior Prediction (FBP) System. The modern FWI System was +#' first issued in 1970 and is the result of work by numerous researchers from +#' across Canada. It evolved from field research which began in the 1930's and +#' regional fire hazard and fire danger tables developed from that early +#' research. +#' +#' The modern System (Van Wagner 1987) provides six output indices which +#' represent fuel moisture and potential fire behavior in a standard pine +#' forest fuel type. Inputs are a daily noon observation of fire weather, which +#' consists of screen-level air temperature and relative humidity, 10 meter +#' open wind speed and 24 accumulated precipitation. +#' +#' The first three outputs of the system (the Fire Fuel Moisture Code (ffmc), +#' the Duff Moisture Code (dmc), and the Drought Code (dc)) track moisture in +#' different layers of the fuel making up the forest floor. Their calculation +#' relies on the daily fire weather observation and also, importantly, the +#' moisture code value from the previous day as they are in essence bookkeeping +#' systems tracking the amount of moisture (water) in to and out of the layer. +#' It is therefore important that when calculating FWI System outputs over an +#' entire fire season, an uninterrupted daily weather stream is provided; one +#' day is the assumed time step in the models and thus missing data must be +#' filled in. +#' +#' The next three outputs of the System are relative (unitless) indicators of +#' aspects of fire behavior potential: spread rate (the Initial Spread Index, +#' isi), fuel consumption (the Build-up Index, bui) and fire intensity per unit +#' length of fire front (the Fire Weather Index, fwi). This final index, the +#' fwi, is the component of the System used to establish the daily fire danger +#' level for a region and communicated to the public. This final index can be +#' transformed to the Daily Severity Rating (dsr) to provide a more +#' reasonably-scaled estimate of fire control difficulty. +#' +#' Both the Duff Moisture Code (dmc) and Drought Code (dc) are influenced by +#' day length (see Van Wagner 1987). Day length adjustments for different +#' ranges in latitude can be used (as described in Lawson and Armitage 2008 +#' (\url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/29152.pdf})) and are included +#' in this R function; latitude must be positive in the northern hemisphere and +#' negative in the southern hemisphere. +#' +#' The default initial (i.e., "start-up") fuel moisture code values (FFMC=85, +#' DMC=6, DC=15) provide a reasonable set of conditions for most springtime +#' conditions in Canada, the Northern U.S., and Alaska. They are not suitable +#' for particularly dry winters and are presumably not appropriate for +#' different parts of the world. +#' +#' @param input A dataframe containing input variables of daily weather +#' observations taken at noon LST. Variable names have to be the same as in the +#' following list, but they are case insensitive. The order in which the input +#' variables are entered is not important. +#' +#' \tabular{lll}{ +#' \var{id} \tab (optional) +#' \tab Unique identifier of a weather\cr +#' \tab\tab station or spatial point (no restriction on\cr +#' \tab\tab data type); required when \code{batch=TRUE}\cr +#' \var{lat} \tab (recommended) \tab Latitude (decimal degree, default=55)\cr +#' \var{long} \tab (optional) \tab Longitude (decimal degree)\cr +#' \var{yr} \tab (optional) \tab Year of observation; +#' required when \code{batch=TRUE}\cr +#' \var{mon} \tab (recommended) \tab Month of the year +#' (integer 1-12, default=7)\cr +#' \var{day} \tab (optional) \tab Dayof the month (integer); required when +#' \code{batch=TRUE}\cr +#' \var{temp} \tab (required) \tab Temperature (centigrade)\cr +#' \var{rh} \tab (required) \tab Relative humidity (\%)\cr +#' \var{ws} \tab (required) \tab 10-m height wind speed (km/h)\cr +#' \var{prec} \tab (required) \tab 24-hour rainfall (mm)\cr } +#' +#' @param init A data.frame or vector contains either the initial values for +#' FFMC, DMC, and DC or the same variables that were calculated for the +#' previous day and will be used for the current day's calculation. The +#' function also accepts a vector if the initial or previous day FWI values is +#' for only one weather station (a warning message comes up if a single set of +#' initial values is used for multiple weather stations). Defaults are the +#' standard initial values for FFMC, DMC, and DC defined as the following: +#' \tabular{lll}{ +#' \bold{Variable} \tab \bold{Description} \tab \bold{Default} \cr +#' \var{ffmc} +#' \tab Previous day Fine Fuel Moisture Code (FFMC; unitless) \tab 85 \cr +#' \var{dmc} \tab Previous day Duff Moisture Code (DMC; unitless)\tab 6 \cr +#' \var{dc} \tab Previous Day Drought Code (DC; unitless) \tab 15\cr +#' \var{lat} \tab Latitude of the weather station (\emph{Optional}) \tab 55 \cr} +#' +#' @param batch Whether the computation is iterative or single step, default is +#' TRUE. When \code{batch=TRUE}, the function will calculate daily FWI System +#' outputs for one weather station over a period of time chronologically with +#' the initial conditions given (\code{init}) applied only to the first day of +#' calculation. If multiple weather stations are processed, an additional "id" +#' column is required in the input to label different stations, and the data +#' needs to be sorted by date/time and "id". If \code{batch=FALSE}, the +#' function calculates only one time step (1 day) base on either the initial +#' start values or the previous day's FWI System variables, which should also +#' be assigned to \code{init} argument. +#' +#' @param out The function offers two output options, \code{out="all"} will +#' produce a data frame that includes both the input and the FWI System +#' outputs; \code{out="fwi"} will generate a data frame with only the FWI +#' system components. +#' +#' @param lat.adjust The function offers options for whether day length +#' adjustments should be applied to the calculations. The default value is +#' "TRUE". +#' +#' @param uppercase Output in upper cases or lower cases would be decided by +#' this argument. Default is TRUE. +#' +#' @return \code{fwi} returns a dataframe which includes both the input and the +#' FWI System variables as described below: +#' \item{Input Variables }{Including temp, rh, ws, and prec with id, long, lat, +#' yr, mon, or day as optional.} +#' \item{ffmc }{Fine Fuel Moisture Code} +#' \item{dmc }{Duff Moisture Code} +#' \item{dc }{Drought Code} +#' \item{isi }{Initial Spread Index} +#' \item{bui }{Buildup Index} +#' \item{fwi }{Fire Weather Index} +#' \item{dsr }{Daily Severity Rating} +#' +#' @author Xianli Wang, Alan Cantin, Marc-André Parisien, Mike Wotton, Kerry +#' Anderson, and Mike Flannigan +#' +#' @seealso \code{\link{fbp}}, \code{\link{fwiRaster}}, \code{\link{gfmc}}, +#' \code{\link{hffmc}}, \code{\link{hffmcRaster}}, \code{\link{sdmc}}, +#' \code{\link{overwinter_drought_code}}, \code{\link{fire_season}} +#' +#' @references 1. Van Wagner, C.E. and T.L. Pickett. 1985. Equations and +#' FORTRAN program for the Canadian Forest Fire Weather Index System. Can. For. +#' Serv., Ottawa, Ont. For. Tech. Rep. 33. 18 p. +#' \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/19973.pdf} +#' +#' 2. Van Wagner, C.E. 1987. Development and structure of the Canadian forest +#' fire weather index system. Forest Technology Report 35. (Canadian Forestry +#' Service: Ottawa). \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/19927.pdf} +#' +#' 3. Lawson, B.D. and O.B. Armitage. 2008. Weather guide for the Canadian +#' Forest Fire Danger Rating System. Nat. Resour. Can., Can. For. Serv., North. +#' For. Cent., Edmonton, AB. +#' \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/29152.pdf} +#' +#' @keywords methods +#' +#' @examples +#' +#' library(cffdrs) +#' # The test data is a standard test +#' # dataset for FWI system (Van Wagner and Pickett 1985) +#' data("test_fwi") +#' # Show the data, which is already sorted by time: +#' head(test_fwi) +#' # long lat yr mon day temp rh ws prec +#' # -100 40 1985 4 13 17 42 25 0 +#' # -100 40 1985 4 14 20 21 25 2.4 +#' # -100 40 1985 4 15 8.5 40 17 0 +#' # -100 40 1985 4 16 6.5 25 6 0 +#' # -100 40 1985 4 17 13 34 24 0 +#' +#' ## (1) FWI System variables for a single weather station: +#' # Using the default initial values and batch argument, +#' # the function calculate FWI variables chronically: +#' fwi.out1 <- fwi(test_fwi) +#' # Using a different set of initial values: +#' fwi.out2 <- fwi( +#' test_fwi, +#' init = data.frame(ffmc = 80, dmc = 10, dc = 16, lat = 50) +#' ) +#' # This could also be done as the following: +#' fwi.out2 <- fwi(test_fwi, init = data.frame(80, 10, 6, 50)) +#' # Or: +#' fwi.out2 <- fwi(test_fwi, init = c(80, 10, 6, 50)) +#' # Latitude could be ignored, and the default value (55) will +#' # be used: +#' fwi.out2 <- fwi(test_fwi, init = data.frame(80, 10, 6)) +#' +#' ## (2) FWI for one or multiple stations in a single day: +#' # Change batch argument to FALSE, fwi calculates FWI +#' # components based on previous day's fwi outputs: +#' +#' fwi.out3 <- fwi(test_fwi, init = fwi.out1, batch = FALSE) +#' # Using a suite of initials, assuming variables from fwi.out1 +#' # are the initial values for different records. +#' init_suite <- fwi.out1[, c("FFMC", "DMC", "DC", "LAT")] +#' # Calculating FWI variables for one day but with multiple +#' # stations. Because the calculations is for one time step, +#' # batch=FALSE: +#' fwi.out4 <- fwi(test_fwi, init = init_suite, batch = FALSE) +#' +#' ## (3) FWI for multiple weather stations over a period of time: +#' # Assuming there are 4 weather stations in the test dataset, and they are +#' # ordered by day: +#' test_fwi$day <- rep(1:(nrow(test_fwi) / 4), each = 4) +#' test_fwi$id <- rep(1:4, length(unique(test_fwi$day))) +#' # Running the function with the same default initial inputs, will receive a +#' # warning message, but that is fine: +#' fwi(test_fwi) +#' +#' ## (4) Daylength adjustment: +#' # Change latitude values where the monthly daylength adjustments +#' # are different from the standard ones +#' test_fwi$lat <- 22 +#' # With daylength adjustment +#' fwi(test_fwi)[1:3, ] +#' # Without daylength adjustment +#' fwi(test_fwi, lat.adjust = FALSE)[1:3, ] +#' +#' @export fwi +#' +fwi <- function( + input, + init = data.frame(ffmc = 85, dmc = 6, dc = 15, lat = 55), + batch = TRUE, + out = "all", + lat.adjust = TRUE, + uppercase = TRUE) { + ############################################################################# + # Description: Canadian Forest Fire Weather Index Calculations. All code + # is based on a C code library that was written by Canadian + # Forest Service Employees, which was originally based on + # the Fortran code listed in the reference below. All equations + # in this code refer to that document, unless otherwise noted. + # + # Equations and FORTRAN program for the Canadian Forest Fire + # Weather Index System. 1985. Van Wagner, C.E.; Pickett, T.L. + # Canadian Forestry Service, Petawawa National Forestry + # Institute, Chalk River, Ontario. Forestry Technical Report 33. + # 18 p. + # + # Additional reference on FWI system + # + # Development and structure of the Canadian Forest Fire Weather + # Index System. 1987. Van Wagner, C.E. Canadian Forestry Service, + # Headquarters, Ottawa. Forestry Technical Report 35. 35 p. + # + # Args: input: View Documentation (fwi.Rd) for full description + # of input data frame + # init: Initializing moisture values + # ffmc: Fine Fuel Moisture Code (default 85) + # dmc: Duff Moisture Code (default 6) + # dc: Drought Code (default 15) + # lat: Latitude (decimal degrees, default 55) + # batch: Function can be run in a batch mode, where multiple + # weather stations or points can be calculated at once. + # (TRUE/FALSE, default TRUE) + # out: Display the calculated FWI values, with or without the + # inputs. (all/fwi, default all) + # lat.adjust: Option to adjust day length in the calculations + # (TRUE/FALSE, default TRUE) + # uppercase: Output names in upper or lower case - a commonly + # asked for feature, as dataset naming conventions vary + # considerably. (TRUE/FALSE, default TRUE) + # + # + # Returns: A data.frame of the calculated FWI values with or without + # the input data attached to it. + # + ############################################################################# */ -export default fwi; +interface InputType { + id?: number | string; + lat?: number; + long?: number; + yr?: number; + mon?: number; + day?: number; + temp: number; + rh: number; + ws: number; + prec: number; + [key: string]: any; + } + + interface InitType { + ffmc: number; + dmc: number; + dc: number; + lat?: number; + } + + interface FWIOutput { + ffmc: number; + dmc: number; + dc: number; + isi: number; + bui: number; + fwi: number; + dsr: number; + [key: string]: any; + } + +/*************************************************************************************** */ + + function fineFuelMoistureCode( + ffmc_yda: number, temp: number, rh: number, ws: number, prec: number + ): number { + return Math.random() * 100; + } + + function duffMoistureCode( + dmc_yda: number, temp: number, rh: number, prec: number, lat: number, mon: number, lat_adjust: boolean + ): number { + return Math.random() * 100; + } + + function droughtCode( + dc_yda: number, temp: number, rh: number, prec: number, lat: number, mon: number, lat_adjust: boolean + ): number { + return Math.random() * 100; + } + + function initialSpreadIndex(ffmc: number, ws: number, flag: boolean): number { + return Math.random() * 100; + } + + function buildupIndex(dmc: number, dc: number): number { + return Math.random() * 100; + } + + function fireWeatherIndex(isi: number, bui: number): number { + return Math.random() * 100; + } + +/*************************************************************************************** */ + + function fwi( + input: InputType[], + init: InitType = { ffmc: 85, dmc: 6, dc: 15, lat: 55 }, + batch: boolean = true, + out: string = "all", + latAdjust: boolean = true, + uppercase: boolean = true + ): FWIOutput[] { + input.forEach(entry => { + Object.keys(entry).forEach(key => { + entry[key.toLowerCase()] = entry[key]; + }); + }); + + if (Array.isArray(init)) { + init = { ffmc: init[0], dmc: init[1], dc: init[2], lat: init[3] || 55 }; + } + + let { ffmc: ffmc_yda, dmc: dmc_yda, dc: dc_yda, lat: defaultLat } = init; + + let lat = input.map(entry => entry.lat ?? defaultLat); + let long = input.map(entry => entry.long ?? -120); + let yr = input.map(entry => entry.yr ?? 5000); + let mon = input.map(entry => entry.mon ?? 7); + let day = input.map(entry => entry.day ?? -99); + + if (batch) { + if (input.some(entry => entry.id !== undefined)) { + input.sort((a, b) => (a.yr! - b.yr!) || (a.mon! - b.mon!) || (a.day! - b.day!) || (a.id! > b.id! ? 1 : -1)); + } + } + + let temp = input.map(entry => entry.temp); + let prec = input.map(entry => entry.prec); + let ws = input.map(entry => entry.ws); + let rh = input.map(entry => entry.rh); + + let n = batch ? (new Set(input.map(entry => entry.id))).size : 1; + + let ffmc: number[] = []; + let dmc: number[] = []; + let dc: number[] = []; + let isi: number[] = []; + let bui: number[] = []; + let fwi: number[] = []; + let dsr: number[] = []; + + let n0 = temp.length / n; + + for (let i = 0; i < n0; i++) { + let k = [...Array(n).keys()].map(j => j + (i * n)); + k.forEach(j => rh[j] = rh[j] >= 100 ? 99.9999 : rh[j]); + + let ffmc1 = k.map(j => fineFuelMoistureCode(ffmc_yda, temp[j], rh[j], ws[j], prec[j])); + let dmc1 = k.map(j => duffMoistureCode(dmc_yda, temp[j], rh[j], prec[j], lat[j]!, mon[j], latAdjust)); + let dc1 = k.map(j => droughtCode(dc_yda, temp[j], rh[j], prec[j], lat[j]!, mon[j], latAdjust)); + let isi1 = k.map((j, idx) => initialSpreadIndex(ffmc1[idx], ws[j], false)); + let bui1 = k.map((j, idx) => buildupIndex(dmc1[idx], dc1[idx])); + let fwi1 = k.map((j, idx) => fireWeatherIndex(isi1[idx], bui1[idx])); + let dsr1 = fwi1.map(value => 0.0272 * Math.pow(value, 1.77)); + + ffmc.push(...ffmc1); + dmc.push(...dmc1); + dc.push(...dc1); + isi.push(...isi1); + bui.push(...bui1); + fwi.push(...fwi1); + dsr.push(...dsr1); + + ffmc_yda = ffmc1[ffmc1.length - 1]; + dmc_yda = dmc1[dmc1.length - 1]; + dc_yda = dc1[dc1.length - 1]; + } + + let newFWI: FWIOutput[] = []; + if (out === "fwi") { + newFWI = ffmc.map((value, index) => ({ + ffmc: value, + dmc: dmc[index], + dc: dc[index], + isi: isi[index], + bui: bui[index], + fwi: fwi[index], + dsr: dsr[index] + })); + } else { + if (out === "all") { + newFWI = input.map((entry, index) => ({ + ...entry, + ffmc: ffmc[index], + dmc: dmc[index], + dc: dc[index], + isi: isi[index], + bui: bui[index], + fwi: fwi[index], + dsr: dsr[index] + })); + } + } + + if (uppercase) { + newFWI = newFWI.map(entry => { + let newEntry: any = {}; + Object.keys(entry).forEach(key => { + newEntry[key.toUpperCase()] = entry[key]; + }); + return newEntry; + }); + } + + return newFWI; + } + + // Example usage + const testData: InputType[] = [ + { id: 1, long: -100, lat: 40, yr: 1985, mon: 4, day: 13, temp: 17, rh: 42, ws: 25, prec: 0 }, + { id: 1, long: -100, lat: 40, yr: 1985, mon: 4, day: 14, temp: 20, rh: 21, ws: 25, prec: 2.4 }, + { id: 1, long: -100, lat: 40, yr: 1985, mon: 4, day: 15, temp: 8.5, rh: 40, ws: 17, prec: 0 }, + { id: 1, long: -100, lat: 40, yr: 1985, mon: 4, day: 16, temp: 6.5, rh: 25, ws: 6, prec: 0 }, + { id: 1, long: -100, lat: 40, yr: 1985, mon: 4, day: 17, temp: 13, rh: 34, ws: 24, prec: 0 }, + // More data entries... + ]; + + console.log(fwi(testData, { ffmc: 85, dmc: 6, dc: 15, lat: 55 }, true, "all", true, true)); + \ No newline at end of file diff --git a/src/fwi/fwiRaster.ts b/src/fwi/fwiRaster.ts new file mode 100644 index 0000000..ddf2f20 --- /dev/null +++ b/src/fwi/fwiRaster.ts @@ -0,0 +1,279 @@ +/* #' Raster-based Fire Weather Index System +#' +#' @description \code{fwiRaster} is used to calculate the outputs of the +#' Canadian Forest Fire Weather Index (FWI) System for one day based on noon +#' local standard time (LST) weather observations of temperature, relative +#' humidity, wind speed, and 24-hour rainfall, as well as the previous day's +#' fuel moisture conditions. This function takes rasterized input and generates +#' raster maps as outputs. +#' +#' The Canadian Forest Fire Weather Index (FWI) System is a major subsystem of +#' the Canadian Forest Fire Danger Rating System, which also includes Canadian +#' Forest Fire Behavior Prediction (FBP) System. The modern FWI System was +#' first issued in 1970 and is the result of work by numerous researchers from +#' across Canada. It evolved from field research which began in the 1930's and +#' regional fire hazard and fire danger tables developed from that early +#' research. +#' +#' The modern System (Van Wagner 1987) provides six output indices which +#' represent fuel moisture and potential fire behavior in a standard pine +#' forest fuel type. Inputs are a daily noon observation of fire weather, which +#' consists of screen-level air temperature and relative humidity, 10 meter +#' open wind speed and 24 accumulated precipitation. +#' +#' The first three outputs of the system (the Fire Fuel Moisture Code, the Duff +#' Moisture Code, and the Drought Code) track moisture in different layers of +#' the fuel making up the forest floor. Their calculation relies on the daily +#' fire weather observation and also, importantly, the code value from the +#' previous day as they are in essence bookkeeping systems tracking the amount +#' of moisture (water) in to and out of the layer. It is therefore important +#' that when calculating FWI System outputs over an entire fire season, an +#' uninterrupted daily weather stream is provided; one day is the assumed time +#' step in the models and thus missing data must be filled in. +#' +#' The next three outputs of the System are relative (unitless) indicators of +#' aspects of fire behavior potential: spread rate (the Initial Spread Index), +#' fuel consumption (the Build-up Index) and fire intensity per unit length of +#' fire front (the Fire Weather Index). This final index, the fwi, is the +#' component of the System used to establish the daily fire danger level for a +#' region and communicated to the public. This final index can be transformed +#' to the Daily Severity Rating (dsr) to provide a more reasonably-scaled +#' estimate of fire control difficulty. +#' +#' Both the Duff Moisture Code (dmc) and Drought Code (dc) are influenced by +#' day length (see Van Wagner, 1987). Day length adjustments for different +#' ranges in latitude can be used (as described in Lawson and Armitage 2008 +#' (\url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/29152.pdf})) and are included +#' in this R function; latitude must be positive in the northern hemisphere and +#' negative in the southern hemisphere. +#' +#' The default initial (i.e., "start-up") fuel moisture code values (FFMC=85, +#' DMC=6, DC=15) provide a reasonable set of conditions for most springtime +#' conditions in Canada, the Northern U.S., and Alaska. They are not suitable +#' for particularly dry winters and are presumably not appropriate for +#' different parts of the world. +#' +#' @param input A stack or brick containing rasterized daily weather +#' observations taken at noon LST. Variable names have to be the same as in the +#' following list, but they are case insensitive. The order in which the inputs +#' are entered is not important. +#' +#' \tabular{lll}{ +#' \var{lat} \tab (recommended) \tab Latitude (decimal degree, +#' __default=55__)\cr +#' \var{temp} \tab (required) \tab Temperature (centigrade)\cr +#' \var{rh} \tab (required) \tab Relative humidity (\%)\cr +#' \var{ws} \tab +#' (required) \tab 10-m height wind speed (km/h)\cr +#' \var{prec} \tab (required) +#' \tab 24-hour rainfall (mm)\cr } +#' @param init A vector that contains the initial values for FFMC, DMC, and DC +#' or a stack that contains raster maps of the three moisture codes calculated +#' for the previous day, which will be used for the current day's calculation. +#' Defaults are the standard initial values for FFMC, DMC, and DC defined as +#' the following: +#' +#' \tabular{lll}{ +#' \bold{Variable} \tab \bold{Description} \tab \bold{Default} \cr +#' \var{ffmc} +#' \tab Previous day Fine Fuel Moisture Code (FFMC; unitless) \tab 85 \cr +#' \var{dmc} \tab Previous day Duff Moisture Code (DMC; unitless)\tab 6 \cr +#' \var{dc} \tab Previous Day Drought Code (DC; unitless) \tab 15\cr +#' \var{lat} \tab Latitude of the weather station (\emph{Optional})\tab 55 \cr} +#' +#' @param mon Month of the year (integer 1~12, default=7). Month is used in +#' latitude adjustment (\code{lat.adjust}), it is therefore recommended when +#' \code{lat.adjust=TRUE} was chosen. +#' @param out The function offers two output options, \code{out="all"} will +#' produce a raster stack include both the input and the FWI System outputs; +#' \code{out="fwi"} will generate a stack with only the FWI system components. +#' @param lat.adjust The function offers options for whether latitude +#' adjustments to day lengths should be applied to the calculations. The +#' default value is "TRUE". +#' @param uppercase Output in upper cases or lower cases would be decided by +#' this argument. Default is TRUE. +#' @return By default, \code{fwi} returns a raster stack which includes both +#' the input and the FWI System variables, as describe below: \item{Inputs +#' }{Including \code{temp}, \code{rh}, \code{ws}, and \code{prec} with +#' \code{lat} as optional.} \item{ffmc }{Fine Fuel Moisture Code} \item{dmc +#' }{Duff Moisture Code} \item{dc }{Drought Code} \item{isi }{Initial Spread +#' Index} \item{bui }{Buildup Index} \item{fwi }{Fire Weather Index} \item{dsr +#' }{Daily Severity Rating} +#' +#' @author Xianli Wang, Alan Cantin, Marc-André Parisien, Mike Wotton, Kerry +#' Anderson, and Mike Flannigan +#' +#' @seealso \code{\link{fbp}}, \code{\link{fbpRaster}}, \code{\link{fwi}}, +#' \code{\link{hffmc}}, \code{\link{hffmcRaster}} +#' +#' @references 1. Van Wagner, C.E. and T.L. Pickett. 1985. Equations and +#' FORTRAN program for the Canadian Forest Fire Weather Index System. Can. For. +#' Serv., Ottawa, Ont. For. Tech. Rep. 33. 18 p. +#' \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/19973.pdf} +#' +#' 2. Van Wagner, C.E. 1987. Development and structure of the Canadian forest +#' fire weather index system. Forest Technology Report 35. (Canadian Forestry +#' Service: Ottawa). \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/19927.pdf} +#' +#' 3. Lawson, B.D. and O.B. Armitage. 2008. Weather guide for the Canadian +#' Forest Fire Danger Rating System. Nat. Resour. Can., Can. For. Serv., North. +#' For. Cent., Edmonton, AB. +#' \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/29152.pdf} */ + +interface RasterInput { + [key: string]: number[] | undefined; + lat?: number[]; + temp: number[]; + rh: number[]; + ws: number[]; + prec: number[]; + } + + interface InitValues { + ffmc: number; + dmc: number; + dc: number; + lat?: number; + } + + interface FWIOutput { + ffmc: number[]; + dmc: number[]; + dc: number[]; + isi: number[]; + bui: number[]; + fwi: number[]; + dsr: number[]; + } + + function fineFuelMoistureCode( + ffmc_yda: number, temp: number, rh: number, ws: number, prec: number + ): number { + return Math.random() * 100; // Placeholder + } + + function duffMoistureCode( + dmc_yda: number, temp: number, rh: number, prec: number, lat: number, mon: number, latAdjust: boolean + ): number { + return Math.random() * 100; // Placeholder + } + + function droughtCode( + dc_yda: number, temp: number, rh: number, prec: number, lat: number, mon: number, latAdjust: boolean + ): number { + return Math.random() * 100; // Placeholder + } + + function initialSpreadIndex(ffmc: number, ws: number, fbpMod: boolean): number { + return Math.random() * 100; // Placeholder + } + + function buildupIndex(dmc: number, dc: number): number { + return Math.random() * 100; // Placeholder + } + + function fireWeatherIndex(isi: number, bui: number): number { + return Math.random() * 100; // Placeholder + } + + function fwiRaster( + input: RasterInput, + init: InitValues = { ffmc: 85, dmc: 6, dc: 15, lat: 55 }, + mon: number = 7, + out: string = "all", + latAdjust: boolean = true, + uppercase: boolean = true + ): FWIOutput | (FWIOutput & RasterInput) { + const requiredVars = ["temp", "rh", "ws", "prec"]; + for (const varName of requiredVars) { + if (!input[varName]) { + throw new Error(`${varName} is missing!`); + } + } + + if (input.prec.some((val: number) => val < 0)) { + throw new Error("precipitation (prec) cannot be negative!"); + } + if (input.ws.some((val: number) => val < 0)) { + throw new Error("wind speed (ws) cannot be negative!"); + } + if (input.rh.some((val: number) => val < 0)) { + throw new Error("relative humidity (rh) cannot be negative!"); + } + + const lat = input.lat ?? Array(input.temp.length).fill(init.lat ?? 55); + const ffmc_yda = Array(input.temp.length).fill(init.ffmc); + const dmc_yda = Array(input.temp.length).fill(init.dmc); + const dc_yda = Array(input.temp.length).fill(init.dc); + + input.rh = input.rh.map((val: number) => (val >= 100 ? 99.9999 : val)); + + const ffmc = input.temp.map((_, i) => + fineFuelMoistureCode(ffmc_yda[i], input.temp[i], input.rh[i], input.ws[i], input.prec[i]) + ); + + const dmc = input.temp.map((_, i) => + duffMoistureCode(dmc_yda[i], input.temp[i], input.rh[i], input.prec[i], lat[i], mon, latAdjust) + ); + + const dc = input.temp.map((_, i) => + droughtCode(dc_yda[i], input.temp[i], input.rh[i], input.prec[i], lat[i], mon, latAdjust) + ); + + const isi = input.temp.map((_, i) => + initialSpreadIndex(ffmc[i], input.ws[i], false) + ); + + const bui = input.temp.map((_, i) => + buildupIndex(dmc[i], dc[i]) + ); + + const fwi = input.temp.map((_, i) => + fireWeatherIndex(isi[i], bui[i]) + ); + + const dsr = fwi.map(val => 0.0272 * Math.pow(val, 1.77)); + + const newFWI: FWIOutput = { + ffmc, + dmc, + dc, + isi, + bui, + fwi, + dsr + }; + + if (uppercase) { + return { + ffmc: newFWI.ffmc, + dmc: newFWI.dmc, + dc: newFWI.dc, + isi: newFWI.isi, + bui: newFWI.bui, + fwi: newFWI.fwi, + dsr: newFWI.dsr + }; + } + + if (out === "fwi") { + return newFWI; + } + + if (out === "all") { + return { ...input, ...newFWI }; + } + + return newFWI; + } + + // Example usage + /* const testInput: RasterInput = { + temp: [17, 20, 8.5, 6.5, 13], + rh: [42, 21, 40, 25, 34], + ws: [25, 25, 17, 6, 24], + prec: [0, 2.4, 0, 0, 0] + }; */ + + //console.log(fwiRaster(testInput, { ffmc: 85, dmc: 6, dc: 15, lat: 55 }, 7, "all", true, true)); + \ No newline at end of file diff --git a/src/fwi/gfmc.ts b/src/fwi/gfmc.ts new file mode 100644 index 0000000..4c94217 --- /dev/null +++ b/src/fwi/gfmc.ts @@ -0,0 +1,277 @@ +/* #' Grass Fuel Moisture Code +#' +#' @description \code{gfmc} calculates both the moisture content of the surface +#' of a fully cured matted grass layer and also an equivalent Grass Fuel +#' Moisture Code (gfmc) (Wotton, 2009) to create a parallel with the hourly ffmc +#' (see the \code{\link{fwi}} and \code{\link{hffmc}}functions). The calculation +#' is based on hourly (or sub-hourly) weather observations of temperature, +#' relative humidity, wind speed, rainfall, and solar radiation. The user must +#' also estimate an initial value of the gfmc for the layer. This function +#' could be used for either one weather station or multiple weather stations. +#' +#' The Canadian Forest Fire Danger Rating System (CFFDRS) is used throughout +#' Canada, and in a number of countries throughout the world, for estimating +#' fire potential in wildland fuels. This new Grass Fuel Moisture Code (GFMC) +#' is an addition (Wotton 2009) to the CFFDRS and retains the structure of that +#' System's hourly Fine Fuel Moisture Code (HFFMC) (Van Wagner 1977). It tracks +#' moisture content in the top 5 cm of a fully-cured and fully-matted layer of +#' grass and thus is representative of typical after winter conditions in areas +#' that receive snowfall. This new moisture calculation method outputs both +#' the actual moisture content of the layer and also the transformed moisture +#' Code value using the FFMC's FF-scale. In the CFFDRS the moisture codes are +#' in fact relatively simple transformations of actual moisture content such +#' that decreasing moisture content (increasing dryness) is indicated by an +#' increasing Code value. This moisture calculation uses the same input weather +#' observations as the hourly FFMC, but also requires an estimate of solar +#' radiation incident on the fuel. +#' +#' @param input A dataframe containing input variables of daily noon weather +#' observations. Variable names have to be the same as in the following list, +#' but they are case insensitive. The order in which the input variables are +#' entered is not important. +#' +#' \tabular{lll}{ +#' \var{id} \tab (optional) \tab Batch Identification\cr +#' \var{temp} \tab (required) \tab Temperature (centigrade)\cr +#' \var{rh} \tab (required) \tab Relative humidity (\%)\cr +#' \var{ws} \tab (required) \tab 10-m height wind speed (km/h)\cr +#' \var{prec} \tab (required) \tab 1-hour rainfall (mm)\cr +#' \var{isol} \tab (required) \tab Solar radiation (kW/m^2)\cr +#' \var{mon} \tab (recommended) \tab Month of the year (integer' 1-12)\cr +#' \var{day} \tab (optional) \tab Day of the month (integer)\cr } +#' @param GFMCold Previous value of GFMC (i.e. value calculated at the previous +#' time step)[default is 85 (which corresponds to a moisture content of about +#' 16\%)]. On the first calculation this is the estimate of the GFMC value at +#' the start of the time step. The \code{GFMCold} argument can accept a single +#' initial value for multiple weather stations, and also accept a vector of +#' initial values for multiple weather stations. NOTE: this input represents +#' the CODE value, not a direct moisture content value. The CODE values in the +#' Canadian FWI System increase within decreasing moisture content. To roughly +#' convert a moisture content value to a CODE value on the FF-scale (used in +#' the FWI Systems FFMC) use \code{GFMCold} =101-gmc (where gmc is moisture +#' content in \%) +#' +#' @param time.step Time step (hour) [default 1 hour] +#' @param roFL The nominal fuel load of the fine fuel layer, default is 0.3 +#' kg/m^2 +#' @param batch Whether the computation is iterative or single step, default is +#' TRUE. When \code{batch=TRUE}, the function will calculate hourly or +#' sub-hourly GFMC for one weather station over a period of time iteratively. +#' If multiple weather stations are processed, an additional "id" column is +#' required in the input to label different stations, and the data needs to be +#' sorted by time sequence and "id". If \code{batch=FALSE}, the function +#' calculates only one time step (1 hour) base on either the previous hourly +#' GFMC or the initial start value. +#' @param out Output format, default is "GFMCandMC", which contains both GFMC +#' and moisture content (MC) in a data.frame format. Other choices include: +#' "GFMC", "MC", and "ALL", which include both the input and GFMC and MC. +#' @return \code{gfmc} returns GFMC and moisture content (MC) values +#' collectively (default) or separately. +#' @author Xianli Wang, Mike Wotton, Alan Cantin, and Mike Flannigan +#' @seealso \code{\link{fwi}}, \code{\link{hffmc}} +#' @references Wotton, B.M. 2009. A grass moisture model for the Canadian +#' Forest Fire Danger Rating System. In: Proceedings 8th Fire and Forest +#' Meteorology Symposium, Kalispell, MT Oct 13-15, 2009. Paper 3-2. +#' \url{https://ams.confex.com/ams/pdfpapers/155930.pdf} +#' +#' Van Wagner, C.E. 1977. A method of computing fine fuel moisture content +#' throughout the diurnal cycle. Environment Canada, Canadian Forestry Service, +#' Petawawa Forest Experiment Station, Chalk River, Ontario. Information Report +#' PS-X-69. \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/25591.pdf} +#' @keywords methods +#' @importFrom data.table data.table +#' @export gfmc +#' @examples +#' +#' library(cffdrs) +#' # load the test data +#' data("test_gfmc") +#' # show the data format: +#' head(test_gfmc) +#' # yr mon day hr temp rh ws prec isol +#' # 1 2006 5 17 10 15.8 54.6 5.0 0 0.340 +#' # 2 2006 5 17 11 16.3 52.9 5.0 0 0.380 +#' # 3 2006 5 17 12 18.8 45.1 5.0 0 0.626 +#' # 4 2006 5 17 13 20.4 40.8 9.5 0 0.656 +#' # 5 2006 5 17 14 20.1 41.7 8.7 0 0.657 +#' # 6 2006 5 17 15 18.6 45.8 13.5 0 0.629 +#' # (1) gfmc default: +#' # Re-order the data by year, month, day, and hour: +#' dat <- test_gfmc[with(test_gfmc, order(yr, mon, day, hr)), ] +#' # Because the test data has 24 hours input variables +#' # it is possible to calculate the hourly GFMC continuously +#' # through multiple days(with the default initial GFMCold=85): +#' dat$gfmc_default <- gfmc(dat,out="GFMC") +#' # two variables will be added to the input, GFMC and MC +#' head(dat) +#' # (2) For multiple weather stations: +#' # One time step (1 hour) with default initial value: +#' foo <- gfmc(dat, batch = FALSE) +#' # Chronological hourly GFMC with only one initial +#' # value (GFMCold=85), but multiple weather stations. +#' # Note: data is ordered by date/time and the station id. Subset +#' # the data by keeping only the first 10 hours of observations +#' # each day: +#' dat1 <- subset(dat, hr %in% c(0:9)) +#' # assuming observations were from the same day but with +#' # 9 different weather stations: +#' dat1$day <- NULL +#' dat1 <- dat1[with(dat1, order(yr, mon, hr)), ] +#' dat1$id <- rep(1:8, nrow(dat1) / 8) +#' # check the data: +#' head(dat1) +#' # Calculate GFMC for multiple stations: +#' dat1$gfmc01 <- gfmc(dat1, batch = TRUE) +#' # We can provide multiple initial GFMC (GFMCold) as a vector: +#' dat1$gfmc02 <- gfmc( +#' dat1, +#' GFMCold = sample(70:100, 8, replace = TRUE), +#' batch = TRUE +#' ) +#' # (3)output argument +#' ## include all inputs and outputs: +#' dat0 <- dat[with(dat, order(yr, mon, day, hr)), ] +#' foo <- gfmc(dat, out = "ALL") +#' ## subhourly time step: +#' gfmc(dat0, time.step = 1.5) */ + +interface GFMInput { + id?: number[]; + temp: number[]; + rh: number[]; + ws: number[]; + prec: number[]; + isol: number[]; + mon?: number[]; + day?: number[]; + } + + interface GFMOutput { + GFMC: number[]; + MC: number[]; + } + + function grassFuelMoisture( + temp: number, + rh: number, + ws: number, + prec: number, + isol: number, + GFMCold: number, + timeStep: number, + roFL: number + ): number { + // Placeholder for the actual calculation logic + return Math.random() * 100; // Placeholder + } + + function grassFuelMoistureCode(MC: number): number { + // Placeholder for the actual calculation logic + return 101 - MC; // Example conversion logic + } + + function gfmc( + input: GFMInput, + GFMCold: number | number[] = 85, + batch: boolean = true, + timeStep: number = 1, + roFL: number = 0.3, + out: string = "GFMCandMC" + ): GFMOutput | (GFMInput & GFMOutput) | { GFMC: number[] } | { MC: number[] } { + const requiredCols = [ + { full: "temperature", short: "temp" }, + { full: "precipitation", short: "prec" }, + { full: "wind speed", short: "ws" }, + { full: "relative humidity", short: "rh" }, + { full: "insolation", short: "isol" } + ]; + + for (const col of requiredCols) { + if (!(col.short in input)) { + throw new Error(`${col.full} is missing!`); + } + } + + let n: number; + if (batch) { + if (input.id) { + n = new Set(input.id).size; + if (n !== input.id.length) { + throw new Error( + "Multiple stations have to start and end at the same dates/time, and input data must be sorted by date/time and id" + ); + } + } else { + n = 1; + } + } else { + n = input.temp.length; + } + + if (input.temp.length % n !== 0) { + console.warn("Input data do not match with the number of weather stations"); + } + + if (typeof GFMCold === "number") { + GFMCold = Array(n).fill(GFMCold); + } + + if (GFMCold.length !== n) { + throw new Error("Number of GFMCold doesn't match number of weather stations"); + } + + const validOutTypes = ["GFMCandMC", "MC", "GFMC", "ALL"]; + if (!validOutTypes.includes(out)) { + throw new Error(`'${out}' is an invalid 'out' type.`); + } + + const n0 = Math.floor(input.temp.length / n); + let GFMC_out: number[] = []; + let MC_out: number[] = []; + + for (let i = 0; i < n0; i++) { + const k = Array.from({ length: n }, (_, index) => n * i + index); + + const MCStep = k.map(index => + grassFuelMoisture( + input.temp[index], + input.rh[index], + input.ws[index], + input.prec[index], + input.isol[index], + GFMCold[index], + timeStep, + roFL + ) + ); + + const GFMCStep = MCStep.map(moistureContent => grassFuelMoistureCode(moistureContent)); + + GFMCold = GFMCStep; + GFMC_out = GFMC_out.concat(GFMCStep); + MC_out = MC_out.concat(MCStep); + } + + switch (out) { + case "GFMC": + return { GFMC: GFMC_out }; + case "MC": + return { MC: MC_out }; + case "ALL": + return { ...input, GFMC: GFMC_out, MC: MC_out }; + default: + return { GFMC: GFMC_out, MC: MC_out }; + } + } + + // Example usage + const testInput: GFMInput = { + temp: [15.8, 16.3, 18.8, 20.4, 20.1, 18.6], + rh: [54.6, 52.9, 45.1, 40.8, 41.7, 45.8], + ws: [5.0, 5.0, 5.0, 9.5, 8.7, 13.5], + prec: [0, 0, 0, 0, 0, 0], + isol: [0.340, 0.380, 0.626, 0.656, 0.657, 0.629] + }; + + console.log(gfmc(testInput, 85, true, 1, 0.3, "GFMCandMC")); + \ No newline at end of file diff --git a/src/fwi/gfmcRaster.ts b/src/fwi/gfmcRaster.ts new file mode 100644 index 0000000..70fcd15 --- /dev/null +++ b/src/fwi/gfmcRaster.ts @@ -0,0 +1,145 @@ +/* #' @title Grass Fuel Moisture Raster Calculation +#' +#' @description Calculation of the Grass Fuel Moisture Code. This calculates the +#' moisture content of both the surface of a fully cured matted grass layer and +#' also an equivalent Grass Fuel Moisture Code. All equations come from Wotton +#' (2009) as cited below unless otherwise specified. +#' +#' @references +#' Wotton, B.M. 2009. A grass moisture model for the Canadian +#' Forest Fire Danger Rating System. In: Proceedings 8th Fire and +#' Forest Meteorology Symposium, Kalispell, MT Oct 13-15, 2009. +#' Paper 3-2. \url{https://ams.confex.com/ams/pdfpapers/155930.pdf} +#' +#' @param input [SpatRast stack] +#' \tabular{lll}{ +#' \var{temp} \tab (required) \tab Temperature (centigrade)\cr +#' \var{rh} \tab (required) \tab Relative humidity (\%)\cr +#' \var{ws} \tab (required) \tab 10-m height wind speed (km/h)\cr +#' \var{prec} \tab (required) \tab 1-hour rainfall (mm)\cr +#' \var{isol} \tab (required) \tab Solar radiation (kW/m^2)\cr } */ + +import * as terra from "terra"; // Assuming terra library is available in JS/TS +import { dataTable } from "data.table"; // Assuming a suitable data table library + +interface GFMInput { + temp: terra.Raster; + rh: terra.Raster; + ws: terra.Raster; + prec: terra.Raster; + isol: terra.Raster; +} + +type OutputType = "GFMCandMC" | "MC" | "GFMC" | "ALL"; + +function grassFuelMoisture( + temp: number, + rh: number, + ws: number, + prec: number, + isol: number, + GFMCold: number, + timeStep: number, + roFL: number +): number { + // Placeholder for the actual calculation logic + return Math.random() * 100; // Placeholder +} + +function grassFuelMoistureCode(MC: number): number { + // Placeholder for the actual calculation logic + return 101 - MC; // Example conversion logic +} + +function gfmcRaster( + input: GFMInput, + GFMCold: number | terra.Raster = 85, + timeStep: number = 1, + roFL: number | terra.Raster = 0.3, + out: OutputType = "GFMCandMC" +): terra.Raster | terra.Raster[] { + // Convert input to SpatRaster if necessary + if (!(input.temp instanceof terra.Raster)) { + input = { + temp: terra.rast(input.temp), + rh: terra.rast(input.rh), + ws: terra.rast(input.ws), + prec: terra.rast(input.prec), + isol: terra.rast(input.isol), + }; + } + const inputNames = Object.keys(input); + const requiredCols = { + full: ["temperature", "precipitation", "wind speed", "relative humidity", "insolation"], + short: ["temp", "prec", "ws", "rh", "isol"], + }; + + const missingCols = requiredCols.short.filter(col => !inputNames.includes(col)); + if (missingCols.length > 0) { + throw new Error(`${missingCols.join(", ")} is missing!`); + } + + // Convert GFMCold to SpatRaster if necessary + if (typeof GFMCold === "number") { + GFMCold = terra.setValues(input.temp, GFMCold); + } + if (GFMCold instanceof terra.Raster) { + terra.names(GFMCold, "GFMCold"); + } + + // Convert roFL to SpatRaster if necessary + if (typeof roFL === "number") { + roFL = terra.setValues(input.temp, roFL); + } + if (roFL instanceof terra.Raster) { + terra.names(roFL, "roFL"); + } + + // Ensure valid output type + const validOutTypes = ["GFMCandMC", "MC", "GFMC", "ALL"]; + if (!validOutTypes.includes(out)) { + throw new Error(`'${out}' is an invalid 'out' type.`); + } + + // Calculate MC and GFMC rasters + const mcRaster = terra.lapp( + [ + input.temp, + input.rh, + input.ws, + input.prec, + input.isol, + GFMCold, + roFL + ], + (temp, rh, ws, prec, isol, GFMCold, roFL) => grassFuelMoisture(temp, rh, ws, prec, isol, GFMCold, timeStep, roFL) + ); + terra.names(mcRaster, "MC"); + + const gfmcRaster = terra.lapp(mcRaster, (mc) => grassFuelMoistureCode(mc)); + terra.names(gfmcRaster, "GFMC"); + + // Return requested output + switch (out) { + case "ALL": + return [...Object.values(input), gfmcRaster, mcRaster]; + case "GFMC": + return gfmcRaster; + case "MC": + return mcRaster; + case "GFMCandMC": + default: + return [gfmcRaster, mcRaster]; + } +} + +// Example usage +const testInput: GFMInput = { + temp: terra.rast({ nrows: 25, ncols: 25, vals: Array(25 * 25).fill(20) }), + rh: terra.rast({ nrows: 25, ncols: 25, vals: Array(25 * 25).fill(50) }), + ws: terra.rast({ nrows: 25, ncols: 25, vals: Array(25 * 25).fill(10) }), + prec: terra.rast({ nrows: 25, ncols: 25, vals: Array(25 * 25).fill(0) }), + isol: terra.rast({ nrows: 25, ncols: 25, vals: Array(25 * 25).fill(0.5) }) +}; + +console.log(gfmcRaster(testInput, 85, 1, 0.3, "GFMCandMC")); diff --git a/src/fwi/grass_fuel_moisture.ts b/src/fwi/grass_fuel_moisture.ts new file mode 100644 index 0000000..4179e98 --- /dev/null +++ b/src/fwi/grass_fuel_moisture.ts @@ -0,0 +1,75 @@ +// Define the constant used in the calculation +const FFMC_COEFFICIENT = 0.5; // Adjust this value to match the one used in R + +/** + * Calculates the grass fuel moisture content. + * + * @param temp - Temperature (in degrees Celsius) + * @param rh - Relative Humidity (percentage) + * @param ws - Wind Speed (km/h) + * @param prec - Precipitation (mm) + * @param isol - Insolation (kW/m^2) + * @param GFMCold - Yesterday's Grass Foliar Moisture Content + * @param roFL - The nominal fuel load of the fine fuel layer (default 0.3 kg/m^2) + * @param timeStep - Time step (default 1 hour) + * @returns The calculated moisture content + */ +function grassFuelMoisture( + temp: number, + rh: number, + ws: number, + prec: number, + isol: number, + GFMCold: number, + roFL: number = 0.3, + timeStep: number = 1 +): number { + // Calculate previous moisture code + let MCold = FFMC_COEFFICIENT * ((101 - GFMCold) / (59.5 + GFMCold)); + + // Calculate the moisture content of the layer in % after rainfall + let MCr = prec > 0 ? MCold + 100 * (prec / roFL) : MCold; + MCr = Math.min(MCr, 250); // Constrain to 250 + MCold = MCr; + + // Calculate Fuel temperature + const Tf = temp + 35.07 * isol * Math.exp(-0.06215 * ws); + + // Calculate Saturation Vapour Pressure (Baumgartner et al. 1982) + const eS_T = 6.107 * Math.pow(10, (7.5 * temp) / (237 + temp)); + const eS_Tf = 6.107 * Math.pow(10, (7.5 * Tf) / (237 + Tf)); + + // Calculate Fuel Level Relative Humidity + const RH_f = rh * (eS_T / eS_Tf); + + // Calculate Equilibrium Moisture Content for Drying phase + const EMC_D = ((1.62 * Math.pow(RH_f, 0.532) + 13.7 * Math.exp((RH_f - 100) / 13.0)) + + 0.27 * (26.7 - Tf) * (1 - Math.exp(-0.115 * RH_f))); + + // Calculate Equilibrium Moisture Content for Wetting phase + const EMC_W = ((1.42 * Math.pow(RH_f, 0.512) + 12.0 * Math.exp((RH_f - 100) / 18.0)) + + 0.27 * (26.7 - Tf) * (1 - Math.exp(-0.115 * RH_f))); + + // RH in terms of RH/100 for desorption + let Rf = MCold > EMC_D ? RH_f / 100 : rh; + // RH in terms of 1-RH/100 for absorption + Rf = MCold < EMC_W ? (100 - RH_f) / 100 : Rf; + + // Calculate Inverse Response time of grass (hours) + const K_GRASS = 0.389633 * Math.exp(0.0365 * Tf) * (0.424 * (1 - Math.pow(Rf, 1.7)) + 0.0694 * + Math.sqrt(ws) * (1 - Math.pow(Rf, 8))); + + // Fuel is drying, calculate Moisture Content + let MC0 = MCold > EMC_D + ? EMC_D + (MCold - EMC_D) * Math.exp(-1.0 * Math.log(10.0) * K_GRASS * timeStep) + : MCold; + + // Fuel is wetting, calculate moisture content + MC0 = MCold < EMC_W + ? EMC_W + (MCold - EMC_W) * Math.exp(-1.0 * Math.log(10.0) * K_GRASS * timeStep) + : MC0; + + return MC0; +} + +export {grassFuelMoisture}; \ No newline at end of file diff --git a/src/fwi/grass_fuel_moisture_code.ts b/src/fwi/grass_fuel_moisture_code.ts new file mode 100644 index 0000000..713934d --- /dev/null +++ b/src/fwi/grass_fuel_moisture_code.ts @@ -0,0 +1,19 @@ +// Define the constant used in the calculation +const FFMC_COEFFICIENT = 0.5; // Adjust this value to match the one used in R + +/** + * Calculates the Grass Fuel Moisture Code (GFMC). + * + * @param MC0 - An output from the mcCalc function + * @returns The calculated GFMC + */ +function grassFuelMoistureCode(MC0: number): number { + // Calculate GFMC + const GFMC0 = 59.5 * ((250 - MC0) / (FFMC_COEFFICIENT + MC0)); + return GFMC0; +} + +// Example usage +const MC0 = 100; // Example value +const GFMC = grassFuelMoistureCode(MC0); +console.log(GFMC); // Output the GFMC diff --git a/src/fwi/hffmc.ts b/src/fwi/hffmc.ts new file mode 100644 index 0000000..c43abb7 --- /dev/null +++ b/src/fwi/hffmc.ts @@ -0,0 +1,188 @@ +const FFMC_COEFFICIENT = 0.5; // Adjust this value to match the one used in R + +interface WeatherInput { + temp: number[]; + rh: number[]; + ws: number[]; + prec: number[]; + hr?: number[]; + bui?: number[]; + id?: number[]; +} + +/** + * Calculate the hourly fine fuel moisture code. + * + * @param temp - Temperature + * @param ws - Wind Speed + * @param rh - Relative Humidity + * @param prec - Precipitation + * @param Fo - Previous FFMC + * @param t0 - Time step + * @returns The calculated FFMC + */ +function hourlyFineFuelMoistureCode( + temp: number, + ws: number, + rh: number, + prec: number, + Fo: number, + t0: number +): number { + // Implement the detailed calculation logic here + // Placeholder implementation + return Fo; // Placeholder return, replace with actual calculation logic +} + +/** + * Calculate the initial spread index. + * + * @param f - FFMC + * @param ws - Wind Speed + * @param flag - Boolean flag (not used here) + * @returns The calculated ISI + */ +function initialSpreadIndex(f: number[], ws: number[], flag: boolean): number[] { + // Implement the detailed calculation logic here + return f.map((ffmc, index) => 0); // Placeholder return, replace with actual calculation logic +} + +/** + * Calculate the fire weather index. + * + * @param isi - Initial Spread Index + * @param bui - Build Up Index + * @returns The calculated FWI + */ +function fireWeatherIndex(isi: number[], bui: number[]): number[] { + // Implement the detailed calculation logic here + return isi.map((isiValue, index) => 0); // Placeholder return, replace with actual calculation logic +} + +/** + * Main function to calculate the Hourly Fine Fuel Moisture Code (hffmc). + * + * @param input - Input weather data + * @param ffmcOld - Initial FFMC + * @param timeStep - Time step (default is 1 hour) + * @param calcStep - Calculate time step between observations (default is false) + * @param batch - Iterative computation flag (default is true) + * @param hourlyFWI - Compute hourly ISI, FWI, and DSR (default is false) + * @returns A vector of hourly or sub-hourly FFMC values + */ +function hffmc( + input: WeatherInput, + ffmcOld: number | number[] = 85, + timeStep: number = 1, + calcStep: boolean = false, + batch: boolean = true, + hourlyFWI: boolean = false +): number[] | any { + let t0 = timeStep; + const names = Object.keys(input).map(name => name.toLowerCase()); + let n: number; + + if (batch) { + if (names.includes('id')) { + const ids = Array.from(new Set(input.id)); + n = ids.length; + if (ids.length !== n) { + throw new Error( + "Multiple stations have to start and end at the same dates/time, and the data must be sorted by date/time and id" + ); + } + } else { + n = 1; + } + } else { + n = input.temp.length; + } + + let Fo = Array.isArray(ffmcOld) && ffmcOld.length === 1 && n > 1 ? Array(n).fill(ffmcOld[0]) : ffmcOld; + + if (calcStep) { + const hr = input.hr; + if (!hr) { + console.warn("hour value is missing!"); + } + } + + const requiredCols = ['temp', 'prec', 'ws', 'rh']; + + for (const col of requiredCols) { + if (!names.includes(col)) { + throw new Error(`${col} is missing!`); + } + } + + if (input.prec.some(prec => prec < 0)) { + console.warn("precipitation (prec) cannot be negative!"); + } + if (input.ws.some(ws => ws < 0)) { + console.warn("wind speed (ws) cannot be negative!"); + } + if (input.rh.some(rh => rh < 0)) { + console.warn("relative humidity (rh) cannot be negative!"); + } + if (input.rh.length % n !== 0) { + console.warn("input do not match with number of weather stations"); + } + + const n0 = input.rh.length / n; + let f: number[] = []; + + for (let i = 0; i < n0; i++) { + const k = ((i * n)); // zero-based index + + if (calcStep && i > 0 && input.hr) { + t0 = n0 > 1 ? input.hr[k] - input.hr[k - n] : t0; + t0 = t0 === -23 ? 1 : t0; + t0 = t0 < 0 ? -t0 : t0; + } + + const f1 = hourlyFineFuelMoistureCode( + input.temp[k], + input.ws[k], + input.rh[k], + input.prec[k], + Array.isArray(Fo) ? Fo[i % Fo.length] : Fo, // Ensure Fo is correctly indexed + t0 + ); + + Fo = Array.isArray(Fo) ? [...Fo.slice(0, i), f1, ...Fo.slice(i + 1)] : f1; + f.push(f1); + } + + if (hourlyFWI) { + const bui = input.bui; + const ws = input.ws; + + if (!bui) { + console.warn("Daily BUI is required to calculate hourly FWI"); + } else { + const isi = initialSpreadIndex(f, ws, false); + const fwi = fireWeatherIndex(isi, bui); + const dsr = fwi.map(fwiValue => 0.0272 * Math.pow(fwiValue, 1.77)); + const output = { + ...input, + ffmc: f, + isi, + fwi, + dsr + }; + return output; + } + } else { + return f; + } +} + +// Example usage +const weatherData: WeatherInput = { + temp: [20, 21, 22], + rh: [50, 55, 60], + ws: [10, 12, 14], + prec: [0, 0, 0] +}; +const result = hffmc(weatherData); +console.log(result); diff --git a/src/fwi/hffmcRaster.ts b/src/fwi/hffmcRaster.ts new file mode 100644 index 0000000..aa6dd9f --- /dev/null +++ b/src/fwi/hffmcRaster.ts @@ -0,0 +1,89 @@ +import * as geotiff from 'geotiff'; +import { rasterCalculator } from 'geoblaze'; +import { computeFFMC, computeISI, computeFWI, computeDSR } from './fireWeatherIndexFunctions'; + +interface WeatherStream { + temp: number[][]; + rh: number[][]; + ws: number[][]; + prec: number[][]; + bui?: number[][]; +} + +interface Raster { + [key: string]: number[][]; +} + +async function hffmcRaster( + weatherstream: WeatherStream, + ffmc_old: number | number[][] = 85, + timeStep: number = 1, + hourlyFWI: boolean = false +): Promise { + const { temp, rh, ws, prec, bui } = weatherstream; + + if (!temp || !rh || !ws || !prec) { + throw new Error("One of the required inputs (temp, rh, ws, prec) is missing!"); + } + + let ffmcOldRaster: number[][]; + + if (typeof ffmc_old === 'number') { + ffmcOldRaster = Array(temp.length).fill(Array(temp[0].length).fill(ffmc_old)); + } else { + ffmcOldRaster = ffmc_old; + } + + const ffmcRaster = rasterCalculator([temp, rh, ws, prec, ffmcOldRaster], computeFFMC, { timeStep }); + + if (hourlyFWI && bui) { + const isiRaster = rasterCalculator([ffmcRaster, ws], computeISI, { fbpMod: false }); + const fwiRaster = rasterCalculator([isiRaster, bui], computeFWI); + const dsrRaster = fwiRaster.map(row => row.map(value => 0.0272 * Math.pow(value, 1.77))); + + return { + hffmc: ffmcRaster, + hisi: isiRaster, + hfwi: fwiRaster, + hdsr: dsrRaster + }; + } else if (hourlyFWI && !bui) { + throw new Error("Daily BUI is required to calculate hourly FWI"); + } else { + return { hffmc: ffmcRaster }; + } +} + +// Example placeholder functions for fire weather index calculations +function computeFFMC(temp: number, rh: number, ws: number, prec: number, ffmcOld: number, options: { timeStep: number }): number { + // Placeholder function logic here + return ffmcOld; // Replace with actual FFMC calculation +} + +function computeISI(ffmc: number, ws: number, options: { fbpMod: boolean }): number { + // Placeholder function logic here + return ffmc; // Replace with actual ISI calculation +} + +function computeFWI(isi: number, bui: number): number { + // Placeholder function logic here + return isi; // Replace with actual FWI calculation +} + +// Test the function +(async () => { + const weatherData: WeatherStream = { + temp: [[20, 21], [22, 23]], + rh: [[50, 55], [60, 65]], + ws: [[10, 12], [14, 16]], + prec: [[0, 0], [0, 0]], + bui: [[50, 50], [50, 50]] // Optional + }; + + try { + const result = await hffmcRaster(weatherData, 85, 1, true); + console.log(result); + } catch (error) { + console.error(error); + } +})(); diff --git a/src/fwi/hourly_fine_fuel_moisture_code.ts b/src/fwi/hourly_fine_fuel_moisture_code.ts new file mode 100644 index 0000000..4d2ab0d --- /dev/null +++ b/src/fwi/hourly_fine_fuel_moisture_code.ts @@ -0,0 +1,165 @@ +interface WeatherData { + temp: number[]; + rh: number[]; + ws: number[]; + prec: number[]; + hr?: number[]; + bui?: number[]; + id?: number[]; + } + + interface FFMCOptions { + ffmcOld?: number; + timeStep?: number; + calcStep?: boolean; + batch?: boolean; + hourlyFWI?: boolean; + } + + function hourlyFineFuelMoistureCode( + temp: number[], + rh: number[], + ws: number[], + prec: number[], + Fo: number[], + t0: number + ): number[] { + const FFMC_COEFFICIENT = 59.5; + const mo = Fo.map(f => FFMC_COEFFICIENT * (101 - f) / (59.5 + f)); + const rf = prec; + + const mr = mo.map((moVal, i) => { + if (moVal <= 150) { + return moVal + 42.5 * rf[i] * Math.exp(-100 / (251 - moVal)) * (1 - Math.exp(-6.93 / rf[i])); + } else { + return ( + moVal + + 42.5 * rf[i] * Math.exp(-100 / (251 - moVal)) * (1 - Math.exp(-6.93 / rf[i])) + + 0.0015 * Math.pow(moVal - 150, 2) * Math.sqrt(rf[i]) + ); + } + }); + + const mrCapped = mr.map(m => Math.min(m, 250)); + const moUpdated = prec.map((p, i) => (p > 0 ? mrCapped[i] : mo[i])); + + const Ed = rh.map((rhVal, i) => + 0.942 * Math.pow(rhVal, 0.679) + + 11 * Math.exp((rhVal - 100) / 10) + + 0.18 * (21.1 - temp[i]) * (1 - Math.exp(-0.115 * rhVal)) + ); + + const ko = rh.map((rhVal, i) => + 0.424 * (1 - Math.pow(rhVal / 100, 1.7)) + + 0.0694 * Math.sqrt(ws[i]) * (1 - Math.pow(rhVal / 100, 8)) + ); + + const kd = ko.map((k, i) => k * 0.0579 * Math.exp(0.0365 * temp[i])); + const md = Ed.map((ed, i) => ed + (moUpdated[i] - ed) * Math.pow(10, -kd[i] * t0)); + + const Ew = rh.map((rhVal, i) => + 0.618 * Math.pow(rhVal, 0.753) + + 10 * Math.exp((rhVal - 100) / 10) + + 0.18 * (21.1 - temp[i]) * (1 - Math.exp(-0.115 * rhVal)) + ); + + const k1 = rh.map((rhVal, i) => + 0.424 * (1 - Math.pow((100 - rhVal) / 100, 1.7)) + + 0.0694 * Math.sqrt(ws[i]) * (1 - Math.pow((100 - rhVal) / 100, 8)) + ); + + const kw = k1.map((k, i) => k * 0.0579 * Math.exp(0.0365 * temp[i])); + const mw = Ew.map((ew, i) => ew - (ew - moUpdated[i]) * Math.pow(10, -kw[i] * t0)); + + const m = moUpdated.map((moVal, i) => + moVal > Ed[i] ? md[i] : moVal < Ew[i] ? mw[i] : moVal + ); + + const FoNew = m.map(mVal => 59.5 * (250 - mVal) / (FFMC_COEFFICIENT + mVal)).map(FoVal => Math.max(FoVal, 0)); + + return FoNew; + } + + function hffmc( + input: WeatherData, + options: FFMCOptions = {} + ): number[] | (WeatherData & { ffmc: number[]; isi?: number[]; fwi?: number[]; dsr?: number[] }) { + const { + ffmcOld = 85, + timeStep = 1, + calcStep = false, + batch = true, + hourlyFWI = false + } = options; + + const { temp, rh, ws, prec, hr, bui, id } = input; + + const n = batch ? (id ? new Set(id).size : 1) : temp.length; + let Fo = Array(n).fill(ffmcOld); + + const f: number[] = []; + + const n0 = temp.length / n; + + for (let i = 0; i < n0; i++) { + const k = Array.from({ length: n }, (_, j) => i * n + j); + const t0 = calcStep && i > 0 && hr ? hr[k[0]] - hr[k[0] - n] : timeStep; + const f1 = hourlyFineFuelMoistureCode( + k.map(idx => temp[idx]), + k.map(idx => rh[idx]), + k.map(idx => ws[idx]), + k.map(idx => prec[idx]), + Fo, + t0 + ); + Fo = f1; + f.push(...f1); + } + + if (hourlyFWI && bui) { + const isi = f.map((ffmc, i) => initialSpreadIndex(ffmc, ws[i])); + const fwi = isi.map((isi, i) => fireWeatherIndex(isi, bui[i])); + const dsr = fwi.map(fwi => 0.0272 * Math.pow(fwi, 1.77)); + + return { + ...input, + ffmc: f, + isi, + fwi, + dsr + }; + } else if (hourlyFWI && !bui) { + throw new Error("Daily BUI is required to calculate hourly FWI"); + } else { + return f; + } + } + + function initialSpreadIndex(ffmc: number, ws: number): number { + // Placeholder function logic here + return ffmc; // Replace with actual ISI calculation + } + + function fireWeatherIndex(isi: number, bui: number): number { + // Placeholder function logic here + return isi; // Replace with actual FWI calculation + } + + // Example test case +/* const weatherData: WeatherData = { + temp: [20, 21, 22, 23], + rh: [50, 55, 60, 65], + ws: [10, 12, 14, 16], + prec: [0, 0, 0, 0], + hr: [1, 2, 3, 4], + bui: [50, 50, 50, 50], + id: [1, 1, 1, 1] + }; + + try { + const result = hffmc(weatherData, { ffmcOld: 85, timeStep: 1, hourlyFWI: true }); + console.log(result); + } catch (error) { + console.error(error); + } + */ \ No newline at end of file diff --git a/src/fwi/initial_spread_index.ts b/src/fwi/initial_spread_index.ts new file mode 100644 index 0000000..faf209e --- /dev/null +++ b/src/fwi/initial_spread_index.ts @@ -0,0 +1,55 @@ +/** + * Initial Spread Index Calculator + * + * @description Computes the Initial Spread Index From the FWI System. Equations + * are from Van Wagner (1985) as listed below, except for the modification for + * fbp taken from FCFDG (1992). + * + * Equations and FORTRAN program for the Canadian Forest Fire + * Weather Index System. 1985. Van Wagner, C.E.; Pickett, T.L. + * Canadian Forestry Service, Petawawa National Forestry + * Institute, Chalk River, Ontario. Forestry Technical Report 33. + * 18 p. + * + * Forestry Canada Fire Danger Group (FCFDG) (1992). Development and + * Structure of the Canadian Forest Fire Behavior Prediction System. + * Technical Report ST-X-3, Forestry Canada, Ottawa, Ontario. + * + * @param ffmc Fine Fuel Moisture Code + * @param ws Wind Speed (km/h) + * @param fbpMod TRUE/FALSE if using the fbp modification at the extreme end + * + * @returns ISI - Initial Spread Index + */ + +function initialSpreadIndex(ffmc: number, ws: number, fbpMod: boolean = false): number { + const FFMC_COEFFICIENT = 59.5; + + // Eq. 10 - Moisture content + const fm = FFMC_COEFFICIENT * (101 - ffmc) / (59.5 + ffmc); + + // Eq. 24 - Wind Effect + // the ifelse, also takes care of the ISI modification for the fbp functions + // This modification is Equation 53a in FCFDG (1992) + const fW = ws >= 40 && fbpMod + ? 12 * (1 - Math.exp(-0.0818 * (ws - 28))) + : Math.exp(0.05039 * ws); + + // Eq. 25 - Fine Fuel Moisture + const fF = 91.9 * Math.exp(-0.1386 * fm) * (1 + Math.pow(fm, 5.31) / 49300000); + + // Eq. 26 - Spread Index Equation + const isi = 0.208 * fW * fF; + + return isi; + } + + export { initialSpreadIndex }; + + // Example usage +/* const ffmc = 85; + const ws = 30; + const fbpMod = false; + const isi = initialSpreadIndex(ffmc, ws, fbpMod); + console.log(`Initial Spread Index: ${isi}`); */ + \ No newline at end of file diff --git a/src/fwi/length_to_breadth.ts b/src/fwi/length_to_breadth.ts new file mode 100644 index 0000000..c0cfd78 --- /dev/null +++ b/src/fwi/length_to_breadth.ts @@ -0,0 +1,45 @@ +/** + * Length-to-Breadth ratio + * + * @description Computes the Length to Breadth ratio of an elliptically shaped + * fire. Equations are from FCFDG (1992) except for errata 80 from + * Wotton et. al. (2009). + * + * All variables names are laid out in the same manner as Forestry Canada + * Fire Danger Group (FCFDG) (1992). Development and Structure of the + * Canadian Forest Fire Behavior Prediction System. Technical Report + * ST-X-3, Forestry Canada, Ottawa, Ontario. + * + * Wotton, B.M., Alexander, M.E., Taylor, S.W. 2009. Updates and revisions to + * the 1992 Canadian forest fire behavior prediction system. Nat. Resour. + * Can., Can. For. Serv., Great Lakes For. Cent., Sault Ste. Marie, Ontario, + * Canada. Information Report GLC-X-10, 45p. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param WSV The Wind Speed (km/h) + * + * @returns Length to Breadth ratio value + */ + +function lengthToBreadth(FUELTYPE: string, WSV: number): number { + let LB: number; + + if (["O1A", "O1B"].includes(FUELTYPE)) { + // Correction to original Equation 80 is made here + // Eq. 80a / 80b from Wotton 2009 + LB = WSV >= 1.0 ? 1.1 * Math.pow(WSV, 0.464) : 1.0; // Eq. 80/81 + } else { + LB = 1.0 + Math.pow(8.729 * (1 - Math.exp(-0.030 * WSV)), 2.155); // Eq. 79 + } + + return LB; + } + + // Example usage +/* const fuelType = "O1A"; + const windSpeed = 10; + const lbRatio = lengthToBreadth(fuelType, windSpeed); + console.log(`Length-to-Breadth ratio: ${lbRatio}`); */ + + export { lengthToBreadth }; + \ No newline at end of file diff --git a/src/fwi/length_to_breadth_at_time.ts b/src/fwi/length_to_breadth_at_time.ts new file mode 100644 index 0000000..b17230c --- /dev/null +++ b/src/fwi/length_to_breadth_at_time.ts @@ -0,0 +1,47 @@ +/** + * Length-to-Breadth ratio at time t + * + * @description Computes the Length to Breadth ratio of an elliptically shaped + * fire at elapsed time since ignition. Equations are from FCFDG (1992) + * and Wotton et. al. (2009), and are marked as such. + * + * All variables names are laid out in the same manner as Forestry Canada + * Fire Danger Group (FCFDG) (1992). Development and Structure of the + * Canadian Forest Fire Behavior Prediction System. Technical Report + * ST-X-3, Forestry Canada, Ottawa, Ontario. + * + * Wotton, B.M., Alexander, M.E., Taylor, S.W. 2009. Updates and revisions to + * the 1992 Canadian forest fire behavior prediction system. Nat. Resour. + * Can., Can. For. Serv., Great Lakes For. Cent., Sault Ste. Marie, Ontario, + * Canada. Information Report GLC-X-10, 45p. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param LB Length to Breadth ratio + * @param HR Time since ignition (hours) + * @param CFB Crown Fraction Burned + * + * @returns Length to Breadth ratio at time since ignition + */ + +function lengthToBreadthAtTime(FUELTYPE: string, LB: number, HR: number, CFB: number): number { + // Eq. 72 (FCFDG 1992) - alpha constant value, dependent on fuel type + const alpha: number = ["C1", "O1A", "O1B", "S1", "S2", "S3", "D1"].includes(FUELTYPE) + ? 0.115 + : 0.115 - 18.8 * Math.pow(CFB, 2.5) * Math.exp(-8 * CFB); + + // Eq. 81 (Wotton et.al. 2009) - LB at time since ignition + const LBt: number = (LB - 1) * (1 - Math.exp(-alpha * HR)) + 1; + + return LBt; + } + +export { lengthToBreadthAtTime }; + + // Example usage + /* const fuelType = "C1"; + const lengthBreadthRatio = 2.5; + const hoursSinceIgnition = 5; + const crownFractionBurned = 0.3; + const lbAtTime = lengthToBreadthAtTime(fuelType, lengthBreadthRatio, hoursSinceIgnition, crownFractionBurned); + console.log(`Length-to-Breadth ratio at time: ${lbAtTime}`); */ + \ No newline at end of file diff --git a/src/fwi/lros.ts b/src/fwi/lros.ts new file mode 100644 index 0000000..bc22b5a --- /dev/null +++ b/src/fwi/lros.ts @@ -0,0 +1,171 @@ +/* #' Line-based input for Simard Rate of Spread and Direction +#' +#' @description \code{lros} is used to calculate the rate of spread and +#' direction given one set of three point-based observations of fire arrival +#' time. The function requires that the user specify the time that the fire +#' crossed each point, along with the measured lengths between each pair of +#' observational points, and a reference bearing (one specified side of the +#' triangle). This function allows quick input of a dataframe specifying one or +#' many triangles. +#' +#' \code{lros} Allows R users to calculate the rate of spread and direction of +#' a fire across a triangle, given three time measurements and details about +#' the orientation and distance between observational points. The algorithm is +#' based on the description from Simard et al. (1984). See \code{pros} for more +#' information. +#' +#' The functions require the user to arrange the input dataframe so that each +#' triangle of interest is identified based on a new row in the dataframe. The +#' input format forces the user to identify the triangles, one triangle per row +#' of input dataframe. Very complex arrangements of field plot layouts are +#' possible, and the current version of these functions do not attempt to +#' determine each triangle of interest automatically. +#' +#' @param input A dataframe containing input variables of time fire front +#' crossed points 1, 2, 3, and latitude/longitude for those same points. +#' Variable names have to be the same as in the following list, but they are +#' case insensitive. The order in which the input variables are entered is not +#' important. +#' +#' \tabular{lll}{ +#' \var{T1} +#' \tab (required) +#' \tab Time that the fire front crossed point 1.\cr +#' \tab\tab Time entered in fractional format. \cr +#' \tab\tab Output ROS will depend on the level \cr +#' \tab\tab of precision entered \cr +#' \tab\tab (minute, second, decisecond)\cr +#' \var{T2} +#' \tab (required) +#' \tab Time that the fire front crossed point 2.\cr +#' \tab\tab Time entered in fractional format. \cr +#' \tab\tab Output ROS will depend on the level \cr +#' \tab\tab of precision entered \cr +#' \tab\tab (minute, second, decisecond)\cr +#' \var{T3} +#' \tab (required) +#' \tab Time that the fire front crossed point 3. \cr +#' \tab\tab Time entered in fractional format. \cr +#' \tab\tab Output ROS will depend on the level \cr +#' \tab\tab of precision entered \cr +#' \tab\tab (minute, second, decisecond)\cr +#' \var{LengthT1T2} +#' \tab (required) \tab Length between each pair of\cr +#' \tab\tab observation points T1 and T2 (subscripts \cr +#' \tab\tab denote time-ordered pairs). (meters)\cr +#' \var{LengthT2T3} +#' \tab (required) +#' \tab Length between each pair of\cr +#' \tab\tab observation points T2 and T3 (subscripts \cr +#' \tab\tab denote time-ordered pairs). (meters)\cr +#' \var{LengthT1T3} +#' \tab (required) +#' \tab Length between each pair of\cr +#' \tab\tab observation points T1 and T3 (subscripts \cr +#' \tab\tab denote time-ordered pairs). (meters)\cr +#' \var{BearingT1T2} +#' \tab (required) +#' \tab Reference bearing. For reference,\cr +#' \tab\tab North = 0, West = -90, East = 90 (degrees)\cr +#' \var{BearingT1T3} +#' \tab (required) +#' \tab Reference bearing. For reference,\cr +#' \tab\tab North = 0, West = -90, East = 90 (degrees)\cr +#' } +#' @return \code{lros} returns a dataframe which includes the rate of spread +#' and spread direction. Output units depend on the user’s inputs for +#' distance (typically meters) and time (seconds or minutes). +#' @author Tom Schiks, Xianli Wang, Alan Cantin +#' @seealso \code{\link{pros}}, +#' @references 1. Simard, A.J., Eenigenburg, J.E., Adams, K.B., Nissen, R.L., +#' Deacon, and Deacon, A.G. 1984. A general procedure for sampling and +#' analyzing wildland fire spread. +#' +#' 2. Byram, G.M. 1959. Combustion of forest fuels. In: Davis, K.P. Forest Fire +#' Control and Use. McGraw-Hill, New York. +#' +#' 3. Curry, J.R., and Fons, W.L. 1938. Rate of spread of surface fires in the +#' Ponderosa Pine Type of California. Journal of Agricultural Research 57(4): +#' 239-267. +#' +#' 4. Simard, A.J., Deacon, A.G., and Adams, K.B. 1982. Nondirectional sampling +#' wildland fire spread. Fire Technology: 221-228. +#' @keywords ros */ + +interface LrosInput { + T1: number; + LengthT1T2: number; + T2: number; + LengthT1T3: number; + T3: number; + LengthT2T3: number; + BearingT1T2: number; + BearingT1T3: number; + } + + interface LrosOutput { + Ros: number; + Direction: number; + } + + function direction(bearingT1T2: number, bearingT1T3: number, thetaAdeg: number): number { + // Placeholder for the actual direction calculation + // Replace this with the actual logic from the R function .direction + return bearingT1T2 + bearingT1T3 + thetaAdeg; // Adjust this logic accordingly + } + + function lros(input: LrosInput[]): LrosOutput[] { + // Check if input is an array of objects with the required properties + const requiredColumns = [ + "T1", "LengthT1T2", "T2", "LengthT1T3", "T3", + "LengthT2T3", "BearingT1T2", "BearingT1T3" + ]; + + input.forEach(row => { + const missingColumns = requiredColumns.filter(col => !(col in row)); + if (missingColumns.length > 0) { + throw new Error(`cffdrs::lros Column ${missingColumns.join(", ")} is required in column list.`); + } + }); + + return input.map(row => { + const AngleArad = Math.acos( + (Math.pow(row.LengthT1T3, 2) + Math.pow(row.LengthT1T2, 2) - Math.pow(row.LengthT2T3, 2)) / + (2 * row.LengthT1T3 * row.LengthT1T2) + ); + + const AngleAdeg = (AngleArad * 180) / Math.PI; + + const ThetaArad = Math.atan( + ((row.T3 - row.T1) / (row.T2 - row.T1)) * + (row.LengthT1T2 / (row.LengthT1T3 * Math.sin(AngleArad))) - + (1 / Math.tan(AngleArad)) + ); + + const ThetaAdeg = (ThetaArad * 180) / Math.PI; + + const DIR = direction(row.BearingT1T2, row.BearingT1T3, ThetaAdeg); + + const ROS = (row.LengthT1T2 * Math.cos(ThetaArad)) / (row.T2 - row.T1); + + return { Ros: ROS, Direction: DIR }; + }); + } + + // Example usage + const inputData: LrosInput[] = [ + { + T1: 0, + LengthT1T2: 24.5, + T2: 22.6, + LengthT1T3: 120, + T3: 20.0, + LengthT2T3: 35, + BearingT1T2: 90, + BearingT1T3: -90 + } + ]; + + const outputData = lros(inputData); + console.log(outputData); + \ No newline at end of file diff --git a/src/fwi/overwinter_drought_code.ts b/src/fwi/overwinter_drought_code.ts new file mode 100644 index 0000000..8b1367d --- /dev/null +++ b/src/fwi/overwinter_drought_code.ts @@ -0,0 +1,68 @@ +/* #' Overwintering Drought Code +#' +#' @description \code{overwinter_drought_code} calculates an initial or season +#' starting Drought Code (DC) value based on a standard method of overwintering +#' the Drought Code (Lawson and Armitage 2008). This method uses the final DC +#' value from previous year, over winter precipitation and estimates of how much +#' over-winter precipitation 'refills' the moisture in this fuel layer. This +#' function could be used for either one weather station or for multiple weather +#' stations. +#' +#' Of the three fuel moisture codes (i.e. FFMC, DMC and DC) making up the FWI +#' System, only the DC needs to be considered in terms of its values carrying +#' over from one fire season to the next. In Canada both the FFMC and the DMC +#' are assumed to reach moisture saturation from overwinter precipitation at or +#' before spring melt; this is a reasonable assumption and any error in these +#' assumed starting conditions quickly disappears. If snowfall (or other +#' overwinter precipitation) is not large enough however, the fuel layer +#' tracked by the Drought Code may not fully reach saturation after spring snow +#' melt; because of the long response time in this fuel layer (53 days in +#' standard conditions) a large error in this spring starting condition can +#' affect the DC for a significant portion of the fire season. In areas where +#' overwinter precipitation is 200 mm or more, full moisture recharge occurs +#' and DC overwintering is usually unnecessary. More discussion of +#' overwintering and fuel drying time lag can be found in Lawson and Armitage +#' (2008) and Van Wagner (1985). +#' +#' @param DCf Final fall DC value from previous year +#' @param rw Winter precipitation (mm) +#' @param a User selected values accounting for carry-over fraction (view table +#' below) +#' @param b User selected values accounting for wetting efficiency fraction +#' (view table below) +#' @return \code{overwinter_drought_code} returns either a single value or a +#' vector of wDC values. +#' @author Xianli Wang, Mike Wotton, Alan Cantin, and Mike Flannigan +#' @seealso \code{\link{fwi}}, \code{\link{fire_season}} +#' @references Lawson B.D. and Armitage O.B. 2008. Weather Guide for the +#' Canadian Forest Fire Danger Rating System. Natural Resources Canada, +#' Canadian Forest Service, Northern Forestry Centre, Edmonton, Alberta. 84 p. +#' \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/29152.pdf} +#' +#' Van Wagner, C.E. 1985. Drought, timelag and fire danger rating. Pages +#' 178-185 in L.R. Donoghue and R.E. Martin, eds. Proc. 8th Conf. Fire For. +#' Meteorol., 29 Apr.-3 May 1985, Detroit, MI. Soc. Am. For., Bethesda, MD. +#' \url{https://cfs.nrcan.gc.ca/pubwarehouse/pdfs/23550.pdf} */ + +function overwinterDroughtCode(DCf: number = 100, rw: number = 200, a: number = 0.75, b: number = 0.75): number { + // Eq. 3 - Final fall moisture equivalent of the DC + const Qf = 800 * Math.exp(-DCf / 400); + // Eq. 2 - Starting spring moisture equivalent of the DC + const Qs = a * Qf + b * (3.94 * rw); + // Eq. 4 - Spring start-up value for the DC + let DCs = 400 * Math.log(800 / Qs); + // Constrain DC + DCs = DCs < 15 ? 15 : DCs; + return DCs; + } + + // Example usage + const winterDC1 = overwinterDroughtCode(300, 110); + console.log(winterDC1); + + const winterDC2 = overwinterDroughtCode(300, 110, 1.0, 0.9); + console.log(winterDC2); + + const winterDC3 = [400, 300, 250].map((DCf, i) => overwinterDroughtCode(DCf, [99, 110, 200][i], [0.75, 1.0, 0.75][i], [0.75, 0.9, 0.75][i])); + console.log(winterDC3); + \ No newline at end of file diff --git a/src/fwi/pros.ts b/src/fwi/pros.ts new file mode 100644 index 0000000..3becdbe --- /dev/null +++ b/src/fwi/pros.ts @@ -0,0 +1,106 @@ +import * as geolib from 'geolib'; + +interface InputData { + T1: number; + LONG1: number; + LAT1: number; + T2: number; + LONG2: number; + LAT2: number; + T3: number; + LONG3: number; + LAT3: number; +} + +interface OutputData { + Ros: number; + Direction: number; +} + +function direction(bearingT1T2: number, bearingT1T3: number, ThetaAdeg: number): number { + // Calculate the direction based on the provided bearings and angle + // This is a placeholder implementation; the actual logic should be based on the R code's .direction function + return (bearingT1T2 + bearingT1T3 + ThetaAdeg) % 360; // Adjust this logic as needed +} + +function pros(input: InputData[]): OutputData[] { + // Force uppercase to the column names + input.forEach((entry) => { + for (const key in entry) { + if (entry.hasOwnProperty(key)) { + const upperKey = key.toUpperCase(); + if (upperKey !== key) { + entry[upperKey as keyof InputData] = entry[key as keyof InputData]; + delete entry[key as keyof InputData]; + } + } + } + }); + + // Check if the columns exist + const requiredColumns = [ + "T1", "LONG1", "LAT1", "T2", "LONG2", "LAT2", "T3", "LONG3", "LAT3" + ]; + + input.forEach((entry) => { + requiredColumns.forEach((col) => { + if (!(col in entry)) { + throw new Error(`Column ${col} is required in column list.`); + } + }); + }); + + return input.map((entry) => { + const LengthT1T2 = geolib.getDistance( + { latitude: entry.LAT1, longitude: entry.LONG1 }, + { latitude: entry.LAT2, longitude: entry.LONG2 } + ); + const LengthT1T3 = geolib.getDistance( + { latitude: entry.LAT1, longitude: entry.LONG1 }, + { latitude: entry.LAT3, longitude: entry.LONG3 } + ); + const LengthT2T3 = geolib.getDistance( + { latitude: entry.LAT2, longitude: entry.LONG2 }, + { latitude: entry.LAT3, longitude: entry.LONG3 } + ); + + const bearingT1T2 = geolib.getGreatCircleBearing( + { latitude: entry.LAT1, longitude: entry.LONG1 }, + { latitude: entry.LAT2, longitude: entry.LONG2 } + ); + const bearingT1T3 = geolib.getGreatCircleBearing( + { latitude: entry.LAT1, longitude: entry.LONG1 }, + { latitude: entry.LAT3, longitude: entry.LONG3 } + ); + + const AngleArad = Math.acos( + (Math.pow(LengthT1T3, 2) + Math.pow(LengthT1T2, 2) - Math.pow(LengthT2T3, 2)) + / (2 * LengthT1T3 * LengthT1T2) + ); + const AngleAdeg = (AngleArad * 180) / Math.PI; + + const ThetaArad = Math.atan( + ((entry.T3 - entry.T1) / (entry.T2 - entry.T1)) + * (LengthT1T2 / (LengthT1T3 * Math.sin(AngleArad))) + - (1 / Math.tan(AngleArad)) + ); + const ThetaAdeg = (ThetaArad * 180) / Math.PI; + + const DIR = direction(bearingT1T2, bearingT1T3, ThetaAdeg); + const ROS = (LengthT1T2 * Math.cos(ThetaArad)) / (entry.T2 - entry.T1); + + return { Ros: ROS, Direction: DIR }; + }); +} + +// Example usage +const inputData: InputData[] = [ + { + T1: 2, LONG1: -79.701027, LAT1: 43.808872, + T2: 50, LONG2: -79.699650, LAT2: 43.808833, + T3: 120, LONG3: -79.700387, LAT3: 43.809816 + } +]; + +const result = pros(inputData); +console.log(result); diff --git a/src/fwi/rate_of_spread.ts b/src/fwi/rate_of_spread.ts new file mode 100644 index 0000000..d97f69f --- /dev/null +++ b/src/fwi/rate_of_spread.ts @@ -0,0 +1,102 @@ +/* #' Rate of Spread Calculation +#' +#' +#' @description Computes the Rate of Spread prediction based on fuel type and +#' FWI conditions. Equations are from listed FCFDG (1992) and Wotton et. al. +#' (2009), and are marked as such. +#' +#' All variables names are laid out in the same manner as Forestry Canada +#' Fire Danger Group (FCFDG) (1992). Development and Structure of the +#' Canadian Forest Fire Behavior Prediction System." Technical Report +#' ST-X-3, Forestry Canada, Ottawa, Ontario. +#' +#' Wotton, B.M., Alexander, M.E., Taylor, S.W. 2009. Updates and revisions to +#' the 1992 Canadian forest fire behavior prediction system. Nat. Resour. +#' Can., Can. For. Serv., Great Lakes For. Cent., Sault Ste. Marie, Ontario, +#' Canada. Information Report GLC-X-10, 45p. */ + +import { criticalSurfaceIntensity, surfaceFireRateOfSpread, crownFractionBurned } from "./CFBcalc"; +import { intermediateSurfaceRateOfSpreadC6, surfaceRateOfSpreadC6, crownRateOfSpreadC6, crownFractionBurnedC6, rateOfSpreadC6 } from "./C6calc"; +import { buildupEffect } from "./buildup_effect"; + +interface InputParams { + FUELTYPE: string[]; + ISI: number[]; + BUI: number[]; + FMC: number[]; + SFC: number[]; + PC: number[]; + PDF: number[]; + CC: number[]; + CBH: number[]; +} + +interface ROSResult { + ROS: number[]; + CFB: number[]; + CSI: number[]; + RSO: number[]; +} + +function rateOfSpreadExtended(params: InputParams): ROSResult { + const { FUELTYPE, ISI, BUI, FMC, SFC, PC, PDF, CC, CBH } = params; + const NoBUI = Array(ISI.length).fill(-1); + const d = ["C1", "C2", "C3", "C4", "C5", "C6", "C7", "D1", "M1", "M2", "M3", "M4", "S1", "S2", "S3", "O1A", "O1B"]; + const a = [90, 110, 110, 110, 30, 30, 45, 30, 0, 0, 120, 100, 75, 40, 55, 190, 250]; + const b = [0.0649, 0.0282, 0.0444, 0.0293, 0.0697, 0.0800, 0.0305, 0.0232, 0, 0, 0.0572, 0.0404, 0.0297, 0.0438, 0.0829, 0.0310, 0.0350]; + const c0 = [4.5, 1.5, 3.0, 1.5, 4.0, 3.0, 2.0, 1.6, 0, 0, 1.4, 1.48, 1.3, 1.7, 3.2, 1.4, 1.7]; + + const names = Object.fromEntries(d.map((key, i) => [key, i])); + + let RSI = Array(ISI.length).fill(-1); + + RSI = RSI.map((_, i) => { + if (["C1", "C2", "C3", "C4", "C5", "C7", "D1", "S1", "S2", "S3"].includes(FUELTYPE[i])) { + return a[names[FUELTYPE[i]]] * Math.pow(1 - Math.exp(-b[names[FUELTYPE[i]]] * ISI[i]), c0[names[FUELTYPE[i]]]); + } else if (FUELTYPE[i] === "M1") { + return (PC[i] / 100) * rateOfSpread({ FUELTYPE: Array(ISI.length).fill("C2"), ISI, BUI: NoBUI, FMC, SFC, PC, PDF, CC, CBH })[i] + + ((100 - PC[i]) / 100) * rateOfSpread({ FUELTYPE: Array(ISI.length).fill("D1"), ISI, BUI: NoBUI, FMC, SFC, PC, PDF, CC, CBH })[i]; + } else if (FUELTYPE[i] === "M2") { + return (PC[i] / 100) * rateOfSpread({ FUELTYPE: Array(ISI.length).fill("C2"), ISI, BUI: NoBUI, FMC, SFC, PC, PDF, CC, CBH })[i] + + 0.2 * ((100 - PC[i]) / 100) * rateOfSpread({ FUELTYPE: Array(ISI.length).fill("D1"), ISI, BUI: NoBUI, FMC, SFC, PC, PDF, CC, CBH })[i]; + } + return -1; + }); + + const RSI_m3 = RSI.map((rsi, i) => (FUELTYPE[i] === "M3" ? a[names["M3"]] * Math.pow(1 - Math.exp(-b[names["M3"]] * ISI[i]), c0[names["M3"]]) : rsi)); + RSI = RSI.map((rsi, i) => (FUELTYPE[i] === "M3" ? (PDF[i] / 100 * RSI_m3[i] + (1 - PDF[i] / 100) * rateOfSpread({ FUELTYPE: Array(ISI.length).fill("D1"), ISI, BUI: NoBUI, FMC, SFC, PC, PDF, CC, CBH })[i]) : rsi)); + + const RSI_m4 = RSI.map((rsi, i) => (FUELTYPE[i] === "M4" ? a[names["M4"]] * Math.pow(1 - Math.exp(-b[names["M4"]] * ISI[i]), c0[names["M4"]]) : rsi)); + RSI = RSI.map((rsi, i) => (FUELTYPE[i] === "M4" ? (PDF[i] / 100 * RSI_m4[i] + 0.2 * (1 - PDF[i] / 100) * rateOfSpread({ FUELTYPE: Array(ISI.length).fill("D1"), ISI, BUI: NoBUI, FMC, SFC, PC, PDF, CC, CBH })[i]) : rsi)); + + const CF = RSI.map((rsi, i) => (["O1A", "O1B"].includes(FUELTYPE[i]) ? (CC[i] < 58.8 ? 0.005 * (Math.exp(0.061 * CC[i]) - 1) : 0.176 + 0.02 * (CC[i] - 58.8)) : -99)); + RSI = RSI.map((rsi, i) => (["O1A", "O1B"].includes(FUELTYPE[i]) ? a[names[FUELTYPE[i]]] * Math.pow(1 - Math.exp(-b[names[FUELTYPE[i]]] * ISI[i]), c0[names[FUELTYPE[i]]]) * CF[i] : rsi)); + + const CSI = FMC.map((fmc, i) => criticalSurfaceIntensity(fmc, CBH[i])); + const RSO = SFC.map((sfc, i) => surfaceFireRateOfSpread(CSI[i], sfc)); + + RSI = RSI.map((rsi, i) => (FUELTYPE[i] === "C6" ? intermediateSurfaceRateOfSpreadC6(ISI[i]) : rsi)); + const RSC = FUELTYPE.map((ftype, i) => (ftype === "C6" ? crownRateOfSpreadC6(ISI[i], FMC[i]) : NaN)); + + const RSS = RSI.map((rsi, i) => (FUELTYPE[i] === "C6" ? surfaceRateOfSpreadC6(RSI[i], BUI[i]) : buildupEffect(FUELTYPE[i], BUI[i]) * rsi)); + + const CFB = RSI.map((rsi, i) => (FUELTYPE[i] === "C6" ? crownFractionBurnedC6(RSC[i], RSS[i], RSO[i]) : crownFractionBurned(RSS[i], RSO[i]))); + const ROS = RSI.map((rsi, i) => (FUELTYPE[i] === "C6" ? rateOfSpreadC6(RSC[i], RSS[i], CFB[i]) : RSS[i])); + + const constrainedROS = ROS.map(ros => (ros <= 0 ? 0.000001 : ros)); + + return { ROS: constrainedROS, CFB, CSI, RSO }; +} + +function rateOfSpread(params: InputParams): number[] { + const { FUELTYPE, ISI, BUI, FMC, SFC, PC, PDF, CC, CBH } = params; + const rosVars = rateOfSpreadExtended({ FUELTYPE, ISI, BUI, FMC, SFC, PC, PDF, CC, CBH }); + return rosVars.ROS; +} + +function ROScalc(...args: any[]): number[] { + console.warn(".Deprecated: 'rate_of_spread' function is deprecated"); + return rateOfSpread(args[0]); +} + +export { rateOfSpread, InputParams }; \ No newline at end of file diff --git a/src/fwi/rate_of_spread_at_theta.ts b/src/fwi/rate_of_spread_at_theta.ts new file mode 100644 index 0000000..19a971e --- /dev/null +++ b/src/fwi/rate_of_spread_at_theta.ts @@ -0,0 +1,39 @@ +/** + * Rate of Spread at a point along the perimeter calculator + * + * Computes the Rate of Spread at any point along the perimeter of an elliptically shaped fire. + * Equations are from Wotton et. al. (2009). + * + * Wotton, B.M., Alexander, M.E., Taylor, S.W. 2009. Updates and revisions to + * the 1992 Canadian forest fire behavior prediction system. Nat. Resour. + * Can., Can. For. Serv., Great Lakes For. Cent., Sault Ste. Marie, Ontario, + * Canada. Information Report GLC-X-10, 45p. + * + * @param ROS Rate of Spread (m/min) + * @param FROS Flank Fire Rate of Spread (m/min) + * @param BROS Back Fire Rate of Spread (m/min) + * @param THETA Angle in radians + * @returns Rate of spread at point theta (m/min) + */ + +function rateOfSpreadAtTheta(ROS: number, FROS: number, BROS: number, THETA: number): number { + let c1 = Math.cos(THETA); + const s1 = Math.sin(THETA); + c1 = c1 === 0 ? Math.cos(THETA + 0.001) : c1; + + // Eq. 94 - Calculate the Rate of Spread at point THETA + // large equation, view the paper to see a better representation + const ROStheta = (((ROS - BROS) / (2 * c1) + (ROS + BROS) / (2 * c1)) * + ((FROS * c1 * Math.sqrt(FROS * FROS * c1 * c1 + (ROS * BROS) * s1 * s1) - + ((ROS * ROS - BROS * BROS) / 4) * s1 * s1) / + (FROS * FROS * c1 * c1 + ((ROS + BROS) / 2) * ((ROS + BROS) / 2) * s1 * s1))); + + return ROStheta; + } + + function ROSthetaCalc(...args: any[]): number { + console.warn(".Deprecated: 'rate_of_spread_at_theta' function is deprecated"); + return rateOfSpreadAtTheta(args[0], args[1], args[2], args[3]); + } + + export { rateOfSpreadAtTheta }; \ No newline at end of file diff --git a/src/fwi/rate_of_spread_at_time.ts b/src/fwi/rate_of_spread_at_time.ts new file mode 100644 index 0000000..fe72f97 --- /dev/null +++ b/src/fwi/rate_of_spread_at_time.ts @@ -0,0 +1,37 @@ +/** + * Rate of Spread at time t calculation + * + * Computes the Rate of Spread prediction based on fuel type and + * FWI conditions at elapsed time since ignition. Equations are from listed + * FCFDG (1992). + * + * All variables names are laid out in the same manner as Forestry Canada + * Fire Danger Group (FCFDG) (1992). Development and Structure of the + * Canadian Forest Fire Behavior Prediction System." Technical Report + * ST-X-3, Forestry Canada, Ottawa, Ontario. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param ROSeq Equilibrium Rate of Spread (m/min) + * @param HR Time since ignition (hours) + * @param CFB Crown Fraction Burned + * @returns Rate of Spread at time since ignition value + */ + +function rateOfSpreadAtTime(FUELTYPE: string, ROSeq: number, HR: number, CFB: number): number { + // Eq. 72 - alpha constant value, dependent on fuel type + const alpha = ["C1", "O1A", "O1B", "S1", "S2", "S3", "D1"].includes(FUELTYPE) + ? 0.115 + : 0.115 - 18.8 * Math.pow(CFB, 2.5) * Math.exp(-8 * CFB); + + // Eq. 70 - Rate of Spread at time since ignition + const ROSt = ROSeq * (1 - Math.exp(-alpha * HR)); + + return ROSt; + } + + function ROStCalc(...args: any[]): number { + console.warn(".Deprecated: 'rate_of_spread_at_time' function is deprecated"); + return rateOfSpreadAtTime(args[0], args[1], args[2], args[3]); + } + + export { rateOfSpreadAtTime }; \ No newline at end of file diff --git a/src/fwi/sdmc.ts b/src/fwi/sdmc.ts new file mode 100644 index 0000000..c7b46d4 --- /dev/null +++ b/src/fwi/sdmc.ts @@ -0,0 +1,124 @@ +/** + * Sheltered Duff Moisture Code + * + * @description Computes the Sheltered Duff Moisture Code (sDMC) based on daily noon weather + * observations of temperature, relative humidity, wind speed, 24-hour rainfall, + * and a previous day's calculated or estimated value of sDMC. + * + * @param input An array of objects containing input variables of daily noon weather observations. + * Variable names have to be the same as in the following list, but they are case insensitive. + * The order in which the input variables are entered is not important either. + * - temp: Temperature (centigrade) + * - rh: Relative humidity (%) + * - ws: 10-m height wind speed (km/h) + * - prec: 1-hour rainfall (mm) + * - mon: Month of the observations (integer 1-12) + * - day: Day of the observations (integer) + * @param sdmc_old Previous day's value of SDMC. If null, initial SDMC values will be calculated based on the initial DMC. + * @param batch Whether the computation is iterative or single step, default is true. + * @returns A single value or a vector of SDMC values. + */ + +interface WeatherData { + temp: number; + rh: number; + ws: number; + prec: number; + mon?: number; + day?: number; + id?: number; + dmc: number; + } + + function sdmc(input: WeatherData[], sdmc_old: number | number[] | null = null, batch: boolean = true): number[] { + // Convert input array of objects to lowercase keys + const lowerCaseInput = input.map(entry => { + const keys = Object.keys(entry) as (keyof WeatherData)[]; + const newEntry: any = {}; + keys.forEach(key => newEntry[key.toLowerCase()] = entry[key]); + return newEntry as WeatherData; + }); + + // Order dataset by month and day + lowerCaseInput.sort((a, b) => (a.mon! - b.mon!) || (a.day! - b.day!)); + + let n = 1; + if (batch && lowerCaseInput.some(entry => entry.id !== undefined)) { + lowerCaseInput.sort((a, b) => (a.mon! - b.mon!) || (a.day! - b.day!) || (a.id! - b.id!)); + n = new Set(lowerCaseInput.map(entry => entry.id)).size; + const ids = lowerCaseInput.slice(0, n).map(entry => entry.id); + if (ids.length !== n) { + throw new Error("Multiple stations have to start and end at the same dates, and input data must be sorted by date/time and id"); + } + } else { + n = lowerCaseInput.length; + } + + const temp = lowerCaseInput.map(entry => entry.temp); + const prec = lowerCaseInput.map(entry => entry.prec); + const ws = lowerCaseInput.map(entry => entry.ws); + const rh = lowerCaseInput.map(entry => entry.rh); + const mon = lowerCaseInput.map(entry => entry.mon!); + const dmc = lowerCaseInput.map(entry => entry.dmc); + + if (!dmc.length) { + console.warn("dmc is missing!"); + } + if (!temp.length) { + console.warn("temperature (temp) is missing!"); + } + if (!prec.length) { + console.warn("precipitation (prec) is missing!"); + } + if (!ws.length) { + console.warn("wind speed (ws) is missing!"); + } + if (!rh.length) { + console.warn("relative humidity (rh) is missing!"); + } + if (temp.length % n !== 0) { + console.warn("Input data do not match with number of weather stations"); + } + + const el = [6.5, 7.5, 9.0, 12.8, 13.9, 13.9, 12.4, 10.9, 9.4, 8.0, 7.0, 6.0]; + + const constrainedRh = rh.map(val => Math.min(Math.max(val, 0), 99.9)); + const constrainedWs = ws.map(val => Math.max(val, 0)); + const constrainedPrec = prec.map(val => Math.max(val, 0)); + + const n0 = Math.floor(temp.length / n); + let SDMC: number[] = []; + + for (let i = 0; i < n0; i++) { + const k = Array.from({ length: n }, (_, idx) => n * i + idx); + + if (sdmc_old === null) { + sdmc_old = k.map(idx => { + let initial_sdmc = 2.6 + (1.7 * dmc[idx]) - 6.0; + return Math.max(initial_sdmc, 12); + }); + } + + const t0 = k.map(idx => Math.max(temp[idx], -1.1)); + const rk = k.map((idx, j) => 4.91 / 3.57 * 1.894 * (t0[j] + 1.1) * (100 - constrainedRh[idx]) * el[mon[idx] - 1] * 0.0001); + const rw = k.map(idx => constrainedPrec[idx] < 7.69 ? 0.218 * constrainedPrec[idx] - 0.094 : 0.83 * constrainedPrec[idx] - 4.8); + const wmi = k.map((_, j) => 20.0 + 280.0 / Math.exp(0.023 * sdmc_old![j])); + const b = k.map((_, j) => { + let b_val = sdmc_old![j] <= 33 ? 100.0 / (0.5 + 0.3 * sdmc_old![j]) : 14.0 - 1.3 * Math.log(sdmc_old![j]); + if (sdmc_old![j] > 65) { + b_val = 6.2 * Math.log(sdmc_old![j]) - 17.2; + } + return b_val; + }); + const wmr = k.map((_, j) => wmi[j] + 1000.0 * rw[j] / (48.77 + b[j] * rw[j])); + const pr = k.map((idx, j) => constrainedPrec[idx] <= 0.44 ? sdmc_old![j] : 43.43 * (5.6348 - Math.log(wmr[j] - 20))); + const constrainedPr = pr.map(val => Math.max(val, 0)); + const SDMC0 = k.map((_, j) => Math.max(constrainedPr[j] + rk[j], 0)); + + SDMC = SDMC.concat(SDMC0); + sdmc_old = SDMC0; + } + + return SDMC; + } + \ No newline at end of file diff --git a/src/fwi/surface_fuel_consumption.ts b/src/fwi/surface_fuel_consumption.ts new file mode 100644 index 0000000..be74209 --- /dev/null +++ b/src/fwi/surface_fuel_consumption.ts @@ -0,0 +1,89 @@ +/** + * Surface Fuel Consumption Calculator + * + * @description Computes the Surface Fuel Consumption by Fuel Type. + * All variable names are laid out in the same manner as FCFDG (1992) or + * Wotton et. al (2009). + * + * Forestry Canada Fire Danger Group (FCFDG) (1992). "Development and + * Structure of the Canadian Forest Fire Behavior Prediction System." + * Technical Report ST-X-3, Forestry Canada, Ottawa, Ontario. + * + * Wotton, B.M., Alexander, M.E., Taylor, S.W. 2009. Updates and revisions to + * the 1992 Canadian forest fire behavior prediction system. Nat. Resour. + * Can., Can. For. Serv., Great Lakes For. Cent., Sault Ste. Marie, Ontario, + * Canada. Information Report GLC-X-10, 45p. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param FFMC Fine Fuel Moisture Code + * @param BUI Buildup Index + * @param PC Percent Conifer (%) + * @param GFL Grass Fuel Load (kg/m^2) + * + * @returns SFC Surface Fuel Consumption (kg/m^2) + */ + +function surfaceFuelConsumption(FUELTYPE: string[], FFMC: number[], BUI: number[], PC: number[], GFL: number[]): number[] { + const SFC: number[] = Array(FFMC.length).fill(-999); + + for (let i = 0; i < FFMC.length; i++) { + switch (FUELTYPE[i]) { + case "C1": + SFC[i] = FFMC[i] > 84 ? + 0.75 + 0.75 * Math.sqrt(1 - Math.exp(-0.23 * (FFMC[i] - 84))) : + 0.75 - 0.75 * Math.sqrt(1 - Math.exp(-0.23 * (84 - FFMC[i]))); + break; + case "C2": + case "M3": + case "M4": + SFC[i] = 5.0 * (1 - Math.exp(-0.0115 * BUI[i])); + break; + case "C3": + case "C4": + SFC[i] = Math.pow(5.0 * (1 - Math.exp(-0.0164 * BUI[i])), 2.24); + break; + case "C5": + case "C6": + SFC[i] = Math.pow(5.0 * (1 - Math.exp(-0.0149 * BUI[i])), 2.48); + break; + case "C7": + SFC[i] = (FFMC[i] > 70 ? 2 * (1 - Math.exp(-0.104 * (FFMC[i] - 70))) : 0) + + 1.5 * (1 - Math.exp(-0.0201 * BUI[i])); + break; + case "D1": + SFC[i] = 1.5 * (1 - Math.exp(-0.0183 * BUI[i])); + break; + case "M1": + case "M2": + SFC[i] = (PC[i] / 100 * 5.0 * (1 - Math.exp(-0.0115 * BUI[i]))) + + ((100 - PC[i]) / 100 * 1.5 * (1 - Math.exp(-0.0183 * BUI[i]))); + break; + case "O1A": + case "O1B": + SFC[i] = GFL[i]; + break; + case "S1": + SFC[i] = 4.0 * (1 - Math.exp(-0.025 * BUI[i])) + 4.0 * (1 - Math.exp(-0.034 * BUI[i])); + break; + case "S2": + SFC[i] = 10.0 * (1 - Math.exp(-0.013 * BUI[i])) + 6.0 * (1 - Math.exp(-0.060 * BUI[i])); + break; + case "S3": + SFC[i] = 12.0 * (1 - Math.exp(-0.0166 * BUI[i])) + 20.0 * (1 - Math.exp(-0.0210 * BUI[i])); + break; + default: + SFC[i] = -999; + } + } + + // Constrain SFC values to be greater than 0 + return SFC.map(value => Math.max(value, 0.000001)); + } + + function deprecatedSurfaceFuelConsumption(FUELTYPE: string[], FFMC: number[], BUI: number[], PC: number[], GFL: number[]) { + console.warn("surfaceFuelConsumption is deprecated. Use surfaceFuelConsumption instead."); + return surfaceFuelConsumption(FUELTYPE, FFMC, BUI, PC, GFL); + } + + export { surfaceFuelConsumption }; + \ No newline at end of file diff --git a/src/fwi/total_fuel_consumption.ts b/src/fwi/total_fuel_consumption.ts new file mode 100644 index 0000000..0b2bf0a --- /dev/null +++ b/src/fwi/total_fuel_consumption.ts @@ -0,0 +1,58 @@ +/** + * Total Fuel Consumption calculation + * + * @description Computes the Total (Surface + Crown) Fuel Consumption by Fuel + * Type. + * All variable names are laid out in the same manner as FCFDG (1992) or + * Wotton et. al (2009). + * + * Forestry Canada Fire Danger Group (FCFDG) (1992). "Development and + * Structure of the Canadian Forest Fire Behavior Prediction System." + * Technical Report ST-X-3, Forestry Canada, Ottawa, Ontario. + * + * Wotton, B.M., Alexander, M.E., Taylor, S.W. 2009. Updates and revisions to + * the 1992 Canadian forest fire behavior prediction system. Nat. Resour. + * Can., Can. For. Serv., Great Lakes For. Cent., Sault Ste. Marie, Ontario, + * Canada. Information Report GLC-X-10, 45p. + * + * @param FUELTYPE The Fire Behaviour Prediction FuelType + * @param CFL Crown Fuel Load (kg/m^2) + * @param CFB Crown Fraction Burned (0-1) + * @param SFC Surface Fuel Consumption (kg/m^2) + * @param PC Percent Conifer (%) + * @param PDF Percent Dead Balsam Fir (%) + * @param option Type of output (TFC, CFC, default="TFC") + * + * @returns TFC Total (Surface + Crown) Fuel Consumption (kg/m^2) OR + * CFC Crown Fuel Consumption (kg/m^2) + */ + +type FuelType = "C1" | "C2" | "C3" | "C4" | "C5" | "C6" | "C7" | "D1" | "M1" | "M2" | "M3" | "M4" | "O1A" | "O1B" | "S1" | "S2" | "S3"; +type Option = "TFC" | "CFC"; + +function crownFuelConsumption(FUELTYPE: FuelType[], CFL: number[], CFB: number[], PC: number[], PDF: number[]): number[] { + return CFL.map((cfl, index) => { + let cfc = cfl * CFB[index]; + if (["M1", "M2"].includes(FUELTYPE[index])) { + cfc = (PC[index] / 100) * cfc; + } else if (["M3", "M4"].includes(FUELTYPE[index])) { + cfc = (PDF[index] / 100) * cfc; + } + return cfc; + }); +} + +function totalFuelConsumption(FUELTYPE: FuelType[], CFL: number[], CFB: number[], SFC: number[], PC: number[], PDF: number[], option: Option = "TFC"): number[] { + const CFC = crownFuelConsumption(FUELTYPE, CFL, CFB, PC, PDF); + + if (option === "CFC") { + return CFC; + } + + return SFC.map((sfc, index) => sfc + CFC[index]); +} + +function deprecatedTotalFuelConsumption(FUELTYPE: FuelType[], CFL: number[], CFB: number[], SFC: number[], PC: number[], PDF: number[], option: Option = "TFC") { + console.warn("totalFuelConsumption is deprecated. Use totalFuelConsumption instead."); + return totalFuelConsumption(FUELTYPE, CFL, CFB, SFC, PC, PDF, option); +}