SFMS: Daily FFMC

api/app/auto_spatial_advisory/sfms.py

@@ -1,6 +1,7 @@
-from cffdrs import bui, dc, dmc
+from cffdrs import bui, dc, dmc, ffmc
 from numba import vectorize
 
 vectorized_bui = vectorize(bui)
 vectorized_dc = vectorize(dc)
 vectorized_dmc = vectorize(dmc)
+vectorized_ffmc = vectorize(ffmc)

api/app/jobs/sfms_calculations.py

@@ -6,7 +6,7 @@
 
 from app import configure_logging
 from app.rocketchat_notifications import send_rocketchat_notification
-from app.sfms.date_range_processor import BUIDateRangeProcessor
+from app.sfms.daily_fwi_processor import DailyFWIProcessor
 from app.sfms.raster_addresser import RasterKeyAddresser
 from app.utils.s3_client import S3Client
 from app.utils.time import get_utc_now
@@ -19,7 +19,7 @@
 
 
 class SFMSCalcJob:
-    async def calculate_bui(self, start_time: datetime):
+    async def calculate_daily_fwi(self, start_time: datetime):
         """
         Entry point for processing SFMS DMC/DC/BUI rasters. To run from a specific date manually in openshift,
         see openshift/sfms-calculate/README.md
@@ -28,17 +28,17 @@
 
         start_exec = get_utc_now()
 
-        bui_processor = BUIDateRangeProcessor(start_time, DAYS_TO_CALCULATE, RasterKeyAddresser())
+        daily_processor = DailyFWIProcessor(start_time, DAYS_TO_CALCULATE, RasterKeyAddresser())
 
         async with S3Client() as s3_client:
-            await bui_processor.process_bui(s3_client, multi_wps_dataset_context, multi_wps_dataset_context)
+            await daily_processor.process(s3_client, multi_wps_dataset_context, multi_wps_dataset_context)
 
         # calculate the execution time.
         execution_time = get_utc_now() - start_exec
         hours, remainder = divmod(execution_time.seconds, 3600)
         minutes, seconds = divmod(remainder, 60)
 
-        logger.info(f"BUI processing finished -- time elapsed {hours} hours, {minutes} minutes, {seconds:.2f} seconds")
+        logger.info(f"Daily FWI processing finished -- time elapsed {hours} hours, {minutes} minutes, {seconds:.2f} seconds")
 
 
 def main():
@@ -56,9 +56,9 @@
         job = SFMSCalcJob()
         loop = asyncio.new_event_loop()
         asyncio.set_event_loop(loop)
-        loop.run_until_complete(job.calculate_bui(start_time))
+        loop.run_until_complete(job.calculate_daily_fwi(start_time))
     except Exception as e:
-        logger.error("An exception occurred while processing DMC/DC/BUI raster calculations", exc_info=e)
+        logger.error("An exception occurred while processing SFMS raster calculations", exc_info=e)
         rc_message = ":scream: Encountered an error while processing SFMS raster data."
         send_rocketchat_notification(rc_message, e)
         sys.exit(os.EX_SOFTWARE)

api/app/sfms/date_range_processor.py → api/app/sfms/daily_fwi_processor.py

@@ -8,7 +8,7 @@
 
 from app.geospatial.wps_dataset import WPSDataset
 from app.sfms.raster_addresser import FWIParameter, RasterKeyAddresser
-from app.sfms.fwi_processor import calculate_bui, calculate_dc, calculate_dmc
+from app.sfms.fwi_processor import calculate_bui, calculate_dc, calculate_dmc, calculate_ffmc
 from app.utils.geospatial import GDALResamplingMethod
 from app.utils.s3 import set_s3_gdal_config
 from app.utils.s3_client import S3Client
@@ -20,48 +20,49 @@
 MultiDatasetContext = Callable[[List[str]], Iterator[List["WPSDataset"]]]
 
 
-class BUIDateRangeProcessor:
+class DailyFWIProcessor:
     """
