EPPT-2585: Fire Severity Index: Add a DuffMoistureCode class to IMPROVER

improver/api/__init__.py

@@ -51,6 +51,7 @@
     "DayNightMask": "improver.utilities.solar",
     "DifferenceBetweenAdjacentGridSquares": "improver.utilities.spatial",
     "DroughtCode": "improver.fire_weather.drought_code",
+    "DuffMoistureCode": "improver.fire_weather.duff_moisture_code",
     "EnforceConsistentForecasts": "improver.utilities.forecast_reference_enforcement",
     "EnsembleReordering": "improver.ensemble_copula_coupling.ensemble_copula_coupling",
     "EstimateCoefficientsForEnsembleCalibration": "improver.calibration.emos_calibration",

improver/fire_weather/__init__.py

@@ -104,18 +104,20 @@ class FireWeatherIndexBase(BasePlugin):
     }
 
     def load_input_cubes(self, cubes: tuple[Cube] | CubeList, month: int | None = None):
-        """Loads the required input cubes for the calculation. These
-        are stored internally as Cube objects.
+        """Loads the required input cubes for the calculation. These are stored
+        internally as Cube objects.
 
         Args:
-            cubes (tuple[iris.cube.Cube] | iris.cube.CubeList): Input cubes containing the necessary data.
-            month (int | None): Month of the year (1-12), required only if REQUIRES_MONTH is True.
+            cubes:
+                Input cubes containing the necessary data.
+            month:
+                Month of the year (1-12), required only if REQUIRES_MONTH is True.
                 Defaults to None.
 
         Raises:
-            ValueError: If the number of cubes does not match the expected
-                number, if month is required but not provided, or if month
-                is out of range.
+            ValueError:
+                If the number of cubes does not match the expected number, if
+                month is required but not provided, or if month is out of range.
         """
         if len(cubes) != len(self.INPUT_CUBE_NAMES):
             raise ValueError(
@@ -152,10 +154,11 @@ def _get_attribute_name(self, standard_name: str) -> str:
         """Convert a cube standard name to an attribute name.
 
         Args:
-            standard_name (str): The cube's standard name
+            standard_name:
+                The cube's standard name
 
         Returns:
-            str: The attribute name to use for storing the cube
+            The attribute name to use for storing the cube
 
         Examples:
             "air_temperature" -> "temperature"
@@ -176,11 +179,14 @@ def _validate_input_range(self, cube: Cube, attr_name: str) -> None:
         """Validate that input data falls within expected physical ranges.
 
         Args:
-            cube (iris.cube.Cube): The input cube to validate
-            attr_name (str): The attribute name for the cube
+            cube:
+                The input cube to validate
+            attr_name:
+                The attribute name for the cube
 
         Raises:
-            ValueError: If any values fall outside the valid range for this input type,
+            ValueError:
+                If any values fall outside the valid range for this input type,
                 or if data contains NaN or Inf values
         """
         if attr_name not in self._VALID_RANGES:
