# (C) Crown Copyright, Met Office. All rights reserved.
#
# This file is part of 'IMPROVER' and is released under the BSD 3-Clause license.
# See LICENSE in the root of the repository for full licensing details.
import os
from typing import Union
import numpy as np
from iris.cube import Cube, CubeList
from improver.fire_weather import IterativeFireWeatherBase
FFMC_START_VALUE = os.environ.get("FFMC_START_VALUE", 85)
FFMC_LAG_TIME = os.environ.get("FFMC_LAG_TIME", 3)
[docs]
class FineFuelMoistureCode(IterativeFireWeatherBase):
"""
Plugin to calculate the Fine Fuel Moisture Code (FFMC).
The FFMC is a numerical rating of the moisture content of litter and other
fine fuels, representing the relative ease of ignition and flammability of fine fuel.
Values range from 0-101, with higher values indicating drier conditions.
This process is adapted directly from:
Equations and FORTRAN Program for the
Canadian Forest Fire Weather Index System
(C.E. Van Wagner and T.L. Pickett, 1985).
Page 5, Equations 1-10.
Expected input units:
- Temperature: degrees Celsius
- Precipitation: mm (24-hour accumulation)
- Relative humidity: percentage (0-100)
- Wind speed: km/h
- Previous FFMC: dimensionless (0-101)
"""
STARTING_VALUE = FFMC_START_VALUE
LAG_TIME = FFMC_LAG_TIME
METADATA_SOURCE_CUBE = "fine_fuel_moisture_code"
INPUT_CUBE_NAMES = [
"air_temperature",
"lwe_thickness_of_precipitation_amount",
"relative_humidity",
"wind_speed",
METADATA_SOURCE_CUBE,
]
OUTPUT_CUBE_NAME = "fine_fuel_moisture_code"
# Valid output ranges for warning checks (output_name: (min, max))
# Minimum and maximum feasible values for each output index are drawn from
# values reported in:
# Wang, X., Oliver, J., Swystun, T., Hanes, C.C., Erni, S. and Flannigan,
# M.D., 2023. Critical fire weather conditions during active fire spread
# days in Canada. Science of the total environment, 869, p.161831.
VALID_OUTPUT_RANGE = (0.0, 101.0)
# Disambiguate input FFMC (yesterday's value) from output FFMC (today's calculated value)
INPUT_ATTRIBUTE_MAPPINGS = {"fine_fuel_moisture_code": "input_ffmc"}
temperature: Cube
precipitation: Cube
relative_humidity: Cube
wind_speed: Cube
input_ffmc: Cube
initial_moisture_content: np.ndarray
moisture_content: np.ndarray
[docs]
def process(
self,
*cubes: Union[Cube, CubeList],
month: int | None = None,
initialise: bool = False,
clip_ffmc: bool = False,
) -> Cube:
"""
Args:
*cubes:
One or more input cubes as specified by INPUT_CUBE_NAMES. When initialise is True `cubes` should
exclude the OUTPUT_CUBE_NAME, which should otherwise be given as the iterative input.
month:
Month parameter (1-12), required only if REQUIRES_MONTH is True
initialise:
True when starting the iterative process else False
clip_ffmc:
If true fine fuel moisture code values will be clipped to
a minimum of 0 and a maximum of 101.
Returns:
The calculated output cube.
"""
self.clip_ffmc = clip_ffmc
return super().process(cubes, month=month, initialise=initialise)
[docs]
def _calculate(self) -> np.ndarray:
"""Calculate the Fine Fuel Moisture Code (FFMC).
Returns:
The calculated FFMC values for the current day.
"""
# Step 1 & 2: Calculate the previous day's moisture content
self._calculate_moisture_content()
# Step 3: Perform rainfall adjustment
self._perform_rainfall_adjustment()
# Step 4: Calculate Equilibrium Moisture Content for the drying phase
E_d = self._calculate_EMC_for_drying_phase()
# VECTORIZATION NOTE: The original algorithm in Van Wagner and Pickett (1985)
# is applied to each point linearly. We are applying this to Cubes. This means
# that the order of Steps 5-7 is not identical to the original algorithm.
# Step 5a & 5b: Calculate moisture content through drying rate
# VECTORIZATION NOTE: We calculate drying for all points, then apply selectively.
# This differs from the scalar algorithm which only calculates when moisture content > E_d.
moisture_content_from_drying = (
self._calculate_moisture_content_through_drying_rate(E_d)
)
mask_wetting = self.moisture_content > E_d
# Step 6: Calculate Equilibrium Moisture Content for the wetting phase
# VECTORIZATION NOTE: We calculate E_w for all points for efficiency,
# though it's only needed where the moisture content < E_d
E_w = self._calculate_EMC_for_wetting_phase()
# Step 7a & 7b: Calculate moisture content through wetting equilibrium
# VECTORIZATION NOTE: We calculate wetting for all points, then apply selectively.
# This differs from the scalar algorithm which only calculates when moisture content < E_d
# AND moisture content < E_w.
moisture_content_from_wetting = (
self._calculate_moisture_content_through_wetting_equilibrium(E_w)
)
mask_drying = (
self.moisture_content < E_d
) # Identifies where wetting might apply
mask_apply_wetting = np.logical_and(mask_drying, self.moisture_content < E_w)
# Step 5b, 7b, 8: Apply the appropriate transformation using mutually exclusive logic.
# VECTORIZATION NOTE: We use nested np.where to implement the three-way
# if-elseif-else structure from the scalar algorithm, ensuring only one
# transformation is applied per grid point.
self.moisture_content = np.where(
mask_wetting,
# If moisture content > E_d: apply drying (Step 5b, Equation 8)
moisture_content_from_drying,
np.where(
mask_apply_wetting,
# Else if (moisture_content < E_d) AND (moisture_content < E_w): apply wetting (Step 7b, Equation 9)
moisture_content_from_wetting,
# Else: no change, (Step 8: E_d >= moisture_content >= E_w)
self.moisture_content,
),
)
# Step 9: Calculate fine fuel moisture code (FFMC) from moisture content
ffmc = self._calculate_ffmc_from_moisture_content(self.clip_ffmc)
return ffmc
[docs]
def _calculate_moisture_content(self):
"""Calculates the previous day's moisture content for a given input value
of the fine fuel moisture code, and initialises the moisture_content
attribute to that value.
From Van Wagner and Pickett (1985), Page 5: Equation 1, Steps 1 & 2.
"""
# Steps 1 & 2: Calculate the previous day's moisture content
self.initial_moisture_content = (
147.2 * (101.0 - self.input_ffmc.data) / (59.5 + self.input_ffmc.data)
)
# Initialise today's moisture content to the previous day's value
self.moisture_content = self.initial_moisture_content.copy()
[docs]
def _calculate_EMC_for_drying_phase(self) -> np.ndarray:
"""Calculates the Equilibrium Moisture Content (EMC) for the drying phase
under current environmental conditions (relative humidity, and temperature)
From Van Wagner and Pickett (1985), Page 5: Equation 4, and Step 4.
Returns:
The Equilibrium Moisture Content for the drying phase (E_d). Array
shape matches the input cube data shape. Values are in moisture
content units (dimensionless).
"""
# Equation 4: Calculate EMC for drying phase (E_d)
E_d = (
0.942 * self.relative_humidity.data**0.679
+ 11 * np.exp((self.relative_humidity.data - 100) / 10)
+ 0.18
* (21.1 - self.temperature.data)
* (1 - np.exp(-0.115 * self.relative_humidity.data))
)
return E_d
[docs]
def _calculate_moisture_content_through_drying_rate(
self,
E_d: np.ndarray,
) -> np.ndarray:
"""Calculates the moisture content through the drying rate.
From Van Wagner and Pickett (1985), Page 5: Equations 6a, 6b, 8, and Step 5.
Args:
E_d:
The Equilibrium Moisture Content for the drying phase.
Returns:
Array of moisture content (dimensionless) with drying applied at all
grid points. Shape matches input cube data shape.
"""
# Equation 6a: Calculate the log drying rate intermediate step
k_o = 0.424 * (
1 - (self.relative_humidity.data / 100.0) ** 1.7
) + 0.0694 * np.sqrt(self.wind_speed.data) * (
1 - (self.relative_humidity.data / 100.0) ** 8
)
# Equation 6b: Calculate the log drying rate
k_d = k_o * 0.581 * np.exp(0.0365 * self.temperature.data)
# Equation 8: Calculate the new moisture content via drying
new_moisture_content = E_d + (self.moisture_content - E_d) * 10 ** (-k_d)
return new_moisture_content
[docs]
def _calculate_EMC_for_wetting_phase(self) -> np.ndarray:
"""Calculates the Equilibrium Moisture Content (EMC) for the wetting phase
under current environmental conditions (relative humidity, and temperature)
From Van Wagner and Pickett (1985), Page 5: Equation 5, and Step 6.
Returns:
The Equilibrium Moisture Content for the wetting phase (E_w). Array
shape matches the input cube data shape. Values are in moisture
content units (dimensionless).
"""
# Equation 5: Calculate the EMC for the wetting phase (E_w)
E_w = (
0.618 * self.relative_humidity.data**0.753
+ 10.0 * np.exp((self.relative_humidity.data - 100.0) / 10.0)
+ 0.18
* (21.1 - self.temperature.data)
* (1 - np.exp(-0.115 * self.relative_humidity.data))
)
return E_w
[docs]
def _calculate_moisture_content_through_wetting_equilibrium(
self,
E_w: np.ndarray,
) -> np.ndarray:
"""Calculates the moisture content through the wetting equilibrium.
From Van Wagner and Pickett (1985), Page 5: Equations 7a, 7b, 9, and Step 7.
Args:
E_w:
The Equilibrium Moisture Content for the wetting phase.
Returns:
Array of moisture content (dimensionless) with wetting applied at all
grid points. Shape matches input cube data shape.
"""
# Equation 7a: Calculate the log wetting rate intermediate step
k_l = 0.424 * (
1 - ((100.0 - self.relative_humidity.data) / 100.0) ** 1.7
) + 0.0694 * np.sqrt(self.wind_speed.data) * (
1 - ((100.0 - self.relative_humidity.data) / 100.0) ** 8
)
# Equation 7b: Calculate the log wetting rate intermediate step
k_w = k_l * 0.581 * np.exp(0.0365 * self.temperature.data)
# Equation 9: Calculate the new moisture content via wetting
new_moisture_content = E_w - (E_w - self.moisture_content) * 10 ** (-k_w)
return new_moisture_content
[docs]
def _calculate_ffmc_from_moisture_content(
self, clip_ffmc: bool = False
) -> np.ndarray:
"""Calculates the fine fuel moisture code (FFMC) from the moisture
content.
From Van Wagner and Pickett (1985), Page 5: Equation 10, and Step 9.
Args:
clip_ffmc:
If true fine fuel moisture code values will be clipped to
a minimum of 0 and a maximum of 101.
Returns:
The calculated FFMC values (dimensionless, clipped range 0-101).
Return array will be clipped only if clip_ffmc is True.
Array shape matches input cube data shape.
"""
# Equation 10: Calculate FFMC from moisture content
ffmc = 59.5 * (250.0 - self.moisture_content) / (147.2 + self.moisture_content)
if clip_ffmc:
return np.clip(ffmc, 0, 101)
return ffmc
