[+] Add starter files

2021-10-22 20:36:12 -04:00
parent e128848a4e
commit 9ac593288d
25 changed files with 9051 additions and 0 deletions
@@ -0,0 +1,291 @@
+"""CSC110 Fall 2021 Assignment 3: Forest Fire Weather Index System
+
+Copyright and Usage Information
+===============================
+
+This file is provided solely for the personal and private use of students
+taking CSC110 at the University of Toronto St. George campus. All forms of
+distribution of this code, whether as given or with any changes, are
+expressly prohibited. For more information on copyright for CSC110 materials,
+please consult our Course Syllabus.
+
+This file is Copyright (c) 2021 Mario Badr and Tom Fairgrieve.
+"""
+import math
+from dataclasses import dataclass
+
+# Initial values that can be used for FFMC, DMC, and DC
+INITIAL_FFMC = 85.0
+INITIAL_DMC = 6.0
+INITIAL_DC = 15.0
+
+# Per-month lookup tables
+DMC_DAY_LENGTH_EFFECTIVE = {
+    1: 6.5, 2: 7.5, 3: 9.0, 4: 12.8, 5: 13.9, 6: 13.9,
+    7: 12.4, 8: 10.9, 9: 9.4, 10: 8.0, 11: 7.0, 12: 6.0
+}
+DC_DAY_LENGTH_FACTORS = {
+    1: -1.6, 2: -1.6, 3: -1.6, 4: 0.9, 5: 3.8, 6: 5.8,
+    7: 6.4, 8: 5.0, 9: 2.4, 10: 0.4, 11: -1.6, 12: -1.6
+}
+
+
+@dataclass
+class WeatherMetrics:
+    """A bundle of metrics that are measured by weather stations.
+
+    Instance Attributes:
+        - month: The month of the year (e.g., January is 1, December is 12)
+        - day: the day of the month
+        - temperature: The noon temperature in degrees Celsius
+        - humidity: The noon relative humidity, in %
+        - wind_speed: The noon wind speed, in km/h
+        - precipitation: The rainfall at noon, in mm
+
+    Representation Invariants:
+        - self.month in {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}
+        - self.humidity >= 0.0
+        - self.wind_speed >= 0.0
+        - self.precipitation >= 0.0
+    """
+    month: int
+    day: int
+    temperature: float
+    humidity: float
+    wind_speed: float
+    precipitation: float
+
+
+@dataclass
+class FfwiOutput:
+    """A bundle of the output metrics in the Canadian Forest Fire Weather Index System.
+
+    Instance Attributes:
+        - ffmc: the Fine Fuel Moisture Code
+        - dmc: the Duff Moisture Code
+        - dc: the Drought Code
+        - isi: the Initial Spread Index
+        - bui: the Buildup Index
+        - fwi: the Fire Weather Index
+
+    Representation Invariants:
+        - 0.0 <= self.ffmc <= 101.0
+        - self.dmc >= 1.0
+        - self.bui >= 0.0
+    """
+    ffmc: float
+    dmc: float
+    dc: float
+    isi: float
+    bui: float
+    fwi: float
+
+
+def calculate_mr(precipitation: float, mo: float) -> float:
+    """Return the fine fuel moisture content after rain (mr) based on the effective
+    rainfall in precipitation and the fine fuel moisture content from the previous day.
+
+    Preconditions:
+        - precipitation > 0.5
+    """
+    # Equation 2
+    rf = precipitation - 0.5
+
+    if mo > 150.0:
+        # Equation 3b
+        mr = (mo + 42.5 * rf * math.exp(-100.0 / (251.0 - mo)) * (1.0 - math.exp(-6.93 / rf))) + (
+            0.0015 * (mo - 150.0) ** 2) * math.sqrt(rf)
+    else:
+        # Equation 3a
+        mr = mo + 42.5 * rf * math.exp(-100.0 / (251.0 - mo)) * (1.0 - math.exp(-6.93 / rf))
+
+    return mr
+
+
+def calculate_m(wm: WeatherMetrics, ed: float, mo: float) -> float:
+    """Return the fine fuel moisture content after drying (m) based on the measurements in wm, the
+    EMC for drying in ed, and the fine fuel moisture content from the previous day in mo.
+
+    Preconditions:
+        - mo <= 250.0
+    """
+    if mo == ed:
+        return mo
+
+    if mo < ed:
+        # Equation 5
+        ew = 0.618 * (wm.humidity ** .753) + (10.0 * math.exp((wm.humidity - 100.0) / 10.0)) \
+             + 0.18 * (21.1 - wm.temperature) * (1.0 - 1.0 / math.exp(0.115 * wm.humidity))
+
+        if mo <= ew:
+            # Use log wetting rate
+            # Equation 7a
+            k1 = 0.424 * (1.0 - ((100.0 - wm.humidity) / 100.0) ** 1.7) + \
+                 (.0694 * math.sqrt(wm.wind_speed)) * (1.0 - ((100.0 - wm.humidity) / 100.0) ** 8)
+            kw = k1 * (0.581 * math.exp(0.0365 * wm.temperature))  # Equation 7b
+            return ew - (ew - mo) / 10.0 ** kw  # Equation 9
+        else:
+            return mo
+
+    # Use log drying rate
+    # Equation 6a
+    k0 = 0.424 * (1.0 - (wm.humidity / 100.0) ** 1.7) + (
+            (0.0694 * math.sqrt(wm.wind_speed)) * (1.0 - (wm.humidity / 100.0) ** 8))
+    kd = k0 * (0.581 * math.exp(0.0365 * wm.temperature))  # Equation 6b
+    m = ed + (mo - ed) / 10.0 ** kd  # Equation 8
+
+    return m
+
+
+def calculate_ffmc(wm: WeatherMetrics, f0: float) -> float:
+    """Return the Fine Fuel Moisture Code (FFMC) based on the measurements in wm and the previous
+    day's FFMC in f0.
+    """
+    # Calculate the fine fuel moisture content from the previous day
+    mo = (147.2 * (101.0 - f0)) / (59.5 + f0)  # Equation 1
+    if wm.precipitation > 0.5:
+        mo = calculate_mr(wm.precipitation, mo)
+
+    if mo > 250.0:
+        mo = 250.0
+
+    # Equation 4 - Fine Fuel equilibrium moisture content (EMC) for drying
+    ed = 0.942 * (wm.humidity ** .679) + (11.0 * math.exp((wm.humidity - 100.0) / 10.0)) + (
+            0.18 * (21.1 - wm.temperature) * (1.0 - 1.0 / math.exp(0.1150 * wm.humidity)))
+
+    m = calculate_m(wm, ed, mo)
+
+    # Equation 10
+    f = (59.5 * (250.0 - m)) / (147.2 + m)
+
+    if f > 101.0:
+        f = 101.0
+    elif f <= 0.0:
+        f = 0.0
+
+    return f
+
+
+def calculate_dmr(precipitation: float, dm0: float) -> float:
+    """Calculate the Duff moisture content after rain based on the current precipitation and the
+    previous day's DMC in dm0.
+
+    Preconditions:
+        - precipitation > 1.5
+    """
+    rw = 0.92 * precipitation - 1.27  # Equation 11
+    wmi = 20.0 + 280.0 / math.exp(0.023 * dm0)  # Equation 12
+
+    if dm0 <= 33.0:
+        b = 100.0 / (0.5 + 0.3 * dm0)  # Equation 13a
+    else:
+        if dm0 <= 65.0:
+            b = 14.0 - 1.3 * math.log(dm0)  # Equation 13b
+        else:
+            b = 6.2 * math.log(dm0) - 17.2  # Equation 13c
+
+    return wmi + (1000 * rw) / (48.77 + b * rw)  # Equation 14
+
+
+def calculate_dmc_k(temperature: float, humidity: float, month: int) -> float:
+    """Return the log drying rate in DMC based on the temperature, humidity, and month."""
+    if temperature < -1.1:
+        # Cannot use temperatures less than -1.1 in Equation 16
+        temperature = -1.1
+
+    # Equations 16 and 17
+    return 1.894 * (temperature + 1.1) * (100.0 - humidity) * (
+            DMC_DAY_LENGTH_EFFECTIVE[month] * 0.0001)
+
+
+def calculate_dmc(wm: WeatherMetrics, dm0: float) -> float:
+    """Return the Duff Moisture Code (DMC) based on the measurements in wm and the previous day's
+    DMC in dm0.
+    """
+    if wm.precipitation <= 1.5:
+        pr = dm0
+    else:
+        dmr = calculate_dmr(wm.precipitation, dm0)
+        pr = 43.43 * (5.6348 - math.log(dmr - 20.0))  # Equation 15
+
+    # The DMC after rain cannot, theoretically, be less than 0; ensure it is at least 0
+    pr = max(0.0, pr)
+
+    rk = calculate_dmc_k(wm.temperature, wm.humidity, wm.month)
+    d = pr + rk
+
+    if d <= 1.0:
+        d = 1.0
+    return d
+
+
+def calculate_qr(precipitation: float, dc0: float) -> float:
+    """Return the moisture equivalent after rain based on the current precipitation and the
+    previous day's DC in dc0.
+
+    Preconditions:
+        - precipitation > 2.8
+    """
+    rd = 0.83 * precipitation - 1.27  # Equation 18
+    qo = 800.0 * math.exp(-dc0 / 400.0)  # Equation 19
+    qr = qo + 3.937 * rd  # Equation 20
+
+    return qr
+
+
+def calculate_dc(wm: WeatherMetrics, dc0: float) -> float:
+    """Return the Drought Code (DC) based on the measurements in wm and the previous day's DC in
+    dc0.
+    """
+    temperature = wm.temperature
+    if wm.temperature < -2.8:
+        temperature = -2.8
+
+    v = 0.36 * (temperature + 2.8) + DC_DAY_LENGTH_FACTORS[wm.month]  # Equation 22
+    # The potential evapotranspiration cannot, theoretically, be less than 0; ensure it's at least 0
+    v = max(0.0, v)
+
+    if wm.precipitation > 2.8:
+        qr = calculate_qr(wm.precipitation, dc0)
+        dr = 400 * math.log(800 / qr)  # Equation 21
+        return dr + 0.5 * v  # Equation 23
+    else:
+        return dc0 + 0.5 * v  # Equation 23
+
+
+def calculate_isi(wm: WeatherMetrics, ffmc: float) -> float:
+    """Return the Initial Spread Index (ISI) based on the measurements in wm and the current ffmc.
+    """
+    mo = 147.2 * (101.0 - ffmc) / (59.5 + ffmc)
+    ff = 19.115 * math.exp(mo * -0.1386) * (1.0 + (mo ** 5.31) / 49300000.0)
+
+    return ff * math.exp(0.05039 * wm.wind_speed)
+
+
+def calculate_bui(dmc: float, dc: float) -> float:
+    """Return the Buildup Index (BUI) based on the current dmc and the current dc.
+    """
+    if dmc <= 0.4 * dc:
+        b = (0.8 * dc * dmc) / (dmc + 0.4 * dc)
+    else:
+        b = dmc - (1.0 - 0.8 * dc / (dmc + 0.4 * dc)) * (0.92 + (0.0114 * dmc) ** 1.7)
+
+    if b < 0.0:
+        b = 0.0
+
+    return b
+
+
+def calculate_fwi(isi: float, bui: float) -> float:
+    """Return the Fire Weather Index (FWI) based on the current isi and the current bui.
+    """
+    if bui <= 80.0:
+        f_d = 0.626 * bui ** 0.809 + 2.0  # Equation 28a
+    else:
+        f_d = 1000.0 / (25. + 108.64 * math.exp(-0.023 * bui))  # Equation 28b
+
+    bb = 0.1 * isi * f_d  # Equation 29
+    if bb <= 1.0:
+        return bb  # Equation 30b
+    else:
+        return math.exp(2.72 * (0.434 * math.log(bb)) ** 0.647)  # Equation 30a