@@ -223,13 +229,15 @@ def _make_output_cube(
         24-hour accumulation period.
 
         Args:
-            data (np.ndarray): The output data array
-            template_cube (iris.cube.Cube | None): The cube to use as a template for metadata.
+            data:
+                The output data array
+            template_cube:
+                The cube to use as a template for metadata.
                 If None, uses the first input cube. Defaults to None
 
         Returns:
-            iris.cube.Cube: The output cube containing the output data with proper
-                metadata and coordinates.
+            The output cube containing the output data with proper metadata
+            and coordinates.
         """
         if template_cube is None:
             # Use first input cube as template
@@ -267,23 +275,27 @@ def _calculate(self) -> np.ndarray:
         the specific calculation logic for that component.
 
         Raises:
-            NotImplementedError: This method must be implemented by subclasses.
+            NotImplementedError:
+                This method must be implemented by subclasses.
         """
         raise NotImplementedError("Subclasses must implement the _calculate method.")
 
     def process(self, cubes: tuple[Cube] | CubeList, month: int | None = None) -> Cube:
         """Calculate the fire weather index component.
 
         Args:
-            cubes (tuple[iris.cube.Cube] | iris.cube.CubeList): Input cubes as specified by INPUT_CUBE_NAMES
-            month (int | None): Month parameter (1-12), required only if REQUIRES_MONTH is True
+            cubes:
+                Input cubes as specified by INPUT_CUBE_NAMES
+            month:
+                Month parameter (1-12), required only if REQUIRES_MONTH is True
                 Defaults to None.
 
         Returns:
-            iris.cube.Cube: The calculated output cube.
+            The calculated output cube.
 
         Warns:
-            UserWarning: If output values fall outside typical expected ranges
+            UserWarning:
+                If output values fall outside typical expected ranges
         """
         self.load_input_cubes(cubes, month)
         output_data = self._calculate()
@@ -298,10 +310,12 @@ def _validate_output_range(self, output_cube: Cube) -> None:
         """Check if output values fall within expected ranges and issue warnings if not.
 
         Args:
-            output_cube (iris.cube.Cube): The output cube to validate
+            output_cube:
+                The output cube to validate
 
         Warns:
-            UserWarning: If output contains NaN, Inf, or values outside expected ranges
+            UserWarning:
+                If output contains NaN, Inf, or values outside expected ranges
         """
         output_name = output_cube.name()
 

improver/fire_weather/duff_moisture_code.py

@@ -0,0 +1,197 @@
+# (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 numpy as np
+from iris.cube import Cube
+
+from improver.fire_weather import FireWeatherIndexBase
+
+
+class DuffMoistureCode(FireWeatherIndexBase):
+    """
+    Plugin to calculate the Duff Moisture Code (DMC) following
+    the Canadian Forest Fire Weather Index System.
+
+    The DMC is a numerical rating of the average moisture content of loosely
+    compacted organic layers of moderate depth. It indicates fuel consumption
+    in moderate duff layers and medium-size woody material.
+    Higher values indicate 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 6, Equations 11-16.
+
+    Expected input units:
+        - Temperature: degrees Celsius
+        - Precipitation: mm (24-hour accumulation)
+        - Relative humidity: percentage (0-100)
+        - Previous DMC: dimensionless
+        - Month: integer (1-12) for day length factor lookup
+    """
+
+    INPUT_CUBE_NAMES = [
+        "air_temperature",
+        "lwe_thickness_of_precipitation_amount",
+        "relative_humidity",
+        "duff_moisture_code",
+    ]
+    OUTPUT_CUBE_NAME = "duff_moisture_code"
+    REQUIRES_MONTH = True
+    # Disambiguate input DMC (yesterday's value) from output DMC (today's calculated value)
+    INPUT_ATTRIBUTE_MAPPINGS = {"duff_moisture_code": "input_dmc"}
+
+    temperature: Cube
+    precipitation: Cube
+    relative_humidity: Cube
+    input_dmc: Cube
+    month: int
+    previous_dmc: np.ndarray
+
+    # Day length factors for DMC calculation (L_e values from Table 3)
+    # Index 0 is unused, indices 1-12 correspond to months January-December
+    DMC_DAY_LENGTH_FACTORS = [
+        0.0,  # Placeholder for index 0
+        6.5,  # January
+        7.5,  # February
+        9.0,  # March
+        12.8,  # April
+        13.9,  # May
+        13.9,  # June
+        12.4,  # July
+        10.9,  # August
+        9.4,  # September
+        8.0,  # October
+        7.0,  # November
+        6.0,  # December
+    ]
+
+    def _calculate(self) -> np.ndarray:
+        """Calculate the Duff Moisture Code (DMC).
+
+        Returns:
+            np.ndarray: The calculated DMC values for the current day.
+        """
+        # Step 1: Set today's DMC value to the previous day's DMC value
+        self.previous_dmc = self.input_dmc.data.copy()
+
+        # Step 2: Perform rainfall adjustment, if precipitation > 1.5 mm
+        self._perform_rainfall_adjustment()
+
+        # Steps 3 & 4: Calculate drying rate
+        drying_rate = self._calculate_drying_rate()
+
+        # Step 5: Calculate DMC from adjusted previous DMC and drying rate
+        dmc = self._calculate_dmc(drying_rate)
+
+        return dmc
+
+    def _perform_rainfall_adjustment(self):
+        """Updates the previous DMC value based on available precipitation
+        accumulation data for the previous 24 hours. This is done element-wise
+        for each grid point.
+
+        From Van Wagner and Pickett (1985), Page 6: Equations 11-15,
+        and Steps 2a-2e corresponding to rainfall adjustment.
+        """
+        # Only adjust if precipitation > 1.5 mm
+        precip_mask = self.precipitation.data > 1.5
+
+        # Step 2a: Calculate effective rainfall via Equation 11
+        effective_rain = 0.92 * self.precipitation.data - 1.27
+
+        # Step 2b: Calculate initial moisture content from previous DMC via Equation 12
+        moisture_content_initial = 20.0 + np.exp(5.6348 - self.previous_dmc / 43.43)
+
+        # Step 2c: Calculate the slope variable based on previous DMC via Equations 13a, 13b, 13c
+        # In the original algorithm, the slope variable is referred to as 'b'
+        # VECTORIZATION NOTE: This structure matches the original algorithm while being vectorized
+        # Clip previous_dmc to avoid log(0) warnings in equations 13b and 13c
+        dmc_clipped = np.maximum(self.previous_dmc, 1e-10)
+        slope_variable = np.where(
+            self.previous_dmc <= 33.0,
+            # Equation 13a: DMC <= 33
+            100.0 / (0.5 + 0.3 * self.previous_dmc),
+            np.where(
+                self.previous_dmc <= 65.0,
+                # Equation 13b: 33 < DMC <= 65
+                14.0 - 1.3 * np.log(dmc_clipped),
+                # Equation 13c: DMC > 65
+                6.2 * np.log(dmc_clipped) - 17.2,
+            ),
+        )
+
+        # Step 2d: Calculate moisture content after rain via Equation 14
+        # Protect against division by zero (though mathematically unlikely)
+        denominator = 48.77 + slope_variable * effective_rain
+        moisture_content_after_rain = moisture_content_initial + (
+            1000.0 * effective_rain
+        ) / np.maximum(denominator, 1e-10)
+
+        # Step 2e: Calculate DMC after rain via Equation 15
+        # This is modified to avoid log of zero or negative values
+        log_arg = np.clip(moisture_content_after_rain - 20.0, 1e-10, None)
+        dmc_after_rain = 244.72 - 43.43 * np.log(log_arg)
+
+        # Apply lower bound of 0
+        dmc_after_rain = np.maximum(dmc_after_rain, 0.0)
+
+        # Update previous_dmc where precipitation > 1.5
+        self.previous_dmc = np.where(precip_mask, dmc_after_rain, self.previous_dmc)
+
+    def _calculate_drying_rate(self) -> np.ndarray:
+        """Calculates the drying rate for DMC. This is multiplied by 100 for
+        computational efficiency in the final DMC calculation. The original
+        algorithm calculates K and then multiplies it by 100 in the DMC equation.
+
+        Temperature is bounded to a minimum of -1.1°C. Uses the day length factor
+        for the current month from DMC_DAY_LENGTH_FACTORS.
+
+        From Van Wagner and Pickett (1985), Page 6: Equation 16, Steps 3 & 4.
+
+        Returns:
+            The drying rate (dimensionless). Shape matches input cube data
+            shape. This value is added to previous DMC to get today's DMC.
+        """
+        # Apply temperature lower bound of -1.1°C
+        temp_adjusted = np.maximum(self.temperature.data, -1.1)
+
+        # Step 3: Get day length factor for current month
+        day_length_factor = self.DMC_DAY_LENGTH_FACTORS[self.month]
+
+        # Step 4: Calculate drying rate via Equation 16
+        drying_rate = (
+            1.894
+            * (temp_adjusted + 1.1)
+            * (100.0 - self.relative_humidity.data)
+            * day_length_factor
+            * 1e-4
+        )
+
+        return drying_rate
+
+    def _calculate_dmc(self, drying_rate: np.ndarray) -> np.ndarray:
+        """Calculates the Duff Moisture Code from previous DMC and drying rate.
+        Note that the drying rate is expected to be pre-multiplied by 100
+        for computational efficiency. This mathematically matches the original
+        algorithm, but is more efficient to implement this way.
+
+        From Van Wagner and Pickett (1985), Page 6: Equation 16.
+
+        Args:
+            drying_rate:
+                The drying rate, represented by 'K' in the original equations.
+
+        Returns:
+            The calculated DMC value.
+        """
+        # Equation 16: Calculate DMC
+        dmc = self.previous_dmc + drying_rate
+
+        # Apply lower bound of 0
+        dmc = np.maximum(dmc, 0.0)
+
+        return dmc