-    Class for calculating/generating forecasted DMC/DC/BUI rasters for a date range
+    Class for calculating/generating forecasted daily FWI rasters for a date range
     """
 
     def __init__(self, start_datetime: datetime, days: int, addresser: RasterKeyAddresser):
         self.start_datetime = start_datetime
         self.days = days
         self.addresser = addresser
 
-    async def process_bui(self, s3_client: S3Client, input_dataset_context: MultiDatasetContext, new_dmc_dc_context: MultiDatasetContext):
+    async def process(self, s3_client: S3Client, input_dataset_context: MultiDatasetContext, new_dmc_dc_context: MultiDatasetContext):
         set_s3_gdal_config()
 
         for day in range(self.days):
             datetime_to_calculate_utc, previous_fwi_datetime, prediction_hour = self._get_calculate_dates(day)
             logger.info(f"Calculating DMC/DC/BUI for {datetime_to_calculate_utc.isoformat()}")
 
             # Get and check existence of weather s3 keys
-            temp_key, rh_key, _, precip_key = self.addresser.get_weather_data_keys(self.start_datetime, datetime_to_calculate_utc, prediction_hour)
+            temp_key, rh_key, wind_speed_key, precip_key = self.addresser.get_weather_data_keys(self.start_datetime, datetime_to_calculate_utc, prediction_hour)
             weather_keys_exist = await s3_client.all_objects_exist(temp_key, rh_key, precip_key)
             if not weather_keys_exist:
                 logging.warning(f"No weather keys found for {model_run_for_hour(self.start_datetime.hour):02} model run")
                 break
 
             # get and check existence of fwi s3 keys
-            dc_key, dmc_key = self._get_previous_fwi_keys(day, previous_fwi_datetime)
+            dc_key, dmc_key, ffmc_key = self._get_previous_fwi_keys(day, previous_fwi_datetime)
             fwi_keys_exist = await s3_client.all_objects_exist(dc_key, dmc_key)
             if not fwi_keys_exist:
                 logging.warning(f"No previous DMC/DC keys found for {previous_fwi_datetime.date().isoformat()}")
                 break
 
-            temp_key, rh_key, precip_key = self.addresser.gdal_prefix_keys(temp_key, rh_key, precip_key)
+            temp_key, rh_key, wind_speed_key, precip_key = self.addresser.gdal_prefix_keys(temp_key, rh_key, wind_speed_key, precip_key)
             dc_key, dmc_key = self.addresser.gdal_prefix_keys(dc_key, dmc_key)
 
             with tempfile.TemporaryDirectory() as temp_dir:
-                with input_dataset_context([temp_key, rh_key, precip_key, dc_key, dmc_key]) as input_datasets:
+                with input_dataset_context([temp_key, rh_key, wind_speed_key, precip_key, dc_key, dmc_key, ffmc_key]) as input_datasets:
                     input_datasets = cast(List[WPSDataset], input_datasets)  # Ensure correct type inference
-                    temp_ds, rh_ds, precip_ds, dc_ds, dmc_ds = input_datasets
+                    temp_ds, rh_ds, wind_speed_ds, precip_ds, dc_ds, dmc_ds, ffmc_ds = input_datasets
 
                     # Warp weather datasets to match fwi
                     warped_temp_ds = temp_ds.warp_to_match(dmc_ds, f"{temp_dir}/{os.path.basename(temp_key)}", GDALResamplingMethod.BILINEAR)
                     warped_rh_ds = rh_ds.warp_to_match(dmc_ds, f"{temp_dir}/{os.path.basename(rh_key)}", GDALResamplingMethod.BILINEAR)
+                    warped_wind_speed_ds = wind_speed_ds.warp_to_match(dmc_ds, f"{temp_dir}/{os.path.basename(wind_speed_key)}", GDALResamplingMethod.BILINEAR)
                     warped_precip_ds = precip_ds.warp_to_match(dmc_ds, f"{temp_dir}/{os.path.basename(precip_key)}", GDALResamplingMethod.BILINEAR)
 
                     # close unneeded datasets to reduce memory usage
@@ -96,6 +97,18 @@ async def process_bui(self, s3_client: S3Client, input_dataset_context: MultiDat
                         dc_nodata_value,
                     )
 
+                    # Create and store FFMC dataset
+                    ffmc_values, ffmc_no_data_value = calculate_ffmc(ffmc_ds, warped_temp_ds, warped_rh_ds, warped_wind_speed_ds, warped_precip_ds)
+                    new_ffmc_key = self.addresser.get_calculated_index_key(datetime_to_calculate_utc, FWIParameter.FFMC)
+                    await s3_client.persist_raster_data(
+                        temp_dir,
+                        new_ffmc_key,
+                        dc_ds.as_gdal_ds().GetGeoTransform(),
+                        dc_ds.as_gdal_ds().GetProjection(),
+                        ffmc_values,
+                        ffmc_no_data_value,
+                    )
+
                     # Open new DMC and DC datasets and calculate BUI
                     new_bui_key = self.addresser.get_calculated_index_key(datetime_to_calculate_utc, FWIParameter.BUI)
                     with new_dmc_dc_context([new_dmc_path, new_dc_path]) as new_dmc_dc_datasets:
@@ -137,7 +150,10 @@ def _get_previous_fwi_keys(self, day_to_calculate: int, previous_fwi_datetime: d
         if day_to_calculate == 0:  # if we're running the first day of the calculation, use previously uploaded actuals
             dc_key = self.addresser.get_uploaded_index_key(previous_fwi_datetime, FWIParameter.DC)
             dmc_key = self.addresser.get_uploaded_index_key(previous_fwi_datetime, FWIParameter.DMC)
+            ffmc_key = self.addresser.get_uploaded_index_key(previous_fwi_datetime, FWIParameter.FFMC)
         else:  # otherwise use the last calculated key
             dc_key = self.addresser.get_calculated_index_key(previous_fwi_datetime, FWIParameter.DC)
             dmc_key = self.addresser.get_calculated_index_key(previous_fwi_datetime, FWIParameter.DMC)
-        return dc_key, dmc_key
+            ffmc_key = self.addresser.get_calculated_index_key(previous_fwi_datetime, FWIParameter.FFMC)
+
+        return dc_key, dmc_key, ffmc_key

api/app/sfms/fwi_processor.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 from app.geospatial.wps_dataset import WPSDataset
-from app.auto_spatial_advisory.sfms import vectorized_dmc, vectorized_dc, vectorized_bui
+from app.auto_spatial_advisory.sfms import vectorized_dmc, vectorized_dc, vectorized_bui, vectorized_ffmc
 
 logger = logging.getLogger(__name__)
 
@@ -55,3 +55,21 @@ def calculate_bui(dmc_ds: WPSDataset, dc_ds: WPSDataset):
         bui_values[nodata_mask] = nodata_value
 
     return bui_values, nodata_value
+
+
+def calculate_ffmc(previous_ffmc_ds: WPSDataset, temp_ds: WPSDataset, rh_ds: WPSDataset, precip_ds: WPSDataset, wind_speed_ds: WPSDataset):
+    previous_ffmc_array, _ = previous_ffmc_ds.replace_nodata_with(0)
+    temp_array, _ = temp_ds.replace_nodata_with(0)
+    rh_array, _ = rh_ds.replace_nodata_with(0)
+    precip_array, _ = precip_ds.replace_nodata_with(0)
+    wind_speed_array, _ = wind_speed_ds.replace_nodata_with(0)
+
+    start = perf_counter()
+    ffmc_values = vectorized_ffmc(previous_ffmc_array, temp_array, rh_array, precip_array, wind_speed_array)
+    logger.info("%f seconds to calculate vectorized ffmc", perf_counter() - start)
+
+    nodata_mask, nodata_value = previous_ffmc_ds.get_nodata_mask()
+    if nodata_mask is not None:
+        ffmc_values[nodata_mask] = nodata_value
+
+    return ffmc_values, nodata_value

api/app/sfms/raster_addresser.py

@@ -1,6 +1,6 @@
 import os
 import enum
-from datetime import datetime, timezone, timedelta
+from datetime import datetime, timezone
 from zoneinfo import ZoneInfo
 from app import config
 from app.weather_models import ModelEnum

api/app/tests/jobs/test_sfms_calculations.py

@@ -12,7 +12,7 @@ def test_sfms_calc_job_fail_default(monkeypatch, mocker: MockerFixture):
     async def mock_job_error():
         raise OSError("Error")
 
-    monkeypatch.setattr(SFMSCalcJob, "calculate_bui", mock_job_error)
+    monkeypatch.setattr(SFMSCalcJob, "calculate_daily_fwi", mock_job_error)
 
     monkeypatch.setattr("sys.argv", ["sfms_calculations.py"])
 
@@ -27,15 +27,14 @@ async def mock_job_error():
 
 
 def test_sfms_calc_job_cli_arg(monkeypatch, mocker: MockerFixture):
-    calc_spy = mocker.patch.object(SFMSCalcJob, "calculate_bui", return_value=None)
+    daily_fwi_calc_spy = mocker.patch.object(SFMSCalcJob, "calculate_daily_fwi", return_value=None)
 
     test_datetime = "2024-10-10 5"
     monkeypatch.setattr("sys.argv", ["sfms_calculations.py", test_datetime])
 
     sfms_calculations.main()
 
-    called_args, _ = calc_spy.call_args
-    assert called_args[0] == datetime.strptime(test_datetime, "%Y-%m-%d %H").replace(tzinfo=timezone.utc)
+    daily_fwi_calc_spy.assert_called_once_with(datetime.strptime(test_datetime, "%Y-%m-%d %H").replace(tzinfo=timezone.utc))
 
 
 @pytest.mark.anyio

api/app/tests/sfms/test_bui_date_range_processor.py → api/app/tests/sfms/test_daily_fwi_processor.py

@@ -5,8 +5,8 @@
 from datetime import datetime, timezone, timedelta
 from pytest_mock import MockerFixture
 from app.geospatial.wps_dataset import WPSDataset
-from app.sfms import date_range_processor
-from app.sfms.date_range_processor import BUIDateRangeProcessor
+from app.sfms import daily_fwi_processor
+from app.sfms.daily_fwi_processor import DailyFWIProcessor
 from app.sfms.raster_addresser import FWIParameter, RasterKeyAddresser
 from app.tests.dataset_common import create_mock_gdal_dataset, create_mock_wps_dataset
 from app.utils.geospatial import GDALResamplingMethod
@@ -22,7 +22,7 @@ def create_mock_wps_datasets(num: int) -> List[WPSDataset]:
 
 
 def create_mock_input_dataset_context():
-    input_datasets = create_mock_wps_datasets(5)
+    input_datasets = create_mock_wps_datasets(7)
 
     @contextmanager
     def mock_input_dataset_context(_: List[str]):
@@ -62,14 +62,15 @@ async def test_bui_date_range_processor(mocker: MockerFixture):
     get_weather_data_key_spy = mocker.spy(mock_key_addresser, "get_weather_data_keys")
     gdal_prefix_keys_spy = mocker.spy(mock_key_addresser, "gdal_prefix_keys")
     get_calculated_index_key_spy = mocker.spy(mock_key_addresser, "get_calculated_index_key")
-    bui_date_range_processor = BUIDateRangeProcessor(TEST_DATETIME, 2, mock_key_addresser)
+    bui_date_range_processor = DailyFWIProcessor(TEST_DATETIME, 2, mock_key_addresser)
     # mock/spy dataset storage
 
     # mock weather index, param datasets used for calculations
     input_datasets, mock_input_dataset_context = create_mock_input_dataset_context()
-    mock_temp_ds, mock_rh_ds, mock_precip_ds, mock_dc_ds, mock_dmc_ds = input_datasets
+    mock_temp_ds, mock_rh_ds, mock_precip_ds, mock_wind_speed_ds, mock_dc_ds, mock_dmc_ds, mock_ffmc_ds = input_datasets
     temp_ds_spy = mocker.spy(mock_temp_ds, "warp_to_match")
     rh_ds_spy = mocker.spy(mock_rh_ds, "warp_to_match")
+    wind_speed_ds_spy = mocker.spy(mock_wind_speed_ds, "warp_to_match")
     precip_ds_spy = mocker.spy(mock_precip_ds, "warp_to_match")
 
     # mock new dmc and dc datasets
@@ -80,17 +81,18 @@ async def test_bui_date_range_processor(mocker: MockerFixture):
     mocker.patch("osgeo.gdal.Open", return_value=create_mock_gdal_dataset())
 
     # calculation spies
-    calculate_dmc_spy = mocker.spy(date_range_processor, "calculate_dmc")
-    calculate_dc_spy = mocker.spy(date_range_processor, "calculate_dc")
-    calculate_bui_spy = mocker.spy(date_range_processor, "calculate_bui")
+    calculate_dmc_spy = mocker.spy(daily_fwi_processor, "calculate_dmc")
+    calculate_dc_spy = mocker.spy(daily_fwi_processor, "calculate_dc")
+    calculate_bui_spy = mocker.spy(daily_fwi_processor, "calculate_bui")
+    calculate_ffmc_spy = mocker.spy(daily_fwi_processor, "calculate_ffmc")
 
     async with S3Client() as mock_s3_client:
         # mock s3 client
         mock_all_objects_exist = AsyncMock(return_value=True)
         mocker.patch.object(mock_s3_client, "all_objects_exist", new=mock_all_objects_exist)
         persist_raster_spy = mocker.patch.object(mock_s3_client, "persist_raster_data", return_value="test_key.tif")
 
-        await bui_date_range_processor.process_bui(mock_s3_client, mock_input_dataset_context, mock_new_dmc_dc_datasets_context)
+        await bui_date_range_processor.process(mock_s3_client, mock_input_dataset_context, mock_new_dmc_dc_datasets_context)
 
         # Verify weather model keys and actual keys are checked for both days
         assert mock_all_objects_exist.call_count == 4
@@ -107,6 +109,7 @@ async def test_bui_date_range_processor(mocker: MockerFixture):
             mocker.call(
                 "weather_models/rdps/2024-10-10/00/temp/CMC_reg_TMP_TGL_2_ps10km_2024101000_P020.grib2",
                 "weather_models/rdps/2024-10-10/00/rh/CMC_reg_RH_TGL_2_ps10km_2024101000_P020.grib2",
+                "weather_models/rdps/2024-10-10/00/wind_speed/CMC_reg_WIND_TGL_10_ps10km_2024101000_P020.grib2",
                 "weather_models/rdps/2024-10-10/12/precip/COMPUTED_reg_APCP_SFC_0_ps10km_20241010_20z.tif",
             ),
             # first day uploads
@@ -115,6 +118,7 @@ async def test_bui_date_range_processor(mocker: MockerFixture):
             mocker.call(
                 "weather_models/rdps/2024-10-10/00/temp/CMC_reg_TMP_TGL_2_ps10km_2024101000_P044.grib2",
                 "weather_models/rdps/2024-10-10/00/rh/CMC_reg_RH_TGL_2_ps10km_2024101000_P044.grib2",
+                "weather_models/rdps/2024-10-10/00/wind_speed/CMC_reg_WIND_TGL_10_ps10km_2024101000_P044.grib2",
                 "weather_models/rdps/2024-10-11/12/precip/COMPUTED_reg_APCP_SFC_0_ps10km_20241011_20z.tif",
             ),
             # second day uploads
@@ -126,12 +130,15 @@ async def test_bui_date_range_processor(mocker: MockerFixture):
             # first day
             mocker.call(EXPECTED_FIRST_DAY, FWIParameter.DMC),
             mocker.call(EXPECTED_FIRST_DAY, FWIParameter.DC),
+            mocker.call(EXPECTED_FIRST_DAY, FWIParameter.FFMC),
             mocker.call(EXPECTED_FIRST_DAY, FWIParameter.BUI),
             # second day, previous days' dc and dmc are looked up first
             mocker.call(EXPECTED_FIRST_DAY, FWIParameter.DC),
             mocker.call(EXPECTED_FIRST_DAY, FWIParameter.DMC),
+            mocker.call(EXPECTED_FIRST_DAY, FWIParameter.FFMC),
             mocker.call(EXPECTED_SECOND_DAY, FWIParameter.DMC),
             mocker.call(EXPECTED_SECOND_DAY, FWIParameter.DC),
+            mocker.call(EXPECTED_SECOND_DAY, FWIParameter.FFMC),
             mocker.call(EXPECTED_SECOND_DAY, FWIParameter.BUI),
         ]
 
@@ -146,6 +153,11 @@ async def test_bui_date_range_processor(mocker: MockerFixture):
             mocker.call(mock_dmc_ds, mocker.ANY, GDALResamplingMethod.BILINEAR),
         ]
 
+        assert wind_speed_ds_spy.call_args_list == [
+            mocker.call(mock_dmc_ds, mocker.ANY, GDALResamplingMethod.BILINEAR),
+            mocker.call(mock_dmc_ds, mocker.ANY, GDALResamplingMethod.BILINEAR),
+        ]
+
         assert precip_ds_spy.call_args_list == [
             mocker.call(mock_dmc_ds, mocker.ANY, GDALResamplingMethod.BILINEAR),
             mocker.call(mock_dmc_ds, mocker.ANY, GDALResamplingMethod.BILINEAR),
@@ -163,13 +175,19 @@ async def test_bui_date_range_processor(mocker: MockerFixture):
             wps_datasets = dc_calls[0][1:4]  # Extract dataset arguments
             assert all(isinstance(ds, WPSDataset) for ds in wps_datasets)
 
+        for ffmc_calls in calculate_ffmc_spy.call_args_list:
+            ffmc_ds = ffmc_calls[0][0]
+            assert ffmc_ds == mock_ffmc_ds
+            wps_datasets = ffmc_calls[0][1:4]  # Extract dataset arguments
+            assert all(isinstance(ds, WPSDataset) for ds in wps_datasets)
+
         assert calculate_bui_spy.call_args_list == [
             mocker.call(mock_new_dmc_ds, mock_new_dc_ds),
             mocker.call(mock_new_dmc_ds, mock_new_dc_ds),
         ]
 
-        # 3 each day, new dmc, dc and bui rasters
-        assert persist_raster_spy.call_count == 6
+        # 4 each day, new dmc, dc and bui rasters
+        assert persist_raster_spy.call_count == 8
 
 
 @pytest.mark.parametrize(
@@ -191,13 +209,13 @@ async def test_no_weather_keys_exist(side_effect_1: bool, side_effect_2: bool, m
     _, mock_new_dmc_dc_datasets_context = create_mock_new_dmc_dc_context()
 
     # calculation spies
-    calculate_dmc_spy = mocker.spy(date_range_processor, "calculate_dmc")
-    calculate_dc_spy = mocker.spy(date_range_processor, "calculate_dc")
-    calculate_bui_spy = mocker.spy(date_range_processor, "calculate_bui")
+    calculate_dmc_spy = mocker.spy(daily_fwi_processor, "calculate_dmc")
+    calculate_dc_spy = mocker.spy(daily_fwi_processor, "calculate_dc")
+    calculate_bui_spy = mocker.spy(daily_fwi_processor, "calculate_bui")
 
-    bui_date_range_processor = BUIDateRangeProcessor(TEST_DATETIME, 1, RasterKeyAddresser())
+    bui_date_range_processor = DailyFWIProcessor(TEST_DATETIME, 1, RasterKeyAddresser())
 
-    await bui_date_range_processor.process_bui(mock_s3_client, mock_input_dataset_context, mock_new_dmc_dc_datasets_context)
+    await bui_date_range_processor.process(mock_s3_client, mock_input_dataset_context, mock_new_dmc_dc_datasets_context)
 
     calculate_dmc_spy.assert_not_called()
     calculate_dc_spy.assert_not_called()