Source code for improver.fire_weather.duff_moisture_code

# (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

import numpy as np
from iris.cube import Cube

from improver.fire_weather import IterativeFireWeatherBase

DMC_START_VALUE = os.environ.get("DMC_START_VALUE", 6)
DMC_LAG_TIME = os.environ.get("DMC_LAG_TIME", 15)




class DuffMoistureCode(IterativeFireWeatherBase):
    """
    Plugin to calculate the Duff Moisture Code (DMC).

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

    STARTING_VALUE = DMC_START_VALUE
    LAG_TIME = DMC_LAG_TIME

    METADATA_SOURCE_CUBE = "duff_moisture_code"
    INPUT_CUBE_NAMES = [
        "air_temperature",
        "lwe_thickness_of_precipitation_amount",
        "relative_humidity",
        METADATA_SOURCE_CUBE,
    ]
    OUTPUT_CUBE_NAME = "duff_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, 400)
    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