diff --git a/.github/workflows/integration.yml b/.github/workflows/integration.yml index d7d438b44..d723feec1 100644 --- a/.github/workflows/integration.yml +++ b/.github/workflows/integration.yml @@ -115,13 +115,27 @@ jobs: with: fetch-depth: 0 - uses: ./.github/actions/python-setup - with: - install-jdk: "true" - name: Unit Test (wps_weather) working-directory: ./backend run: | uv run pytest packages/wps-weather/src/wps_weather/tests -x -o log_cli=true --disable-warnings -vvv + test-wps-sfms: + name: Python - WPS SFMS Test + runs-on: ubuntu-24.04 + container: + image: ghcr.io/bcgov/wps/wps-api-base:01-29-2026 + options: --user 0 + steps: + - uses: actions/checkout@v6 + with: + fetch-depth: 0 + - uses: ./.github/actions/python-setup + - name: Unit Test (wps_sfms) + working-directory: ./backend + run: | + uv run pytest packages/wps-sfms/src/wps_sfms/tests -x -o log_cli=true --disable-warnings -vvv + test-wps-wf1: name: Python - WF1 Test with coverage runs-on: ubuntu-24.04 diff --git a/backend/packages/wps-api/src/app/jobs/sfms_daily_actuals.py b/backend/packages/wps-api/src/app/jobs/sfms_daily_actuals.py index 18c67ea0e..2d1316cbc 100644 --- a/backend/packages/wps-api/src/app/jobs/sfms_daily_actuals.py +++ b/backend/packages/wps-api/src/app/jobs/sfms_daily_actuals.py @@ -1,6 +1,5 @@ """ -Job for running SFMS daily actuals: temperature interpolation followed by -precipitation interpolation for the current date. +Job for running SFMS daily actuals weather interpolation and FWI processing. Usage: python -m app.jobs.sfms_daily_actuals "YYYY-MM-DD" @@ -8,25 +7,30 @@ """ import asyncio +from dataclasses import dataclass import logging import os import sys from datetime import datetime, timezone +from typing import Awaitable, Callable from aiohttp import ClientSession -from wps_sfms.interpolation.source import ( - StationDCSource, - StationDewPointSource, - StationDMCSource, - StationFFMCSource, - StationPrecipitationSource, - StationTemperatureSource, +from wps_sfms.interpolation.field import ( + build_dc_field, + build_dewpoint_field, + build_dmc_field, + build_ffmc_field, + build_precipitation_field, + build_temperature_field, + build_wind_speed_field, + build_wind_vector_field, ) from wps_sfms.processors.fwi import DCCalculator, DMCCalculator, FFMCCalculator, FWIProcessor -from wps_sfms.processors.idw import Interpolator +from wps_sfms.processors.idw import Interpolator, RasterProcessor from wps_sfms.processors.relative_humidity import RHInterpolator from wps_sfms.processors.temperature import TemperatureInterpolator +from wps_sfms.processors.wind import WindDirectionInterpolator, WindSpeedInterpolator from wps_shared.db.crud.fuel_layer import get_fuel_type_raster_by_year from wps_shared.db.crud.sfms_run import save_sfms_run, track_sfms_run from wps_shared.db.database import get_async_read_session_scope, get_async_write_session_scope @@ -43,6 +47,20 @@ logger = logging.getLogger(__name__) +@dataclass(frozen=True) +class RasterInterpolationJob: + job_name: SFMSRunLogJobName + output_key: str + log_label: str + processor: RasterProcessor + + +@dataclass(frozen=True) +class FWICalculationJob: + job_name: SFMSRunLogJobName + calculator: FFMCCalculator | DMCCalculator | DCCalculator + + def is_fwi_interpolation_day(dt: datetime) -> bool: """Return True if FWI indices should be re-interpolated from station observations. @@ -53,6 +71,37 @@ def is_fwi_interpolation_day(dt: datetime) -> bool: return dt.weekday() == 0 and dt.month in (4, 5) +async def _run_tracked_job( + job_name: SFMSRunLogJobName, + sfms_run_id: int, + session, + action: Callable[[], Awaitable[object]], +): + @track_sfms_run(job_name, sfms_run_id, session) + async def _wrapped(): + return await action() + + return await _wrapped() + + +async def _process_raster_job( + *, + job_name: SFMSRunLogJobName, + sfms_run_id: int, + session, + processor: RasterProcessor, + s3_client: S3Client, + fuel_raster_path: str, + output_key: str, + log_label: str, +) -> None: + async def _run() -> None: + s3_key = await processor.process(s3_client, fuel_raster_path, output_key) + logger.info("%s: %s", log_label, s3_key) + + await _run_tracked_job(job_name, sfms_run_id, session, _run) + + async def run_weather_interpolation( datetime_to_process: datetime, raster_addresser: SFMSNGRasterAddresser, @@ -62,50 +111,69 @@ async def run_weather_interpolation( sfms_run_id: int, session, ) -> None: - """Interpolate temperature, RH, and precipitation from station observations.""" + """Interpolate weather rasters from station observations.""" mask_path = raster_addresser.get_mask_key() dem_path = raster_addresser.get_dem_key() - temp_key = raster_addresser.get_actual_weather_key( + temp_output_key = raster_addresser.get_actual_weather_key( datetime_to_process, SFMSInterpolatedWeatherParameter.TEMP ) - temp_processor = TemperatureInterpolator(mask_path, dem_path) - rh_processor = RHInterpolator(mask_path, dem_path, raster_addresser.gdal_path(temp_key)) - precip_processor = Interpolator(mask_path) - - @track_sfms_run(SFMSRunLogJobName.TEMPERATURE_INTERPOLATION, sfms_run_id, session) - async def run_temperature_interpolation() -> None: - temp_s3_key = await temp_processor.process( - s3_client, fuel_raster_path, StationTemperatureSource(sfms_actuals), temp_key - ) - logger.info("Temperature interpolation raster: %s", temp_s3_key) - - @track_sfms_run(SFMSRunLogJobName.RH_INTERPOLATION, sfms_run_id, session) - async def run_rh_interpolation() -> None: - rh_s3_key = await rh_processor.process( - s3_client, - fuel_raster_path, - StationDewPointSource(sfms_actuals), - raster_addresser.get_actual_weather_key( + temp_raster_path = raster_addresser.gdal_path(temp_output_key) + jobs = [ + RasterInterpolationJob( + job_name=SFMSRunLogJobName.TEMPERATURE_INTERPOLATION, + output_key=temp_output_key, + log_label="Temperature interpolation raster", + processor=TemperatureInterpolator( + mask_path, dem_path, build_temperature_field(sfms_actuals) + ), + ), + RasterInterpolationJob( + job_name=SFMSRunLogJobName.RH_INTERPOLATION, + output_key=raster_addresser.get_actual_weather_key( datetime_to_process, SFMSInterpolatedWeatherParameter.RH ), - ) - logger.info("RH interpolation raster: %s", rh_s3_key) - - @track_sfms_run(SFMSRunLogJobName.PRECIPITATION_INTERPOLATION, sfms_run_id, session) - async def run_precipitation_interpolation() -> None: - precip_s3_key = await precip_processor.process( - s3_client, - fuel_raster_path, - StationPrecipitationSource(sfms_actuals), - raster_addresser.get_actual_weather_key( + log_label="RH interpolation raster", + processor=RHInterpolator( + mask_path, dem_path, temp_raster_path, build_dewpoint_field(sfms_actuals) + ), + ), + RasterInterpolationJob( + job_name=SFMSRunLogJobName.WIND_SPEED_INTERPOLATION, + output_key=raster_addresser.get_actual_weather_key( + datetime_to_process, SFMSInterpolatedWeatherParameter.WIND_SPEED + ), + log_label="Wind speed interpolation raster", + processor=WindSpeedInterpolator(mask_path, build_wind_speed_field(sfms_actuals)), + ), + RasterInterpolationJob( + job_name=SFMSRunLogJobName.WIND_DIRECTION_INTERPOLATION, + output_key=raster_addresser.get_actual_weather_key( + datetime_to_process, SFMSInterpolatedWeatherParameter.WIND_DIRECTION + ), + log_label="Wind direction interpolation raster", + processor=WindDirectionInterpolator(mask_path, build_wind_vector_field(sfms_actuals)), + ), + RasterInterpolationJob( + job_name=SFMSRunLogJobName.PRECIPITATION_INTERPOLATION, + output_key=raster_addresser.get_actual_weather_key( datetime_to_process, SFMSInterpolatedWeatherParameter.PRECIP ), - ) - logger.info("Precip interpolation raster: %s", precip_s3_key) + log_label="Precip interpolation raster", + processor=Interpolator(mask_path, build_precipitation_field(sfms_actuals)), + ), + ] - await run_temperature_interpolation() - await run_rh_interpolation() - await run_precipitation_interpolation() + for job in jobs: + await _process_raster_job( + job_name=job.job_name, + sfms_run_id=sfms_run_id, + session=session, + processor=job.processor, + s3_client=s3_client, + fuel_raster_path=fuel_raster_path, + output_key=job.output_key, + log_label=job.log_label, + ) async def run_fwi_interpolation( @@ -123,23 +191,41 @@ async def run_fwi_interpolation( ) mask_path = raster_addresser.get_mask_key() - fwi_sources = [ - (SFMSRunLogJobName.FFMC_INTERPOLATION, StationFFMCSource(sfms_actuals), FWIParameter.FFMC), - (SFMSRunLogJobName.DMC_INTERPOLATION, StationDMCSource(sfms_actuals), FWIParameter.DMC), - (SFMSRunLogJobName.DC_INTERPOLATION, StationDCSource(sfms_actuals), FWIParameter.DC), + ffmc_output_key = raster_addresser.get_actual_index_key(datetime_to_process, FWIParameter.FFMC) + dmc_output_key = raster_addresser.get_actual_index_key(datetime_to_process, FWIParameter.DMC) + dc_output_key = raster_addresser.get_actual_index_key(datetime_to_process, FWIParameter.DC) + jobs = [ + RasterInterpolationJob( + job_name=SFMSRunLogJobName.FFMC_INTERPOLATION, + output_key=ffmc_output_key, + log_label=f"{SFMSRunLogJobName.FFMC_INTERPOLATION.value} raster", + processor=Interpolator(mask_path, build_ffmc_field(sfms_actuals)), + ), + RasterInterpolationJob( + job_name=SFMSRunLogJobName.DMC_INTERPOLATION, + output_key=dmc_output_key, + log_label=f"{SFMSRunLogJobName.DMC_INTERPOLATION.value} raster", + processor=Interpolator(mask_path, build_dmc_field(sfms_actuals)), + ), + RasterInterpolationJob( + job_name=SFMSRunLogJobName.DC_INTERPOLATION, + output_key=dc_output_key, + log_label=f"{SFMSRunLogJobName.DC_INTERPOLATION.value} raster", + processor=Interpolator(mask_path, build_dc_field(sfms_actuals)), + ), ] - processor = Interpolator(mask_path) - - for job_name, source, fwi_param in fwi_sources: - - @track_sfms_run(job_name, sfms_run_id, session) - async def _run(_source=source, _job_name=job_name, _fwi_param=fwi_param) -> None: - output_key = raster_addresser.get_actual_index_key(datetime_to_process, _fwi_param) - s3_key = await processor.process(s3_client, fuel_raster_path, _source, output_key) - logger.info("%s interpolation raster: %s", _job_name.value, s3_key) - - await _run() + for job in jobs: + await _process_raster_job( + job_name=job.job_name, + sfms_run_id=sfms_run_id, + session=session, + processor=job.processor, + s3_client=s3_client, + fuel_raster_path=fuel_raster_path, + output_key=job.output_key, + log_label=job.log_label, + ) async def run_fwi_calculations( @@ -158,16 +244,24 @@ async def run_fwi_calculations( fwi_processor = FWIProcessor(datetime_to_process) month = datetime_to_process.month - fwi_calculations = [ - (SFMSRunLogJobName.FFMC_CALCULATION, FFMCCalculator()), - (SFMSRunLogJobName.DMC_CALCULATION, DMCCalculator(month)), - (SFMSRunLogJobName.DC_CALCULATION, DCCalculator(month)), + jobs = [ + FWICalculationJob( + job_name=SFMSRunLogJobName.FFMC_CALCULATION, + calculator=FFMCCalculator(), + ), + FWICalculationJob( + job_name=SFMSRunLogJobName.DMC_CALCULATION, + calculator=DMCCalculator(month), + ), + FWICalculationJob( + job_name=SFMSRunLogJobName.DC_CALCULATION, + calculator=DCCalculator(month), + ), ] - for job_name, calculator in fwi_calculations: + for job in jobs: - @track_sfms_run(job_name, sfms_run_id, session) - async def _run(_calculator=calculator) -> None: + async def _run(_calculator=job.calculator) -> None: _fwi_inputs = raster_addresser.get_actual_fwi_inputs( datetime_to_process, _calculator.fwi_param ) @@ -175,11 +269,11 @@ async def _run(_calculator=calculator) -> None: s3_client, multi_wps_dataset_context, _calculator, _fwi_inputs ) - await _run() + await _run_tracked_job(job.job_name, sfms_run_id, session, _run) async def run_sfms_daily_actuals(target_date: datetime) -> None: - """Run temperature then precipitation interpolation for the given date.""" + """Run SFMS daily weather interpolation and FWI updates for the given date.""" logger.info("Starting SFMS daily actuals for %s", target_date.date()) raster_addresser = SFMSNGRasterAddresser() diff --git a/backend/packages/wps-api/src/app/tests/conftest.py b/backend/packages/wps-api/src/app/tests/conftest.py index 3923ee259..ec9ec713d 100644 --- a/backend/packages/wps-api/src/app/tests/conftest.py +++ b/backend/packages/wps-api/src/app/tests/conftest.py @@ -33,6 +33,7 @@ def create_mock_sfms_actuals(): relative_humidity=50.0, precipitation=2.5, wind_speed=10.0, + wind_direction=180.0, ffmc=85.0, dmc=30.0, dc=200.0, @@ -46,6 +47,7 @@ def create_mock_sfms_actuals(): relative_humidity=60.0, precipitation=5.0, wind_speed=8.0, + wind_direction=200.0, ffmc=80.0, dmc=25.0, dc=180.0, diff --git a/backend/packages/wps-api/src/app/tests/jobs/test_sfms_daily_actuals.py b/backend/packages/wps-api/src/app/tests/jobs/test_sfms_daily_actuals.py index 8d697d6ab..98ea83b77 100644 --- a/backend/packages/wps-api/src/app/tests/jobs/test_sfms_daily_actuals.py +++ b/backend/packages/wps-api/src/app/tests/jobs/test_sfms_daily_actuals.py @@ -4,6 +4,7 @@ import sys from contextlib import asynccontextmanager from datetime import datetime, timezone +from types import SimpleNamespace from typing import NamedTuple from unittest.mock import AsyncMock, MagicMock @@ -17,6 +18,7 @@ from wps_sfms.processors.idw import Interpolator from wps_sfms.processors.relative_humidity import RHInterpolator from wps_sfms.processors.temperature import TemperatureInterpolator +from wps_sfms.processors.wind import WindDirectionInterpolator, WindSpeedInterpolator from wps_shared.db.models.sfms_run import SFMSRunLogStatus from wps_shared.sfms.raster_addresser import FWIParameter @@ -27,9 +29,9 @@ "dt,expected", [ # True cases: Mondays in April or May (any week) - (datetime(2024, 4, 1, tzinfo=timezone.utc), True), # Monday, April, day 1 - (datetime(2024, 4, 8, tzinfo=timezone.utc), True), # Monday, April, day 8 (second Monday) - (datetime(2024, 5, 6, tzinfo=timezone.utc), True), # Monday, May + (datetime(2024, 4, 1, tzinfo=timezone.utc), True), # Monday, April, day 1 + (datetime(2024, 4, 8, tzinfo=timezone.utc), True), # Monday, April, day 8 (second Monday) + (datetime(2024, 5, 6, tzinfo=timezone.utc), True), # Monday, May (datetime(2024, 5, 27, tzinfo=timezone.utc), True), # Monday, May, last Monday # False cases: wrong weekday (datetime(2024, 4, 2, tzinfo=timezone.utc), False), # Tuesday, April @@ -51,6 +53,8 @@ class MockDailyActualsDeps(NamedTuple): s3_client: AsyncMock temp_processor: MagicMock rh_processor: MagicMock + wind_speed_processor: MagicMock + wind_direction_processor: MagicMock interpolation_processor: MagicMock fwi_processor: MagicMock wfwx_api: AsyncMock @@ -73,6 +77,12 @@ def mock_dependencies(mocker: MockerFixture, mock_s3_client, mock_wfwx_api) -> M mock_wfwx_api.get_sfms_daily_actuals_all_stations = AsyncMock( return_value=create_mock_sfms_actuals() ) + mock_wfwx_api.get_station_data = AsyncMock( + return_value=[ + SimpleNamespace(code=100, name="Vancouver Harbour"), + SimpleNamespace(code=101, name="Burnaby Mountain"), + ] + ) mocker.patch(f"{MODULE_PATH}.WfwxApi", return_value=mock_wfwx_api) # Mock get_fuel_type_raster_by_year @@ -100,7 +110,6 @@ async def _read_scope(): mock_addresser = MagicMock() mock_addresser.s3_prefix = "/vsis3/test-bucket" mocker.patch(f"{MODULE_PATH}.SFMSNGRasterAddresser", return_value=mock_addresser) - # Mock processors mock_temp_processor = MagicMock(spec=TemperatureInterpolator) mock_temp_processor.process = AsyncMock(return_value="sfms/interpolated/2024/07/04/temp.tif") @@ -116,6 +125,24 @@ async def _read_scope(): return_value=mock_rh_processor, ) + mock_wind_speed_processor = MagicMock(spec=WindSpeedInterpolator) + mock_wind_speed_processor.process = AsyncMock( + return_value="sfms/interpolated/2024/07/04/wind_speed.tif" + ) + mocker.patch( + f"{MODULE_PATH}.WindSpeedInterpolator", + return_value=mock_wind_speed_processor, + ) + + mock_wind_direction_processor = MagicMock(spec=WindDirectionInterpolator) + mock_wind_direction_processor.process = AsyncMock( + return_value="sfms/interpolated/2024/07/04/wind_direction.tif" + ) + mocker.patch( + f"{MODULE_PATH}.WindDirectionInterpolator", + return_value=mock_wind_direction_processor, + ) + mock_interpolation_processor = MagicMock(spec=Interpolator) mock_interpolation_processor.process = AsyncMock( return_value="sfms/interpolated/2024/07/04/precip.tif" @@ -144,6 +171,8 @@ async def _scope(): s3_client=mock_s3_client, temp_processor=mock_temp_processor, rh_processor=mock_rh_processor, + wind_speed_processor=mock_wind_speed_processor, + wind_direction_processor=mock_wind_direction_processor, interpolation_processor=mock_interpolation_processor, fwi_processor=mock_fwi_processor, wfwx_api=mock_wfwx_api, @@ -156,20 +185,20 @@ class TestRunSfmsDailyActuals: @pytest.mark.anyio async def test_runs_all_processors(self, mock_dependencies: MockDailyActualsDeps): - """Test that temperature, RH, and precipitation processors are all called.""" + """Test that all weather interpolation processors are called.""" target_date = datetime(2024, 7, 4, tzinfo=timezone.utc) await run_sfms_daily_actuals(target_date) mock_dependencies.temp_processor.process.assert_called_once() mock_dependencies.rh_processor.process.assert_called_once() + mock_dependencies.wind_speed_processor.process.assert_called_once() + mock_dependencies.wind_direction_processor.process.assert_called_once() mock_dependencies.interpolation_processor.process.assert_called_once() @pytest.mark.anyio - async def test_runs_processors_in_order( - self, mock_dependencies: MockDailyActualsDeps - ): - """Test that processors run in order: temperature, RH, precipitation.""" + async def test_runs_processors_in_order(self, mock_dependencies: MockDailyActualsDeps): + """Test that weather processors run in the expected sequence.""" call_order = [] mock_dependencies.temp_processor.process = AsyncMock( side_effect=lambda *a, **kw: call_order.append("temp") or "temp.tif" @@ -177,6 +206,12 @@ async def test_runs_processors_in_order( mock_dependencies.rh_processor.process = AsyncMock( side_effect=lambda *a, **kw: call_order.append("rh") or "rh.tif" ) + mock_dependencies.wind_speed_processor.process = AsyncMock( + side_effect=lambda *a, **kw: call_order.append("wind_speed") or "wind_speed.tif" + ) + mock_dependencies.wind_direction_processor.process = AsyncMock( + side_effect=lambda *a, **kw: call_order.append("wind_direction") or "wind_direction.tif" + ) mock_dependencies.interpolation_processor.process = AsyncMock( side_effect=lambda *a, **kw: call_order.append("precip") or "precip.tif" ) @@ -184,11 +219,11 @@ async def test_runs_processors_in_order( target_date = datetime(2024, 7, 4, tzinfo=timezone.utc) await run_sfms_daily_actuals(target_date) - assert call_order == ["temp", "rh", "precip"] + assert call_order == ["temp", "rh", "wind_speed", "wind_direction", "precip"] @pytest.mark.anyio async def test_passes_s3_client_to_processors(self, mock_dependencies: MockDailyActualsDeps): - """Test that the S3 client is passed to both processors.""" + """Test that the S3 client is passed to weather processors.""" target_date = datetime(2024, 7, 4, tzinfo=timezone.utc) await run_sfms_daily_actuals(target_date) @@ -196,6 +231,9 @@ async def test_passes_s3_client_to_processors(self, mock_dependencies: MockDaily temp_call_args = mock_dependencies.temp_processor.process.call_args assert temp_call_args[0][0] is mock_dependencies.s3_client + wind_speed_call_args = mock_dependencies.wind_speed_processor.process.call_args + assert wind_speed_call_args[0][0] is mock_dependencies.s3_client + precip_call_args = mock_dependencies.interpolation_processor.process.call_args assert precip_call_args[0][0] is mock_dependencies.s3_client @@ -227,19 +265,19 @@ async def test_writes_run_log_entries(self, mock_dependencies: MockDailyActualsD await run_sfms_daily_actuals(target_date) - # Six tracked runs (temp, rh, precip, ffmc, dmc, dc) means six execute calls - assert mock_dependencies.db_session.execute.call_count == 6 + # Eight tracked runs (temp, rh, ws, wd, precip, ffmc, dmc, dc) means eight execute calls + assert mock_dependencies.db_session.execute.call_count == 8 @pytest.mark.anyio async def test_logs_success_status(self, mock_dependencies: MockDailyActualsDeps): """Test that successful jobs are updated to success status.""" - records = [MagicMock() for _ in range(6)] + records = [MagicMock() for _ in range(8)] mock_dependencies.db_session.get = AsyncMock(side_effect=records) target_date = datetime(2024, 7, 4, tzinfo=timezone.utc) await run_sfms_daily_actuals(target_date) - assert mock_dependencies.db_session.get.call_count == 6 + assert mock_dependencies.db_session.get.call_count == 8 for record in records: assert record.status == SFMSRunLogStatus.SUCCESS assert record.completed_at is not None @@ -259,6 +297,8 @@ async def test_temperature_failure_logs_failed_and_raises( await run_sfms_daily_actuals(target_date) mock_dependencies.rh_processor.process.assert_not_called() + mock_dependencies.wind_speed_processor.process.assert_not_called() + mock_dependencies.wind_direction_processor.process.assert_not_called() mock_dependencies.interpolation_processor.process.assert_not_called() @pytest.mark.anyio @@ -276,6 +316,9 @@ async def test_precipitation_failure_logs_failed_and_raises( await run_sfms_daily_actuals(target_date) mock_dependencies.temp_processor.process.assert_called_once() + mock_dependencies.rh_processor.process.assert_called_once() + mock_dependencies.wind_speed_processor.process.assert_called_once() + mock_dependencies.wind_direction_processor.process.assert_called_once() class TestMondayFWIInterpolation: @@ -355,21 +398,23 @@ async def test_second_monday_april_runs_fwi_interpolation( async def test_monday_april_writes_six_run_log_entries( self, mock_dependencies: MockDailyActualsDeps ): - """Test that a Monday in April produces 6 run log entries (temp + rh + precip + 3 FWI).""" + """Test that a Monday in April produces 8 run log entries (temp + rh + precip + ws + wd + 3 FWI).""" # 2024-04-01 is the first Monday of April 2024 target_date = datetime(2024, 4, 1, tzinfo=timezone.utc) await run_sfms_daily_actuals(target_date) - # 6 tracked runs: temp, rh, precip, ffmc, dmc, dc - assert mock_dependencies.db_session.execute.call_count == 6 + # 8 tracked runs: temp, rh, ws, wd, precip, ffmc, dmc, dc + assert mock_dependencies.db_session.execute.call_count == 8 class TestFWICalculationVsInterpolation: """Regression tests: FWI interpolation and FWI calculation are mutually exclusive.""" @pytest.mark.anyio - async def test_monday_interpolation_skips_fwi_calculation(self, mock_dependencies: MockDailyActualsDeps): + async def test_monday_interpolation_skips_fwi_calculation( + self, mock_dependencies: MockDailyActualsDeps + ): """On a re-interpolation Monday, FWI calculation must not run.""" # 2024-05-06 is the first Monday of May 2024 target_date = datetime(2024, 5, 6, tzinfo=timezone.utc) @@ -379,7 +424,9 @@ async def test_monday_interpolation_skips_fwi_calculation(self, mock_dependencie mock_dependencies.fwi_processor.calculate_index.assert_not_called() @pytest.mark.anyio - async def test_regular_day_runs_fwi_calculation_not_interpolation(self, mock_dependencies: MockDailyActualsDeps): + async def test_regular_day_runs_fwi_calculation_not_interpolation( + self, mock_dependencies: MockDailyActualsDeps + ): """On a regular day, FWI calculation runs and FWI interpolation does not.""" # 2024-07-04 is a Thursday target_date = datetime(2024, 7, 4, tzinfo=timezone.utc) diff --git a/backend/packages/wps-api/src/app/tests/sfms/test_raster_addresser.py b/backend/packages/wps-api/src/app/tests/sfms/test_raster_addresser.py index fda4541f7..51b0c22c3 100644 --- a/backend/packages/wps-api/src/app/tests/sfms/test_raster_addresser.py +++ b/backend/packages/wps-api/src/app/tests/sfms/test_raster_addresser.py @@ -107,6 +107,10 @@ def test_get_calculated_hffmc_index_key(addresser: RasterKeyAddresser): SFMSInterpolatedWeatherParameter.WIND_SPEED, "sfms/interpolated/wind_speed/2024/01/15/wind_speed_20240115.tif", ), + ( + SFMSInterpolatedWeatherParameter.WIND_DIRECTION, + "sfms/interpolated/wind_direction/2024/01/15/wind_direction_20240115.tif", + ), ( SFMSInterpolatedWeatherParameter.PRECIP, "sfms/interpolated/precipitation/2024/01/15/precipitation_20240115.tif", diff --git a/backend/packages/wps-sfms/src/wps_sfms/interpolation/__init__.py b/backend/packages/wps-sfms/src/wps_sfms/interpolation/__init__.py index d7271ccfe..b7b37c250 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/interpolation/__init__.py +++ b/backend/packages/wps-sfms/src/wps_sfms/interpolation/__init__.py @@ -2,22 +2,42 @@ Interpolation modules for SFMS weather data. """ -from wps_sfms.interpolation.source import ( +from wps_sfms.interpolation.field import ( DEW_POINT_LAPSE_RATE, LAPSE_RATE, - StationActualSource, - StationInterpolationSource, - StationTemperatureSource, - StationPrecipitationSource, + ScalarField, + WindVectorField, + build_attribute_field, + build_dc_field, + build_dewpoint_field, + build_dmc_field, + build_ffmc_field, + build_precipitation_field, + build_temperature_field, + build_wind_speed_field, + build_wind_vector_field, + compute_adjusted_values, + compute_rh, + compute_sea_level_values, ) from wps_sfms.interpolation.common import log_interpolation_stats __all__ = [ "DEW_POINT_LAPSE_RATE", "LAPSE_RATE", - "StationActualSource", - "StationInterpolationSource", - "StationTemperatureSource", - "StationPrecipitationSource", + "ScalarField", + "WindVectorField", + "build_attribute_field", + "build_dc_field", + "build_dewpoint_field", + "build_dmc_field", + "build_ffmc_field", + "build_precipitation_field", + "build_temperature_field", + "build_wind_speed_field", + "build_wind_vector_field", + "compute_adjusted_values", + "compute_rh", + "compute_sea_level_values", "log_interpolation_stats", ] diff --git a/backend/packages/wps-sfms/src/wps_sfms/interpolation/field.py b/backend/packages/wps-sfms/src/wps_sfms/interpolation/field.py new file mode 100644 index 000000000..b5684d090 --- /dev/null +++ b/backend/packages/wps-sfms/src/wps_sfms/interpolation/field.py @@ -0,0 +1,158 @@ +from dataclasses import dataclass +import logging +from typing import List + +import numpy as np +from numpy.typing import NDArray +from wps_shared.schemas.sfms import SFMSDailyActual + +logger = logging.getLogger(__name__) + +# Environmental lapse rate: 6.5°C per 1000m elevation (average observed rate) +# This matches the CWFIS implementation. +LAPSE_RATE = 0.0065 + +# Dew point lapse rate: 2.0°C per 1000m. +DEW_POINT_LAPSE_RATE = 0.002 + +_VALID_SFMS_ATTRIBUTES = frozenset(SFMSDailyActual.model_fields.keys()) + + +@dataclass(frozen=True) +class ScalarField: + lats: NDArray[np.float32] + lons: NDArray[np.float32] + values: NDArray[np.float32] + + +@dataclass(frozen=True) +class WindVectorField: + lats: NDArray[np.float32] + lons: NDArray[np.float32] + u: NDArray[np.float32] + v: NDArray[np.float32] + + +def compute_sea_level_values( + values: NDArray[np.float32], elevs: NDArray[np.float32], lapse_rate: float +) -> NDArray[np.float32]: + """Vectorized: V_sea = V_station + elevation * lapse_rate""" + values = np.asarray(values, dtype=np.float32) + elevs = np.asarray(elevs, dtype=np.float32) + return values + elevs * np.float32(lapse_rate) + + +def compute_adjusted_values( + sea: NDArray[np.float32], elev: NDArray[np.float32], lapse_rate: float +) -> NDArray[np.float32]: + """Vectorized: V(z) = V_sea - z * lapse_rate""" + sea = np.asarray(sea, dtype=np.float32) + elev = np.asarray(elev, dtype=np.float32) + return sea - elev * np.float32(lapse_rate) + + +def compute_rh(temp: NDArray[np.float32], dewpoint: NDArray[np.float32]) -> NDArray[np.float32]: + """ + Compute relative humidity from temperature and dew point using the Arden Buck equation. + + Buck (1981): e_s(T) = 6.1121 * exp((18.678 - T/234.5) * (T / (257.14 + T))) + RH = 100 * e_s(Td) / e_s(T) + + :param temp: Temperature array in Celsius + :param dewpoint: Dew point temperature array in Celsius + :return: Relative humidity as percentage (0-100), clamped + """ + e_td = 6.1121 * np.exp((18.678 - dewpoint / 234.5) * (dewpoint / (257.14 + dewpoint))) + e_t = 6.1121 * np.exp((18.678 - temp / 234.5) * (temp / (257.14 + temp))) + rh = 100.0 * e_td / e_t + return np.clip(rh, 0.0, 100.0).astype(np.float32) + + +def build_attribute_field(actuals: List[SFMSDailyActual], attribute: str) -> ScalarField: + if attribute not in _VALID_SFMS_ATTRIBUTES: + raise ValueError( + f"Unknown attribute {attribute!r} on SFMSDailyActual. Valid attributes: {sorted(_VALID_SFMS_ATTRIBUTES)}" + ) + + valid = [s for s in actuals if getattr(s, attribute) is not None] + return ScalarField( + lats=np.array([s.lat for s in valid], dtype=np.float32), + lons=np.array([s.lon for s in valid], dtype=np.float32), + values=np.array([getattr(s, attribute) for s in valid], dtype=np.float32), + ) + + +def build_temperature_field(actuals: List[SFMSDailyActual]) -> ScalarField: + return _build_lapse_rate_field(actuals, "temperature", LAPSE_RATE, "temperature") + + +def build_dewpoint_field(actuals: List[SFMSDailyActual]) -> ScalarField: + return _build_lapse_rate_field(actuals, "dewpoint", DEW_POINT_LAPSE_RATE, "dew point") + + +def build_precipitation_field(actuals: List[SFMSDailyActual]) -> ScalarField: + return build_attribute_field(actuals, "precipitation") + + +def build_wind_speed_field(actuals: List[SFMSDailyActual]) -> ScalarField: + return build_attribute_field(actuals, "wind_speed") + + +def build_ffmc_field(actuals: List[SFMSDailyActual]) -> ScalarField: + return build_attribute_field(actuals, "ffmc") + + +def build_dmc_field(actuals: List[SFMSDailyActual]) -> ScalarField: + return build_attribute_field(actuals, "dmc") + + +def build_dc_field(actuals: List[SFMSDailyActual]) -> ScalarField: + return build_attribute_field(actuals, "dc") + + +def build_wind_vector_field(actuals: List[SFMSDailyActual]) -> WindVectorField: + valid = [s for s in actuals if s.wind_speed is not None and s.wind_direction is not None] + if not valid: + empty = np.array([], dtype=np.float32) + return WindVectorField(empty, empty, empty, empty) + + lats = np.array([s.lat for s in valid], dtype=np.float32) + lons = np.array([s.lon for s in valid], dtype=np.float32) + speed = np.array([s.wind_speed for s in valid], dtype=np.float32) + direction = np.array([s.wind_direction for s in valid], dtype=np.float32) + + direction_radians = np.radians(direction.astype(np.float32)) + u = (-speed * np.sin(direction_radians)).astype(np.float32) + v = (-speed * np.cos(direction_radians)).astype(np.float32) + return WindVectorField(lats=lats, lons=lons, u=u, v=v) + + +def _build_lapse_rate_field( + actuals: List[SFMSDailyActual], attribute: str, lapse_rate: float, label: str +) -> ScalarField: + if attribute not in _VALID_SFMS_ATTRIBUTES: + raise ValueError( + f"Unknown attribute {attribute!r} on SFMSDailyActual. Valid attributes: {sorted(_VALID_SFMS_ATTRIBUTES)}" + ) + lats = np.array([a.lat for a in actuals], dtype=np.float32) + lons = np.array([a.lon for a in actuals], dtype=np.float32) + elevs = np.array( + [a.elevation if a.elevation is not None else np.nan for a in actuals], dtype=np.float32 + ) + values = np.array( + [getattr(a, attribute) if getattr(a, attribute) is not None else np.nan for a in actuals], + dtype=np.float32, + ) + valid = np.isfinite(elevs) & np.isfinite(values) + + if not np.any(valid): + logger.warning( + "No valid %s stations (missing elevation or value) — interpolation will produce no output", + label, + ) + + return ScalarField( + lats=lats[valid], + lons=lons[valid], + values=compute_sea_level_values(values[valid], elevs[valid], lapse_rate), + ) diff --git a/backend/packages/wps-sfms/src/wps_sfms/interpolation/grid.py b/backend/packages/wps-sfms/src/wps_sfms/interpolation/grid.py new file mode 100644 index 000000000..ac4c4e634 --- /dev/null +++ b/backend/packages/wps-sfms/src/wps_sfms/interpolation/grid.py @@ -0,0 +1,122 @@ +"""Shared raster-grid loading for SFMS interpolation workflows. + +This module centralizes the repeated setup work needed before interpolation: +- open the reference raster and read its grid metadata +- apply the BC mask to determine which pixels are valid for interpolation +- precompute valid pixel coordinates in WGS84 +- optionally load auxiliary rasters that must align with the reference grid + +The goal is to let interpolation code focus on parameter-specific math while +this module owns the common raster plumbing and grid consistency checks. +""" + +from dataclasses import dataclass +from typing import Optional + +import numpy as np +from numpy.typing import NDArray +from wps_shared.geospatial.geospatial import rasters_match +from wps_shared.geospatial.wps_dataset import WPSDataset + + +@dataclass(frozen=True) +class GridContext: + """Shared raster grid state for interpolation workflows.""" + + geotransform: tuple[float, ...] + projection: str + x_size: int + y_size: int + valid_mask: NDArray[np.bool_] + valid_lats: NDArray[np.float32] + valid_lons: NDArray[np.float32] + valid_yi: NDArray[np.intp] + valid_xi: NDArray[np.intp] + total_pixels: int + skipped_nodata_count: int + valid_dem_values: Optional[NDArray[np.float32]] = None + temperature_data: Optional[NDArray[np.float32]] = None + + +def build_grid_context( + reference_raster_path: str, + mask_path: str, + *, + dem_path: Optional[str] = None, + temperature_raster_path: Optional[str] = None, +) -> GridContext: + """Load and validate the shared raster state for interpolation. + + The reference raster defines the output grid. This function applies the BC + mask to that grid, computes the valid pixel coordinates used by + IDW interpolation, and optionally loads aligned auxiliary rasters such as + DEM or temperature. + + The optional raster inputs are validated against the reference raster so + downstream interpolation code can assume the arrays are already on the same + grid and can be indexed directly with the valid pixel mask. + + :param reference_raster_path: Raster whose grid defines the interpolation output. + :param mask_path: BC mask raster; zero/nodata pixels are excluded from interpolation. + :param dem_path: Optional DEM raster aligned to the reference grid. + :param temperature_raster_path: Optional temperature raster aligned to the reference grid. + :return: A ``GridContext`` containing metadata, valid-pixel coordinates, and + any requested auxiliary raster data. + """ + with WPSDataset(reference_raster_path) as ref_ds: + geotransform = ref_ds.ds.GetGeoTransform() + if geotransform is None: + raise ValueError( + f"Failed to get geotransform from reference raster: {reference_raster_path}" + ) + projection = ref_ds.ds.GetProjection() + x_size = ref_ds.ds.RasterXSize + y_size = ref_ds.ds.RasterYSize + + with WPSDataset(mask_path) as mask_ds: + valid_mask = ref_ds.apply_mask(mask_ds) + + valid_lats, valid_lons, valid_yi, valid_xi = ref_ds.get_lat_lon_coords(valid_mask) + total_pixels = x_size * y_size + skipped_nodata_count = total_pixels - len(valid_yi) + + valid_dem_values = None + if dem_path is not None: + with WPSDataset(dem_path) as dem_ds: + _validate_matching_grid(ref_ds, dem_ds, "DEM") + dem_band = dem_ds.ds.GetRasterBand(1) + dem_data = dem_band.ReadAsArray() + if dem_data is None: + raise ValueError("Failed to read DEM data") + valid_dem_values = dem_data[valid_mask].astype(np.float32, copy=False) + + temperature_data = None + if temperature_raster_path is not None: + with WPSDataset(temperature_raster_path) as temp_ds: + _validate_matching_grid(ref_ds, temp_ds, "temperature raster") + temp_band = temp_ds.ds.GetRasterBand(1) + temp_data = temp_band.ReadAsArray() + if temp_data is None: + raise ValueError("Failed to read temperature raster data") + temperature_data = temp_data.astype(np.float32, copy=False) + + return GridContext( + geotransform=tuple(geotransform), + projection=projection, + x_size=x_size, + y_size=y_size, + valid_mask=valid_mask.astype(np.bool_, copy=False), + valid_lats=valid_lats.astype(np.float32, copy=False), + valid_lons=valid_lons.astype(np.float32, copy=False), + valid_yi=valid_yi.astype(np.intp, copy=False), + valid_xi=valid_xi.astype(np.intp, copy=False), + total_pixels=total_pixels, + skipped_nodata_count=skipped_nodata_count, + valid_dem_values=valid_dem_values, + temperature_data=temperature_data, + ) + + +def _validate_matching_grid(reference_ds: WPSDataset, candidate_ds: WPSDataset, label: str) -> None: + if not rasters_match(reference_ds.ds, candidate_ds.ds): + raise ValueError(f"{label} grid does not match reference raster") diff --git a/backend/packages/wps-sfms/src/wps_sfms/interpolation/idw.py b/backend/packages/wps-sfms/src/wps_sfms/interpolation/idw.py index b551868b8..ecbc7cd39 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/interpolation/idw.py +++ b/backend/packages/wps-sfms/src/wps_sfms/interpolation/idw.py @@ -14,6 +14,7 @@ SFMS_NO_DATA, log_interpolation_stats, ) +from wps_sfms.interpolation.grid import build_grid_context logger = logging.getLogger(__name__) @@ -37,70 +38,61 @@ def interpolate_to_raster( """ logger.info("Starting interpolation for %d stations", len(station_lats)) - with WPSDataset(reference_raster_path) as ref_ds: - geo_transform = ref_ds.ds.GetGeoTransform() - if geo_transform is None: - raise ValueError( - f"Failed to get geotransform from reference raster: {reference_raster_path}" - ) - projection = ref_ds.ds.GetProjection() - x_size = ref_ds.ds.RasterXSize - y_size = ref_ds.ds.RasterYSize - - # Use BC mask to determine valid pixels - with WPSDataset(mask_path) as mask_ds: - valid_mask = ref_ds.apply_mask(mask_ds.warp_to_match(ref_ds)) - - lats, lons, valid_yi, valid_xi = ref_ds.get_lat_lon_coords(valid_mask) - - total_pixels = x_size * y_size - skipped_nodata_count = total_pixels - len(valid_yi) - - logger.info("Interpolating for raster grid (%d x %d)", x_size, y_size) - logger.info( - "Processing %d valid pixels (skipping %d NoData pixels)", - len(valid_yi), - skipped_nodata_count, + grid = build_grid_context(reference_raster_path, mask_path) + + logger.info("Interpolating for raster grid (%d x %d)", grid.x_size, grid.y_size) + logger.info( + "Processing %d valid pixels (skipping %d NoData pixels)", + len(grid.valid_yi), + grid.skipped_nodata_count, + ) + + logger.info( + "Running batch IDW interpolation for %d pixels and %d stations", + len(grid.valid_lats), + len(station_lats), + ) + station_lats_array = np.array(station_lats) + station_lons_array = np.array(station_lons) + station_values_array = np.array(station_values) + + interpolated_values = idw_interpolation( + grid.valid_lats, + grid.valid_lons, + station_lats_array, + station_lons_array, + station_values_array, + ) + assert isinstance(interpolated_values, np.ndarray) + + output_array = np.full((grid.y_size, grid.x_size), SFMS_NO_DATA, dtype=np.float32) + + interpolation_succeeded = ~np.isnan(interpolated_values) + interpolated_count = int(np.sum(interpolation_succeeded)) + failed_interpolation_count = len(interpolated_values) - interpolated_count + + output_array[ + grid.valid_yi[interpolation_succeeded], grid.valid_xi[interpolation_succeeded] + ] = interpolated_values[interpolation_succeeded] + + log_interpolation_stats( + grid.total_pixels, + interpolated_count, + failed_interpolation_count, + grid.skipped_nodata_count, + ) + + if interpolated_count == 0: + raise RuntimeError( + f"No pixels were successfully interpolated from {len(station_lats)} station(s). " + "Check that station coordinates fall within the raster extent and that at least " + "one station has a valid value for this parameter." ) - logger.info( - "Running batch IDW interpolation for %d pixels and %d stations", - len(lats), - len(station_lats), - ) - station_lats_array = np.array(station_lats) - station_lons_array = np.array(station_lons) - station_values_array = np.array(station_values) - - interpolated_values = idw_interpolation( - lats, lons, station_lats_array, station_lons_array, station_values_array - ) - assert isinstance(interpolated_values, np.ndarray) - - output_array = np.full((y_size, x_size), SFMS_NO_DATA, dtype=np.float32) - - interpolation_succeeded = ~np.isnan(interpolated_values) - interpolated_count = int(np.sum(interpolation_succeeded)) - failed_interpolation_count = len(interpolated_values) - interpolated_count - - for idx in np.nonzero(interpolation_succeeded)[0]: - output_array[valid_yi[idx], valid_xi[idx]] = interpolated_values[idx] - - log_interpolation_stats( - total_pixels, interpolated_count, failed_interpolation_count, skipped_nodata_count - ) - - if interpolated_count == 0: - raise RuntimeError( - f"No pixels were successfully interpolated from {len(station_lats)} station(s). " - "Check that station coordinates fall within the raster extent and that at least " - "one station has a valid value for this parameter." - ) - - return WPSDataset.from_array( - array=output_array, - geotransform=geo_transform, - projection=projection, - nodata_value=SFMS_NO_DATA, - datatype=gdal.GDT_Float32, - ) + return WPSDataset.from_array( + array=output_array, + geotransform=grid.geotransform, + projection=grid.projection, + nodata_value=SFMS_NO_DATA, + datatype=gdal.GDT_Float32, + ) diff --git a/backend/packages/wps-sfms/src/wps_sfms/interpolation/source.py b/backend/packages/wps-sfms/src/wps_sfms/interpolation/source.py deleted file mode 100644 index 24ec4c87c..000000000 --- a/backend/packages/wps-sfms/src/wps_sfms/interpolation/source.py +++ /dev/null @@ -1,201 +0,0 @@ -from abc import ABC, abstractmethod -import logging -from typing import List, Optional, Protocol, Tuple, runtime_checkable -import numpy as np - -from numpy.typing import NDArray -from wps_shared.schemas.sfms import SFMSDailyActual - -logger = logging.getLogger(__name__) - -# Environmental lapse rate: 6.5°C per 1000m elevation (average observed rate) -# This matches the CWFIS implementation -LAPSE_RATE = 0.0065 - -# Dew point lapse rate: 2.0°C per 1000m: https://www.atmos.illinois.edu/~snodgrss/Airflow_over_mtn.html -DEW_POINT_LAPSE_RATE = 0.002 - - -@runtime_checkable -class StationInterpolationSource(Protocol): - def get_interpolation_data(self) -> Tuple[NDArray[np.float32], NDArray[np.float32], NDArray[np.float32]]: ... - - -class LapseRateAdjustedSource(ABC): - """ - Base class for station sources whose values require lapse-rate elevation - adjustment (e.g. temperature, dew point). - - Subclasses only need to implement ``_extract_values`` to pull the relevant - per-station values from the pre-built numpy arrays. - """ - - def __init__(self, sfms_actuals: List[SFMSDailyActual]): - self._sfms_actuals = sfms_actuals - self._lats: Optional[NDArray[np.float32]] = None - self._lons: Optional[NDArray[np.float32]] = None - self._elevs: Optional[NDArray[np.float32]] = None - self._values: Optional[NDArray[np.float32]] = None - self._valid_mask: Optional[NDArray[np.float32]] = None - - @abstractmethod - def _extract_values(self, actuals: List[SFMSDailyActual]) -> NDArray[np.float32]: - """Return a float32 array of values (one per station), ``np.nan`` where unavailable.""" - ... - - @staticmethod - def compute_sea_level_values( - values: NDArray[np.float32], elevs: NDArray[np.float32], lapse_rate: float - ) -> NDArray[np.float32]: - """Vectorized: V_sea = V_station + elevation * lapse_rate""" - values = np.asarray(values, dtype=np.float32) - elevs = np.asarray(elevs, dtype=np.float32) - return values + elevs * np.float32(lapse_rate) - - @staticmethod - def compute_adjusted_values( - sea: NDArray[np.float32], elev: NDArray[np.float32], lapse_rate: float - ) -> NDArray[np.float32]: - """Vectorized: V(z) = V_sea - z * lapse_rate""" - sea = np.asarray(sea, dtype=np.float32) - elev = np.asarray(elev, dtype=np.float32) - return sea - elev * np.float32(lapse_rate) - - def get_station_count(self) -> int: - return len(self._sfms_actuals) - - def get_station_arrays( - self, only_valid: bool = True - ) -> Tuple[ - NDArray[np.float32], - NDArray[np.float32], - NDArray[np.float32], - NDArray[np.float32], - ]: - """ - Returns (lats, lons, elevs, values) as NumPy arrays (float32). - If only_valid=True, filters out stations missing elevation or value. - """ - self._ensure_arrays() - if only_valid: - v = self._valid_mask - return self._lats[v], self._lons[v], self._elevs[v], self._values[v] - return self._lats, self._lons, self._elevs, self._values - - def get_interpolation_data( - self, lapse_rate: float = LAPSE_RATE - ) -> Tuple[NDArray[np.float32], NDArray[np.float32], NDArray[np.float32]]: - """Returns arrays for IDW: (lats, lons, sea_level_values), vectorized.""" - lats, lons, elevs, values = self.get_station_arrays(only_valid=True) - if lats.size == 0: - logger.warning( - "%s has no valid stations (missing elevation or value) — interpolation will produce no output", - type(self).__name__, - ) - return lats, lons, values - - sea = self.compute_sea_level_values(values, elevs, lapse_rate) - return lats, lons, sea - - # ------------------------------------------------------------------ - # Internals - # ------------------------------------------------------------------ - - def _ensure_arrays(self) -> None: - if self._lats is None: - self._materialize_arrays() - - @staticmethod - def _optional_to_array(actuals: List[SFMSDailyActual], attr: str) -> NDArray[np.float32]: - """Extract an optional float attribute from each actual into a float32 array (None → nan).""" - return np.array( - [getattr(a, attr) if getattr(a, attr) is not None else np.nan for a in actuals], - dtype=np.float32, - ) - - def _materialize_arrays(self) -> None: - """Pulls values from sfms_actuals into float32 NumPy arrays and computes valid mask.""" - actuals = self._sfms_actuals - self._lats = np.array([a.lat for a in actuals], dtype=np.float32) - self._lons = np.array([a.lon for a in actuals], dtype=np.float32) - self._elevs = self._optional_to_array(actuals, "elevation") - self._values = self._extract_values(actuals) - self._valid_mask = np.isfinite(self._elevs) & np.isfinite(self._values) - - -class StationTemperatureSource(LapseRateAdjustedSource): - """Station source for temperature values with lapse-rate elevation adjustment.""" - - def __init__(self, sfms_actuals: List[SFMSDailyActual]): - super().__init__(sfms_actuals) - - def _extract_values(self, actuals: List[SFMSDailyActual]) -> NDArray[np.float32]: - return self._optional_to_array(actuals, "temperature") - - -class StationDewPointSource(LapseRateAdjustedSource): - """ - Station source for dew-point values (computed from temperature + RH) - with lapse-rate elevation adjustment. - """ - - def __init__(self, sfms_actuals: List[SFMSDailyActual]): - super().__init__(sfms_actuals) - - def _extract_values(self, actuals: List[SFMSDailyActual]) -> NDArray[np.float32]: - dewpoints = self._optional_to_array(actuals, "dewpoint") - return dewpoints - - @staticmethod - def compute_rh(temp: np.ndarray, dewpoint: np.ndarray) -> np.ndarray: - """ - Compute relative humidity from temperature and dew point using the Arden Buck equation. - - Buck (1981): e_s(T) = 6.1121 * exp((18.678 - T/234.5) * (T / (257.14 + T))) - RH = 100 * e_s(Td) / e_s(T) - - :param temp: Temperature array in Celsius - :param dewpoint: Dew point temperature array in Celsius - :return: Relative humidity as percentage (0-100), clamped - """ - e_td = 6.1121 * np.exp((18.678 - dewpoint / 234.5) * (dewpoint / (257.14 + dewpoint))) - e_t = 6.1121 * np.exp((18.678 - temp / 234.5) * (temp / (257.14 + temp))) - rh = 100.0 * e_td / e_t - return np.clip(rh, 0.0, 100.0).astype(np.float32) - - -_VALID_SFMS_ATTRIBUTES = frozenset(SFMSDailyActual.model_fields.keys()) - - -class StationActualSource(StationInterpolationSource): - """Generic source for interpolating a named attribute from SFMSDailyActual.""" - - def __init__(self, attribute: str, sfms_actuals: List[SFMSDailyActual]): - if attribute not in _VALID_SFMS_ATTRIBUTES: - raise ValueError(f"Unknown attribute {attribute!r} on SFMSDailyActual. Valid attributes: {sorted(_VALID_SFMS_ATTRIBUTES)}") - self._attribute = attribute - self._sfms_actuals = sfms_actuals - - def get_interpolation_data(self) -> Tuple[NDArray[np.float32], NDArray[np.float32], NDArray[np.float32]]: - valid = [s for s in self._sfms_actuals if getattr(s, self._attribute) is not None] - return ( - np.array([s.lat for s in valid], dtype=np.float32), - np.array([s.lon for s in valid], dtype=np.float32), - np.array([getattr(s, self._attribute) for s in valid], dtype=np.float32), - ) - - -def StationPrecipitationSource(sfms_actuals: List[SFMSDailyActual]) -> StationActualSource: - return StationActualSource("precipitation", sfms_actuals) - - -def StationFFMCSource(sfms_actuals: List[SFMSDailyActual]) -> StationActualSource: - return StationActualSource("ffmc", sfms_actuals) - - -def StationDMCSource(sfms_actuals: List[SFMSDailyActual]) -> StationActualSource: - return StationActualSource("dmc", sfms_actuals) - - -def StationDCSource(sfms_actuals: List[SFMSDailyActual]) -> StationActualSource: - return StationActualSource("dc", sfms_actuals) diff --git a/backend/packages/wps-sfms/src/wps_sfms/processors/__init__.py b/backend/packages/wps-sfms/src/wps_sfms/processors/__init__.py index a3ff3e0d2..8c8d35e30 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/processors/__init__.py +++ b/backend/packages/wps-sfms/src/wps_sfms/processors/__init__.py @@ -2,15 +2,33 @@ Processor modules for SFMS interpolation workflows. """ -from wps_sfms.processors.fwi import FWIProcessor, FWIResult +from wps_sfms.processors.fwi import ( + BUICalculator, + DCCalculator, + DMCCalculator, + FFMCCalculator, + FWIFinalCalculator, + FWIProcessor, + FWIResult, + ISICalculator, +) from wps_sfms.processors.idw import Interpolator from wps_sfms.processors.relative_humidity import RHInterpolator from wps_sfms.processors.temperature import TemperatureInterpolator +from wps_sfms.processors.wind import WindDirectionInterpolator, WindSpeedInterpolator __all__ = [ "FWIProcessor", "FWIResult", + "FFMCCalculator", + "DMCCalculator", + "DCCalculator", + "ISICalculator", + "BUICalculator", + "FWIFinalCalculator", "Interpolator", "RHInterpolator", "TemperatureInterpolator", + "WindDirectionInterpolator", + "WindSpeedInterpolator", ] diff --git a/backend/packages/wps-sfms/src/wps_sfms/processors/fwi.py b/backend/packages/wps-sfms/src/wps_sfms/processors/fwi.py index 80642bb71..5e254a455 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/processors/fwi.py +++ b/backend/packages/wps-sfms/src/wps_sfms/processors/fwi.py @@ -1,32 +1,48 @@ """ -FWI processor for calculating FFMC, DMC, DC rasters from weather inputs. +FWI processor for calculating FWI index rasters from weather/index dependencies. -Accepts a FWIInputs dataclass that fully describes the input and output keys, -allowing the same processor to be used for both actuals and forecasts. +Accepts an FWIInputs dataclass that declares all dependency keys and output keys, +allowing each calculator to specify only the inputs it needs. """ +from dataclasses import dataclass import logging -import tempfile from abc import ABC, abstractmethod +from contextlib import contextmanager from datetime import datetime from time import perf_counter -from typing import Callable, ContextManager, List, NamedTuple, Tuple, cast +from typing import Callable, ContextManager, Iterator, List, Mapping, NamedTuple import numpy as np from osgeo import gdal -from wps_shared.fwi import vectorized_dc, vectorized_dmc, vectorized_ffmc -from wps_shared.geospatial.cog import generate_and_store_cog +from wps_shared.fwi import ( + vectorized_bui, + vectorized_dc, + vectorized_dmc, + vectorized_ffmc, + vectorized_fwi, + vectorized_isi, +) from wps_shared.geospatial.geospatial import rasters_match +from wps_sfms.publish import publish_dataset from wps_shared.geospatial.wps_dataset import WPSDataset -from wps_shared.sfms.raster_addresser import FWIParameter -from wps_sfms.sfmsng_raster_addresser import FWIInputs +from wps_shared.sfms.raster_addresser import FWIParameter, SFMSInterpolatedWeatherParameter from wps_shared.utils.s3 import set_s3_gdal_config from wps_shared.utils.s3_client import S3Client +from wps_sfms.sfmsng_raster_addresser import FWIInputs logger = logging.getLogger(__name__) MultiDatasetContext = Callable[[List[str]], ContextManager[List["WPSDataset"]]] +WeatherDatasetMap = dict[SFMSInterpolatedWeatherParameter, WPSDataset] +IndexDatasetMap = dict[FWIParameter, WPSDataset] + + +@dataclass +class FWIDatasets: + index: IndexDatasetMap + weather: WeatherDatasetMap class FWIResult(NamedTuple): @@ -36,11 +52,12 @@ class FWIResult(NamedTuple): class FWICalculator(ABC): fwi_param: FWIParameter + reference_index_param: FWIParameter + required_weather_params: tuple[SFMSInterpolatedWeatherParameter, ...] = () + required_index_params: tuple[FWIParameter, ...] = () @abstractmethod - def calculate( - self, prev_fwi_ds: WPSDataset, temp_ds: WPSDataset, rh_ds: WPSDataset, precip_ds: WPSDataset - ) -> FWIResult: ... + def calculate(self, datasets: FWIDatasets) -> FWIResult: ... class MonthlyFWICalculator(FWICalculator, ABC): @@ -51,29 +68,43 @@ def __init__(self, month: int): raise ValueError(f"month must be 1–12, got {month}") self.month = month - def _lat_month_arrays(self, ds: WPSDataset) -> Tuple[np.ndarray, np.ndarray]: + def _lat_month_arrays(self, ds: WPSDataset) -> tuple[np.ndarray, np.ndarray]: latitude = ds.generate_latitude_array() return latitude, np.full(latitude.shape, self.month) class FFMCCalculator(FWICalculator): fwi_param = FWIParameter.FFMC - - def calculate( - self, prev_fwi_ds: WPSDataset, temp_ds: WPSDataset, rh_ds: WPSDataset, precip_ds: WPSDataset - ) -> FWIResult: - ffmc_yda, nodata_value = prev_fwi_ds.replace_nodata_with(np.nan) + reference_index_param = FWIParameter.FFMC + required_weather_params = ( + SFMSInterpolatedWeatherParameter.TEMP, + SFMSInterpolatedWeatherParameter.RH, + SFMSInterpolatedWeatherParameter.PRECIP, + SFMSInterpolatedWeatherParameter.WIND_SPEED, + ) + required_index_params = (FWIParameter.FFMC,) + + def calculate(self, datasets: FWIDatasets) -> FWIResult: + temp_ds = datasets.weather[SFMSInterpolatedWeatherParameter.TEMP] + rh_ds = datasets.weather[SFMSInterpolatedWeatherParameter.RH] + precip_ds = datasets.weather[SFMSInterpolatedWeatherParameter.PRECIP] + wind_ds = datasets.weather[SFMSInterpolatedWeatherParameter.WIND_SPEED] + ffmc_prev_ds = datasets.index[FWIParameter.FFMC] + + ffmc_prev, nodata_value = ffmc_prev_ds.replace_nodata_with(np.nan) temp, _ = temp_ds.replace_nodata_with(np.nan) rh, _ = rh_ds.replace_nodata_with(np.nan) prec, _ = precip_ds.replace_nodata_with(np.nan) - ws = np.zeros_like(temp) # TODO: implement wind speed + ws, _ = wind_ds.replace_nodata_with(np.nan) - mask = np.isnan(ffmc_yda) | np.isnan(temp) | np.isnan(rh) | np.isnan(prec) + mask = np.isnan(ffmc_prev) | np.isnan(temp) | np.isnan(rh) | np.isnan(prec) | np.isnan(ws) valid = ~mask - values = np.full(ffmc_yda.shape, nodata_value, dtype=ffmc_yda.dtype) + values = np.full(ffmc_prev.shape, nodata_value, dtype=ffmc_prev.dtype) start = perf_counter() - values[valid] = vectorized_ffmc(ffmc_yda[valid], temp[valid], rh[valid], ws[valid], prec[valid]) + values[valid] = vectorized_ffmc( + ffmc_prev[valid], temp[valid], rh[valid], ws[valid], prec[valid] + ) logger.info("%f seconds to calculate vectorized ffmc", perf_counter() - start) return FWIResult(values, nodata_value) @@ -81,22 +112,34 @@ def calculate( class DMCCalculator(MonthlyFWICalculator): fwi_param = FWIParameter.DMC - - def calculate( - self, prev_fwi_ds: WPSDataset, temp_ds: WPSDataset, rh_ds: WPSDataset, precip_ds: WPSDataset - ) -> FWIResult: - lat, mon = self._lat_month_arrays(prev_fwi_ds) - dmc_yda, nodata_value = prev_fwi_ds.replace_nodata_with(np.nan) + reference_index_param = FWIParameter.DMC + required_weather_params = ( + SFMSInterpolatedWeatherParameter.TEMP, + SFMSInterpolatedWeatherParameter.RH, + SFMSInterpolatedWeatherParameter.PRECIP, + ) + required_index_params = (FWIParameter.DMC,) + + def calculate(self, datasets: FWIDatasets) -> FWIResult: + temp_ds = datasets.weather[SFMSInterpolatedWeatherParameter.TEMP] + rh_ds = datasets.weather[SFMSInterpolatedWeatherParameter.RH] + precip_ds = datasets.weather[SFMSInterpolatedWeatherParameter.PRECIP] + dmc_prev_ds = datasets.index[FWIParameter.DMC] + + lat, mon = self._lat_month_arrays(dmc_prev_ds) + dmc_prev, nodata_value = dmc_prev_ds.replace_nodata_with(np.nan) temp, _ = temp_ds.replace_nodata_with(np.nan) rh, _ = rh_ds.replace_nodata_with(np.nan) prec, _ = precip_ds.replace_nodata_with(np.nan) - mask = np.isnan(dmc_yda) | np.isnan(temp) | np.isnan(rh) | np.isnan(prec) + mask = np.isnan(dmc_prev) | np.isnan(temp) | np.isnan(rh) | np.isnan(prec) valid = ~mask - values = np.full(dmc_yda.shape, nodata_value, dtype=dmc_yda.dtype) + values = np.full(dmc_prev.shape, nodata_value, dtype=dmc_prev.dtype) start = perf_counter() - values[valid] = vectorized_dmc(dmc_yda[valid], temp[valid], rh[valid], prec[valid], lat[valid], mon[valid], True) + values[valid] = vectorized_dmc( + dmc_prev[valid], temp[valid], rh[valid], prec[valid], lat[valid], mon[valid], True + ) logger.info("%f seconds to calculate vectorized dmc", perf_counter() - start) return FWIResult(values, nodata_value) @@ -104,33 +147,194 @@ def calculate( class DCCalculator(MonthlyFWICalculator): fwi_param = FWIParameter.DC - - def calculate( - self, prev_fwi_ds: WPSDataset, temp_ds: WPSDataset, rh_ds: WPSDataset, precip_ds: WPSDataset - ) -> FWIResult: - lat, mon = self._lat_month_arrays(prev_fwi_ds) - dc_yda, nodata_value = prev_fwi_ds.replace_nodata_with(np.nan) + reference_index_param = FWIParameter.DC + required_weather_params = ( + SFMSInterpolatedWeatherParameter.TEMP, + SFMSInterpolatedWeatherParameter.RH, + SFMSInterpolatedWeatherParameter.PRECIP, + ) + required_index_params = (FWIParameter.DC,) + + def calculate(self, datasets: FWIDatasets) -> FWIResult: + temp_ds = datasets.weather[SFMSInterpolatedWeatherParameter.TEMP] + rh_ds = datasets.weather[SFMSInterpolatedWeatherParameter.RH] + precip_ds = datasets.weather[SFMSInterpolatedWeatherParameter.PRECIP] + dc_prev_ds = datasets.index[FWIParameter.DC] + + lat, mon = self._lat_month_arrays(dc_prev_ds) + dc_prev, nodata_value = dc_prev_ds.replace_nodata_with(np.nan) temp, _ = temp_ds.replace_nodata_with(np.nan) rh, _ = rh_ds.replace_nodata_with(np.nan) prec, _ = precip_ds.replace_nodata_with(np.nan) - mask = np.isnan(dc_yda) | np.isnan(temp) | np.isnan(rh) | np.isnan(prec) + mask = np.isnan(dc_prev) | np.isnan(temp) | np.isnan(rh) | np.isnan(prec) valid = ~mask - values = np.full(dc_yda.shape, nodata_value, dtype=dc_yda.dtype) + values = np.full(dc_prev.shape, nodata_value, dtype=dc_prev.dtype) start = perf_counter() - values[valid] = vectorized_dc(dc_yda[valid], temp[valid], rh[valid], prec[valid], lat[valid], mon[valid], True) + values[valid] = vectorized_dc( + dc_prev[valid], temp[valid], rh[valid], prec[valid], lat[valid], mon[valid], True + ) logger.info("%f seconds to calculate vectorized dc", perf_counter() - start) return FWIResult(values, nodata_value) +class ISICalculator(FWICalculator): + fwi_param = FWIParameter.ISI + reference_index_param = FWIParameter.FFMC + required_weather_params = (SFMSInterpolatedWeatherParameter.WIND_SPEED,) + required_index_params = (FWIParameter.FFMC,) + + def calculate(self, datasets: FWIDatasets) -> FWIResult: + ffmc_ds = datasets.index[FWIParameter.FFMC] + wind_ds = datasets.weather[SFMSInterpolatedWeatherParameter.WIND_SPEED] + + ffmc, nodata_value = ffmc_ds.replace_nodata_with(np.nan) + ws, _ = wind_ds.replace_nodata_with(np.nan) + + mask = np.isnan(ffmc) | np.isnan(ws) + valid = ~mask + values = np.full(ffmc.shape, nodata_value, dtype=ffmc.dtype) + + start = perf_counter() + values[valid] = vectorized_isi(ffmc[valid], ws[valid], False) + logger.info("%f seconds to calculate vectorized isi", perf_counter() - start) + + return FWIResult(values, nodata_value) + + +class BUICalculator(FWICalculator): + fwi_param = FWIParameter.BUI + reference_index_param = FWIParameter.DMC + required_index_params = (FWIParameter.DMC, FWIParameter.DC) + + def calculate(self, datasets: FWIDatasets) -> FWIResult: + dmc_ds = datasets.index[FWIParameter.DMC] + dc_ds = datasets.index[FWIParameter.DC] + + dmc, nodata_value = dmc_ds.replace_nodata_with(np.nan) + dc, _ = dc_ds.replace_nodata_with(np.nan) + + mask = np.isnan(dmc) | np.isnan(dc) + valid = ~mask + values = np.full(dmc.shape, nodata_value, dtype=dmc.dtype) + + start = perf_counter() + values[valid] = vectorized_bui(dmc[valid], dc[valid]) + logger.info("%f seconds to calculate vectorized bui", perf_counter() - start) + + return FWIResult(values, nodata_value) + + +class FWIFinalCalculator(FWICalculator): + fwi_param = FWIParameter.FWI + reference_index_param = FWIParameter.ISI + required_index_params = (FWIParameter.ISI, FWIParameter.BUI) + + def calculate(self, datasets: FWIDatasets) -> FWIResult: + isi_ds = datasets.index[FWIParameter.ISI] + bui_ds = datasets.index[FWIParameter.BUI] + + isi, nodata_value = isi_ds.replace_nodata_with(np.nan) + bui, _ = bui_ds.replace_nodata_with(np.nan) + + mask = np.isnan(isi) | np.isnan(bui) + valid = ~mask + values = np.full(isi.shape, nodata_value, dtype=isi.dtype) + + start = perf_counter() + values[valid] = vectorized_fwi(isi[valid], bui[valid]) + logger.info("%f seconds to calculate vectorized fwi", perf_counter() - start) + + return FWIResult(values, nodata_value) + + class FWIProcessor: - """Calculates FFMC, DMC, DC rasters from weather inputs described by FWIInputs.""" + """Calculates FWI index rasters from dependency keys described by FWIInputs.""" def __init__(self, datetime_to_process: datetime): self.datetime_to_process = datetime_to_process + @staticmethod + def _get_required_weather_keys( + calculator: FWICalculator, + fwi_inputs: FWIInputs, + ) -> dict[SFMSInterpolatedWeatherParameter, str]: + missing_params = [ + param + for param in calculator.required_weather_params + if param not in fwi_inputs.weather_keys + ] + if missing_params: + missing = ", ".join(param.value for param in missing_params) + raise ValueError(f"FWIInputs missing weather key mappings for: {missing}") + + return { + param: fwi_inputs.weather_keys[param] for param in calculator.required_weather_params + } + + @staticmethod + def _get_required_index_keys( + calculator: FWICalculator, + fwi_inputs: FWIInputs, + ) -> dict[FWIParameter, str]: + missing_params = [ + param + for param in calculator.required_index_params + if param not in fwi_inputs.index_keys + ] + if missing_params: + missing = ", ".join(param.value for param in missing_params) + raise ValueError(f"FWIInputs missing index key mappings for: {missing}") + + return {param: fwi_inputs.index_keys[param] for param in calculator.required_index_params} + + async def _assert_dependency_keys_exist( + self, + s3_client: S3Client, + keys_by_param: Mapping[SFMSInterpolatedWeatherParameter | FWIParameter, str], + dependency_kind: str, + ) -> None: + if not keys_by_param: + return + + if await s3_client.all_objects_exist(*keys_by_param.values()): + return + + details = ", ".join(f"{param.value}={key}" for param, key in keys_by_param.items()) + raise RuntimeError( + f"Missing {dependency_kind} keys for {self.datetime_to_process.date()}: {details}" + ) + + @contextmanager + def _open_required_datasets( + self, + input_dataset_context: MultiDatasetContext, + weather_keys_by_param: Mapping[SFMSInterpolatedWeatherParameter, str], + index_keys_by_param: Mapping[FWIParameter, str], + ) -> Iterator[FWIDatasets]: + """ + Open all required weather/index dependency rasters and yield them as `FWIDatasets`. + + This method opens every key in one context manager and reorganizes the resulting + datasets back into parameter-keyed weather and index maps for calculator consumption. + """ + input_keys = [*weather_keys_by_param.values(), *index_keys_by_param.values()] + with input_dataset_context(input_keys) as input_datasets: + datasets_by_key: dict[str, WPSDataset] = {} + for ds in input_datasets: + datasets_by_key[ds.ds_path] = ds + + weather_datasets: WeatherDatasetMap = { + param: datasets_by_key[key] for param, key in weather_keys_by_param.items() + } + index_datasets: IndexDatasetMap = { + param: datasets_by_key[key] for param, key in index_keys_by_param.items() + } + + yield FWIDatasets(index=index_datasets, weather=weather_datasets) + async def calculate_index( self, s3_client: S3Client, @@ -139,30 +343,24 @@ async def calculate_index( fwi_inputs: FWIInputs, ): """ - Calculate a single FWI index from the provided inputs. + Calculate a single FWI index from the provided dependencies. :param s3_client: S3Client instance for checking keys and persisting results - :param input_dataset_context: Context manager for opening multiple WPSDatasets + :param input_dataset_context: Context manager for opening dependency datasets :param calculator: FWICalculator instance that performs the index calculation - :param fwi_inputs: All S3 keys and metadata for this calculation + :param fwi_inputs: Dependency keys and metadata for this calculation """ set_s3_gdal_config() - weather_keys_exist = await s3_client.all_objects_exist( - fwi_inputs.temp_key, fwi_inputs.rh_key, fwi_inputs.precip_key - ) - if not weather_keys_exist: - raise RuntimeError( - f"Missing weather keys for {self.datetime_to_process.date()}: " - f"temp={fwi_inputs.temp_key}, rh={fwi_inputs.rh_key}, precip={fwi_inputs.precip_key}" - ) + # get only the dependency keys required by this calculator. + # ex: ISI uses FFMC + wind speed, while BUI uses DMC + DC. + weather_keys_by_param = self._get_required_weather_keys(calculator, fwi_inputs) + index_keys_by_param = self._get_required_index_keys(calculator, fwi_inputs) - fwi_key_exists = await s3_client.all_objects_exist(fwi_inputs.prev_fwi_key) - if not fwi_key_exists: - raise RuntimeError( - f"Missing previous {calculator.fwi_param.value} raster for " - f"{self.datetime_to_process.date()}: {fwi_inputs.prev_fwi_key}" - ) + await self._assert_dependency_keys_exist( + s3_client, weather_keys_by_param, "weather dependency" + ) + await self._assert_dependency_keys_exist(s3_client, index_keys_by_param, "index dependency") logger.info( "Calculating %s %s for %s", @@ -171,56 +369,57 @@ async def calculate_index( self.datetime_to_process.date(), ) - with tempfile.TemporaryDirectory() as temp_dir: - with input_dataset_context( - [ - fwi_inputs.temp_key, - fwi_inputs.rh_key, - fwi_inputs.precip_key, - fwi_inputs.prev_fwi_key, - ] - ) as input_datasets: - input_datasets = cast(List[WPSDataset], input_datasets) - temp_ds, rh_ds, precip_ds, prev_fwi_ds = input_datasets - - # Assert weather rasters already match the FWI grid - if not rasters_match(temp_ds.as_gdal_ds(), prev_fwi_ds.as_gdal_ds()): - raise ValueError( - f"Temperature raster does not match FWI grid: {fwi_inputs.temp_key} vs {fwi_inputs.prev_fwi_key}" - ) - if not rasters_match(rh_ds.as_gdal_ds(), prev_fwi_ds.as_gdal_ds()): + with self._open_required_datasets( + input_dataset_context, weather_keys_by_param, index_keys_by_param + ) as datasets: + weather_datasets = datasets.weather + index_datasets = datasets.index + + # use reference index's geotransform and projection for the output dataset, and verify all dependencies match that grid + reference_ds = index_datasets[calculator.reference_index_param] + reference_key = index_keys_by_param[calculator.reference_index_param] + + # every weather raster must match the reference index grid + for param in calculator.required_weather_params: + weather_ds = weather_datasets[param] + weather_key = weather_keys_by_param[param] + if not rasters_match(weather_ds.as_gdal_ds(), reference_ds.as_gdal_ds()): raise ValueError( - f"RH raster does not match FWI grid: {fwi_inputs.rh_key} vs {fwi_inputs.prev_fwi_key}" + f"{param.value} raster does not match FWI grid: {weather_key} vs {reference_key}" ) - if not rasters_match(precip_ds.as_gdal_ds(), prev_fwi_ds.as_gdal_ds()): + + # every index raster must match the reference index grid + for param in calculator.required_index_params: + if param == calculator.reference_index_param: + continue + index_ds = index_datasets[param] + index_key = index_keys_by_param[param] + if not rasters_match(index_ds.as_gdal_ds(), reference_ds.as_gdal_ds()): raise ValueError( - f"Precip raster does not match FWI grid: {fwi_inputs.precip_key} vs {fwi_inputs.prev_fwi_key}" + f"{param.value} raster does not match FWI grid: {index_key} vs {reference_key}" ) - result = calculator.calculate(prev_fwi_ds, temp_ds, rh_ds, precip_ds) - - await s3_client.persist_raster_data( - temp_dir, - fwi_inputs.output_key, - prev_fwi_ds.as_gdal_ds().GetGeoTransform(), - prev_fwi_ds.as_gdal_ds().GetProjection(), - result.values, - result.nodata_value, + result = calculator.calculate(datasets) + + with WPSDataset.from_array( + result.values, + reference_ds.as_gdal_ds().GetGeoTransform(), + reference_ds.as_gdal_ds().GetProjection(), + result.nodata_value, + ) as output_ds: + published = await publish_dataset( + s3_client=s3_client, + dataset=output_ds, + output_key=fwi_inputs.output_key, ) - with WPSDataset.from_array( - result.values, - prev_fwi_ds.as_gdal_ds().GetGeoTransform(), - prev_fwi_ds.as_gdal_ds().GetProjection(), - result.nodata_value, - ) as output_ds: - generate_and_store_cog(src_ds=output_ds.as_gdal_ds(), output_path=fwi_inputs.cog_key) - logger.info( - "Stored %s %s: %s", - calculator.fwi_param.value, - fwi_inputs.run_type.value, - fwi_inputs.output_key, - ) + logger.info( + "Stored %s %s: %s (COG: %s)", + calculator.fwi_param.value, + fwi_inputs.run_type.value, + published.output_key, + published.cog_key, + ) - # Clear gdal virtual file system cache of S3 metadata - gdal.VSICurlClearCache() + # Clear gdal virtual file system cache of S3 metadata. + gdal.VSICurlClearCache() diff --git a/backend/packages/wps-sfms/src/wps_sfms/processors/idw.py b/backend/packages/wps-sfms/src/wps_sfms/processors/idw.py index 74e48b956..28d0f03be 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/processors/idw.py +++ b/backend/packages/wps-sfms/src/wps_sfms/processors/idw.py @@ -10,40 +10,104 @@ """ import logging -import os -import tempfile +from abc import ABC, abstractmethod +from dataclasses import dataclass -import aiofiles +import numpy as np +from numpy.typing import NDArray +from wps_sfms.publish import publish_dataset from wps_shared.geospatial.wps_dataset import WPSDataset +from wps_shared.geospatial.spatial_interpolation import idw_interpolation from wps_shared.utils.s3 import set_s3_gdal_config from wps_shared.utils.s3_client import S3Client +from wps_sfms.interpolation.field import ScalarField from wps_sfms.interpolation.idw import interpolate_to_raster -from wps_sfms.interpolation.source import StationInterpolationSource logger = logging.getLogger(__name__) -class Interpolator: - """Base processor: plain IDW interpolation + S3 upload. +@dataclass(frozen=True) +class ValidPixelIDWResult: + """Result of IDW interpolation performed on valid raster pixels only.""" - Subclasses override interpolate() to add parameter-specific logic - such as elevation adjustment or derived quantities. - """ + interpolated_values: NDArray[np.float32] + succeeded_mask: NDArray[np.bool_] + rows: NDArray[np.intp] + cols: NDArray[np.intp] + values: NDArray[np.float32] + total_pixels: int + interpolated_count: int + failed_interpolation_count: int + skipped_nodata_count: int + + +def idw_on_valid_pixels( + valid_lats: NDArray[np.float32], + valid_lons: NDArray[np.float32], + valid_yi: NDArray[np.intp], + valid_xi: NDArray[np.intp], + station_lats: NDArray[np.float32], + station_lons: NDArray[np.float32], + station_values: NDArray[np.float32], + total_pixels: int, + label: str, +) -> ValidPixelIDWResult: + """Run batch IDW interpolation for valid pixels and return indexed results.""" + skipped_nodata_count = total_pixels - len(valid_yi) + + logger.info( + "Processing %d valid pixels (skipping %d NoData pixels)", + len(valid_yi), + skipped_nodata_count, + ) + logger.info( + "Running batch %s IDW interpolation for %d pixels and %d stations", + label, + len(valid_lats), + len(station_lats), + ) + + raw_interpolated = idw_interpolation( + valid_lats, valid_lons, station_lats, station_lons, station_values + ) + assert isinstance(raw_interpolated, np.ndarray) + interpolated_values = raw_interpolated.astype(np.float32, copy=False) + + succeeded_mask = ~np.isnan(interpolated_values) + interpolated_count = int(np.sum(succeeded_mask)) + failed_interpolation_count = len(interpolated_values) - interpolated_count + + rows = valid_yi[succeeded_mask] + cols = valid_xi[succeeded_mask] + values = interpolated_values[succeeded_mask].astype(np.float32, copy=False) + + return ValidPixelIDWResult( + interpolated_values=interpolated_values, + succeeded_mask=succeeded_mask, + rows=rows, + cols=cols, + values=values, + total_pixels=total_pixels, + interpolated_count=interpolated_count, + failed_interpolation_count=failed_interpolation_count, + skipped_nodata_count=skipped_nodata_count, + ) + + +class RasterProcessor(ABC): + """Shared upload workflow for interpolation processors.""" def __init__(self, mask_path: str): self.mask_path = mask_path - def interpolate( - self, source: StationInterpolationSource, reference_raster_path: str - ) -> WPSDataset: - lats, lons, values = source.get_interpolation_data() - return interpolate_to_raster(lats, lons, values, reference_raster_path, self.mask_path) + @abstractmethod + def interpolate(self, reference_raster_path: str) -> WPSDataset: + """Build an in-memory raster for upload.""" async def process( self, s3_client: S3Client, reference_raster_path: str, - source: StationInterpolationSource, output_key: str, ) -> str: """ @@ -51,20 +115,35 @@ async def process( :param s3_client: S3Client instance for uploading results :param reference_raster_path: Path to reference raster (defines grid properties) - :param source: Station data source providing interpolation inputs :param output_key: S3 key where the resulting raster will be uploaded :return: S3 key of uploaded raster """ set_s3_gdal_config() logger.info("Starting interpolation, output: %s", output_key) - with self.interpolate(source, reference_raster_path) as dataset: - with tempfile.TemporaryDirectory() as tmp_dir: - tmp_path = os.path.join(tmp_dir, os.path.basename(output_key)) - dataset.export_to_geotiff(tmp_path) + with self.interpolate(reference_raster_path) as dataset: + published = await publish_dataset(s3_client, dataset, output_key) - async with aiofiles.open(tmp_path, "rb") as f: - await s3_client.put_object(key=output_key, body=await f.read()) - - logger.info("Interpolation complete: %s", output_key) + logger.info( + "Interpolation complete: %s (COG: %s)", + published.output_key, + published.cog_key, + ) return output_key + + +class Interpolator(RasterProcessor): + """Scalar-field IDW interpolation plus shared upload workflow.""" + + def __init__(self, mask_path: str, field: ScalarField): + super().__init__(mask_path) + self.field = field + + def interpolate(self, reference_raster_path: str) -> WPSDataset: + return interpolate_to_raster( + self.field.lats, + self.field.lons, + self.field.values, + reference_raster_path, + self.mask_path, + ) diff --git a/backend/packages/wps-sfms/src/wps_sfms/processors/relative_humidity.py b/backend/packages/wps-sfms/src/wps_sfms/processors/relative_humidity.py index 5783a50b0..a2164929f 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/processors/relative_humidity.py +++ b/backend/packages/wps-sfms/src/wps_sfms/processors/relative_humidity.py @@ -1,127 +1,98 @@ import logging import numpy as np from osgeo import gdal -from wps_sfms.interpolation.source import DEW_POINT_LAPSE_RATE, StationDewPointSource +from wps_sfms.interpolation.field import ( + DEW_POINT_LAPSE_RATE, + ScalarField, + compute_adjusted_values, + compute_rh, +) from wps_shared.geospatial.wps_dataset import WPSDataset -from wps_shared.geospatial.spatial_interpolation import idw_interpolation from wps_sfms.interpolation.common import ( SFMS_NO_DATA, log_interpolation_stats, ) -from wps_sfms.processors.idw import Interpolator +from wps_sfms.interpolation.grid import build_grid_context +from wps_sfms.processors.idw import RasterProcessor, idw_on_valid_pixels logger = logging.getLogger(__name__) -class RHInterpolator(Interpolator): +class RHInterpolator(RasterProcessor): """Interpolates RH via dew point IDW + elevation adjustment. Requires that temperature interpolation has already been run for this date, as it reads the interpolated temperature raster from S3. """ - def __init__(self, mask_path: str, dem_path: str, temp_raster_path: str): + def __init__( + self, + mask_path: str, + dem_path: str, + temp_raster_path: str, + field: ScalarField, + ): super().__init__(mask_path) self.dem_path = dem_path self.temp_raster_path = temp_raster_path - - def interpolate( - self, source: StationDewPointSource, reference_raster_path: str - ) -> WPSDataset: - with WPSDataset(reference_raster_path) as ref_ds: - geo_transform = ref_ds.ds.GetGeoTransform() - if geo_transform is None: - raise ValueError( - f"Failed to get geotransform from reference raster: {reference_raster_path}" - ) - projection = ref_ds.ds.GetProjection() - x_size = ref_ds.ds.RasterXSize - y_size = ref_ds.ds.RasterYSize - - with WPSDataset(self.dem_path) as dem_ds: - dem_band: gdal.Band = dem_ds.ds.GetRasterBand(1) - dem_data = dem_band.ReadAsArray() - if dem_data is None: - raise ValueError("Failed to read DEM data") - - with WPSDataset(self.temp_raster_path) as temp_ds: - temp_band: gdal.Band = temp_ds.ds.GetRasterBand(1) - temp_data = temp_band.ReadAsArray() - if temp_data is None: - raise ValueError("Failed to read temperature raster data") - - rh_array = np.full((y_size, x_size), SFMS_NO_DATA, dtype=np.float32) - - with WPSDataset(self.mask_path) as mask_ds: - valid_mask = ref_ds.apply_mask(mask_ds.warp_to_match(ref_ds)) - - lats, lons, valid_yi, valid_xi = dem_ds.get_lat_lon_coords(valid_mask) - valid_elevations = dem_data[valid_mask] - - total_pixels = x_size * y_size - skipped_nodata_count = total_pixels - len(valid_yi) - - logger.info( - "Interpolating dew point for RH raster grid (%d x %d)", x_size, y_size - ) - logger.info( - "Processing %d valid pixels (skipping %d NoData pixels)", - len(valid_yi), - skipped_nodata_count, - ) - - station_lats, station_lons, sea_level_dewpoints = ( - source.get_interpolation_data(lapse_rate=DEW_POINT_LAPSE_RATE) - ) - logger.info( - "Running batch dew point IDW interpolation for %d pixels and %d stations", - len(lats), - len(station_lats), - ) - - interpolated_sea_level_dewpoints = idw_interpolation( - lats, lons, station_lats, station_lons, sea_level_dewpoints - ) - assert isinstance(interpolated_sea_level_dewpoints, np.ndarray) - - interpolation_succeeded = ~np.isnan(interpolated_sea_level_dewpoints) - interpolated_count = int(np.sum(interpolation_succeeded)) - failed_interpolation_count = ( - len(interpolated_sea_level_dewpoints) - interpolated_count - ) - - rows = valid_yi[interpolation_succeeded] - cols = valid_xi[interpolation_succeeded] - - sea = interpolated_sea_level_dewpoints[interpolation_succeeded].astype( - np.float32, copy=False - ) - elev = valid_elevations[interpolation_succeeded].astype(np.float32, copy=False) - adjusted_dewpoints = StationDewPointSource.compute_adjusted_values( - sea, elev, DEW_POINT_LAPSE_RATE - ) - - rh_values = StationDewPointSource.compute_rh( - temp_data[rows, cols].astype(np.float32), - adjusted_dewpoints, - ) - rh_array[rows, cols] = rh_values - - log_interpolation_stats( - total_pixels, interpolated_count, failed_interpolation_count, skipped_nodata_count + self.field = field + + def interpolate(self, reference_raster_path: str) -> WPSDataset: + grid = build_grid_context( + reference_raster_path, + self.mask_path, + dem_path=self.dem_path, + temperature_raster_path=self.temp_raster_path, + ) + assert grid.valid_dem_values is not None + assert grid.temperature_data is not None + + rh_array = np.full((grid.y_size, grid.x_size), SFMS_NO_DATA, dtype=np.float32) + + logger.info( + "Interpolating dew point for RH raster grid (%d x %d)", grid.x_size, grid.y_size + ) + + idw_result = idw_on_valid_pixels( + valid_lats=grid.valid_lats, + valid_lons=grid.valid_lons, + valid_yi=grid.valid_yi, + valid_xi=grid.valid_xi, + station_lats=self.field.lats, + station_lons=self.field.lons, + station_values=self.field.values, + total_pixels=grid.total_pixels, + label="dew point", + ) + + sea = idw_result.values + elev = grid.valid_dem_values[idw_result.succeeded_mask].astype(np.float32, copy=False) + adjusted_dewpoints = compute_adjusted_values(sea, elev, DEW_POINT_LAPSE_RATE) + + rh_values = compute_rh( + grid.temperature_data[idw_result.rows, idw_result.cols], + adjusted_dewpoints, + ) + rh_array[idw_result.rows, idw_result.cols] = rh_values + + log_interpolation_stats( + idw_result.total_pixels, + idw_result.interpolated_count, + idw_result.failed_interpolation_count, + idw_result.skipped_nodata_count, + ) + + if idw_result.interpolated_count == 0: + raise RuntimeError( + f"No pixels were successfully interpolated from {len(self.field.lats)} station(s). " + "Check that station coordinates fall within the raster extent and that at least " + "one station has a valid dew point value." ) - if interpolated_count == 0: - raise RuntimeError( - f"No pixels were successfully interpolated from {len(station_lats)} station(s). " - "Check that station coordinates fall within the raster extent and that at least " - "one station has a valid dew point value." - ) - - return WPSDataset.from_array( - array=rh_array, - geotransform=geo_transform, - projection=projection, - nodata_value=SFMS_NO_DATA, - datatype=gdal.GDT_Float32, - ) + return WPSDataset.from_array( + array=rh_array, + geotransform=grid.geotransform, + projection=grid.projection, + nodata_value=SFMS_NO_DATA, + datatype=gdal.GDT_Float32, + ) diff --git a/backend/packages/wps-sfms/src/wps_sfms/processors/temperature.py b/backend/packages/wps-sfms/src/wps_sfms/processors/temperature.py index 07986dc0d..de0936fa3 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/processors/temperature.py +++ b/backend/packages/wps-sfms/src/wps_sfms/processors/temperature.py @@ -1,101 +1,66 @@ import logging import numpy as np from osgeo import gdal -from wps_sfms.interpolation.source import LAPSE_RATE, StationTemperatureSource +from wps_sfms.interpolation.field import LAPSE_RATE, ScalarField, compute_adjusted_values from wps_shared.geospatial.wps_dataset import WPSDataset -from wps_shared.geospatial.spatial_interpolation import idw_interpolation from wps_sfms.interpolation.common import SFMS_NO_DATA, log_interpolation_stats -from wps_sfms.processors.idw import Interpolator +from wps_sfms.interpolation.grid import build_grid_context +from wps_sfms.processors.idw import RasterProcessor, idw_on_valid_pixels logger = logging.getLogger(__name__) -class TemperatureInterpolator(Interpolator): +class TemperatureInterpolator(RasterProcessor): """Interpolates station temperatures using IDW with elevation adjustment.""" - def __init__(self, mask_path: str, dem_path: str): + def __init__(self, mask_path: str, dem_path: str, field: ScalarField): super().__init__(mask_path) self.dem_path = dem_path - - def interpolate( - self, source: StationTemperatureSource, reference_raster_path: str - ) -> WPSDataset: - with WPSDataset(reference_raster_path) as ref_ds: - geo_transform = ref_ds.ds.GetGeoTransform() - if geo_transform is None: - raise ValueError( - f"Failed to get geotransform from reference raster: {reference_raster_path}" - ) - projection = ref_ds.ds.GetProjection() - x_size = ref_ds.ds.RasterXSize - y_size = ref_ds.ds.RasterYSize - - with WPSDataset(self.dem_path) as dem_ds: - dem_band: gdal.Band = dem_ds.ds.GetRasterBand(1) - dem_data = dem_band.ReadAsArray() - if dem_data is None: - raise ValueError("Failed to read DEM data") - - temp_array = np.full((y_size, x_size), SFMS_NO_DATA, dtype=np.float32) - - with WPSDataset(self.mask_path) as mask_ds: - valid_mask = ref_ds.apply_mask(mask_ds.warp_to_match(ref_ds)) - - lats, lons, valid_yi, valid_xi = dem_ds.get_lat_lon_coords(valid_mask) - valid_elevations = dem_data[valid_mask] - - total_pixels = x_size * y_size - skipped_nodata_count = total_pixels - len(valid_yi) - - logger.info("Interpolating temperature for raster grid (%d x %d)", x_size, y_size) - logger.info( - "Processing %d valid pixels (skipping %d NoData pixels)", - len(valid_yi), - skipped_nodata_count, - ) - - station_lats, station_lons, sea_level_temps = source.get_interpolation_data() - logger.info( - "Running batch temperature IDW interpolation for %d pixels and %d stations", - len(lats), - len(station_lats), - ) - - interpolated_sea_level_temps = idw_interpolation( - lats, lons, station_lats, station_lons, sea_level_temps - ) - assert isinstance(interpolated_sea_level_temps, np.ndarray) - - interpolation_succeeded = ~np.isnan(interpolated_sea_level_temps) - interpolated_count = int(np.sum(interpolation_succeeded)) - failed_interpolation_count = len(interpolated_sea_level_temps) - interpolated_count - - rows = valid_yi[interpolation_succeeded] - cols = valid_xi[interpolation_succeeded] - - sea = interpolated_sea_level_temps[interpolation_succeeded].astype( - np.float32, copy=False - ) - elev = valid_elevations[interpolation_succeeded].astype(np.float32, copy=False) - - actual_temps = source.compute_adjusted_values(sea, elev, LAPSE_RATE) - temp_array[rows, cols] = actual_temps - - log_interpolation_stats( - total_pixels, interpolated_count, failed_interpolation_count, skipped_nodata_count + self.field = field + + def interpolate(self, reference_raster_path: str) -> WPSDataset: + grid = build_grid_context(reference_raster_path, self.mask_path, dem_path=self.dem_path) + assert grid.valid_dem_values is not None + + temp_array = np.full((grid.y_size, grid.x_size), SFMS_NO_DATA, dtype=np.float32) + + logger.info("Interpolating temperature for raster grid (%d x %d)", grid.x_size, grid.y_size) + + idw_result = idw_on_valid_pixels( + valid_lats=grid.valid_lats, + valid_lons=grid.valid_lons, + valid_yi=grid.valid_yi, + valid_xi=grid.valid_xi, + station_lats=self.field.lats, + station_lons=self.field.lons, + station_values=self.field.values, + total_pixels=grid.total_pixels, + label="temperature", + ) + + sea = idw_result.values + elev = grid.valid_dem_values[idw_result.succeeded_mask].astype(np.float32, copy=False) + actual_temps = compute_adjusted_values(sea, elev, LAPSE_RATE) + temp_array[idw_result.rows, idw_result.cols] = actual_temps + + log_interpolation_stats( + idw_result.total_pixels, + idw_result.interpolated_count, + idw_result.failed_interpolation_count, + idw_result.skipped_nodata_count, + ) + + if idw_result.interpolated_count == 0: + raise RuntimeError( + f"No pixels were successfully interpolated from {len(self.field.lats)} station(s). " + "Check that station coordinates fall within the raster extent and that at least " + "one station has a valid temperature value." ) - if interpolated_count == 0: - raise RuntimeError( - f"No pixels were successfully interpolated from {len(station_lats)} station(s). " - "Check that station coordinates fall within the raster extent and that at least " - "one station has a valid temperature value." - ) - - return WPSDataset.from_array( - array=temp_array, - geotransform=geo_transform, - projection=projection, - nodata_value=SFMS_NO_DATA, - datatype=gdal.GDT_Float32, - ) + return WPSDataset.from_array( + array=temp_array, + geotransform=grid.geotransform, + projection=grid.projection, + nodata_value=SFMS_NO_DATA, + datatype=gdal.GDT_Float32, + ) diff --git a/backend/packages/wps-sfms/src/wps_sfms/processors/wind.py b/backend/packages/wps-sfms/src/wps_sfms/processors/wind.py new file mode 100644 index 000000000..fc00ef7dd --- /dev/null +++ b/backend/packages/wps-sfms/src/wps_sfms/processors/wind.py @@ -0,0 +1,118 @@ +import logging + +import numpy as np +from osgeo import gdal +from wps_shared.geospatial.wps_dataset import WPSDataset + +from wps_sfms.interpolation.common import SFMS_NO_DATA, log_interpolation_stats +from wps_sfms.interpolation.field import WindVectorField +from wps_sfms.interpolation.grid import build_grid_context +from wps_sfms.processors.idw import Interpolator, RasterProcessor, idw_on_valid_pixels + +logger = logging.getLogger(__name__) + + +class WindSpeedInterpolator(Interpolator): + """Interpolates wind speed using base IDW workflow.""" + + +class WindDirectionInterpolator(RasterProcessor): + """Interpolates wind direction by IDW on u/v components then reconstructing direction.""" + + def __init__(self, mask_path: str, field: WindVectorField): + super().__init__(mask_path) + self.field = field + + @staticmethod + def compute_direction_from_uv(u: np.ndarray, v: np.ndarray) -> np.ndarray: + """ + Match legacy SFMS wind direction reconstruction from interpolated u/v. + https://github.com/cffdrs/sfms/blob/main/src/SfmsWeather.cpp + + :param u: u component of wind vector + :param v: v component of wind vector + :return: wind direction in degrees + """ + u = np.asarray(u, dtype=np.float32) + v = np.asarray(v, dtype=np.float32) + direction = np.zeros(u.shape, dtype=np.float32) + + zero_v = np.abs(v) < np.float32(1e-6) + nonzero_v = ~zero_v + + direction[nonzero_v] = ( + np.degrees(np.arctan2(u[nonzero_v], v[nonzero_v])) + np.float32(180.0) + ).astype(np.float32) + + zero_u = np.abs(u) < np.float32(1e-6) + + # when v is effectively zero, avoid unstable/ambiguous atan2 results by applying legacy + # SFMS rules: u<0 => 90 deg, u>0 => 270 deg, and u~=0 & v~=0 (calm/no directional signal) => 0 deg. + direction[zero_v & (u < 0.0)] = 90.0 + direction[zero_v & (u > 0.0)] = 270.0 + direction[zero_v & zero_u] = 0.0 + return direction + + def interpolate(self, reference_raster_path: str) -> WPSDataset: + grid = build_grid_context(reference_raster_path, self.mask_path) + wind_direction_array = np.full((grid.y_size, grid.x_size), SFMS_NO_DATA, dtype=np.float32) + + logger.info("Interpolating wind direction for raster grid (%d x %d)", grid.x_size, grid.y_size) + + u_result = idw_on_valid_pixels( + valid_lats=grid.valid_lats, + valid_lons=grid.valid_lons, + valid_yi=grid.valid_yi, + valid_xi=grid.valid_xi, + station_lats=self.field.lats, + station_lons=self.field.lons, + station_values=self.field.u, + total_pixels=grid.total_pixels, + label="wind-u component", + ) + v_result = idw_on_valid_pixels( + valid_lats=grid.valid_lats, + valid_lons=grid.valid_lons, + valid_yi=grid.valid_yi, + valid_xi=grid.valid_xi, + station_lats=self.field.lats, + station_lons=self.field.lons, + station_values=self.field.v, + total_pixels=grid.total_pixels, + label="wind-v component", + ) + + wind_success = u_result.succeeded_mask & v_result.succeeded_mask + interpolated_count = int(np.sum(wind_success)) + failed_interpolation_count = len(wind_success) - interpolated_count + + if interpolated_count > 0: + rows = grid.valid_yi[wind_success] + cols = grid.valid_xi[wind_success] + directions = self.compute_direction_from_uv( + u_result.interpolated_values[wind_success], + v_result.interpolated_values[wind_success], + ) + wind_direction_array[rows, cols] = directions + + log_interpolation_stats( + total_pixels=grid.total_pixels, + interpolated_count=interpolated_count, + failed_interpolation_count=failed_interpolation_count, + skipped_nodata_count=u_result.skipped_nodata_count, + ) + + if interpolated_count == 0: + raise RuntimeError( + f"No pixels were successfully interpolated from {len(self.field.lats)} station(s). " + "Check that station coordinates fall within the raster extent and that at least " + "one station has both a valid wind speed and wind direction value." + ) + + return WPSDataset.from_array( + array=wind_direction_array, + geotransform=grid.geotransform, + projection=grid.projection, + nodata_value=SFMS_NO_DATA, + datatype=gdal.GDT_Float32, + ) diff --git a/backend/packages/wps-sfms/src/wps_sfms/publish.py b/backend/packages/wps-sfms/src/wps_sfms/publish.py new file mode 100644 index 000000000..508b4a9b4 --- /dev/null +++ b/backend/packages/wps-sfms/src/wps_sfms/publish.py @@ -0,0 +1,53 @@ +"""Shared helpers for publishing raster outputs and their derived COGs.""" + +import logging +import os +import tempfile +from dataclasses import dataclass + +import aiofiles + +from wps_sfms.sfmsng_raster_addresser import SFMSNGRasterAddresser +from wps_shared.geospatial.cog import generate_web_optimized_cog +from wps_shared.geospatial.wps_dataset import WPSDataset +from wps_shared.sfms.raster_addresser import GDALPath, S3Key +from wps_shared.utils.s3 import set_s3_gdal_config +from wps_shared.utils.s3_client import S3Client + +logger = logging.getLogger(__name__) + + +@dataclass(frozen=True) +class PublishedRaster: + """Storage locations for a published GeoTIFF and its derived COG.""" + + output_key: S3Key + cog_key: GDALPath | None + + +async def publish_dataset( + s3_client: S3Client, + dataset: WPSDataset, + output_key: S3Key | str, + generate_cog: bool = True, +) -> PublishedRaster: + """Upload a GeoTIFF to object storage and optionally generate a matching web COG.""" + + raster_addresser = SFMSNGRasterAddresser() + s3_output_key = S3Key(str(output_key)) + cog_key = raster_addresser.get_cog_key(s3_output_key) if generate_cog else None + + set_s3_gdal_config() + + with tempfile.TemporaryDirectory() as tmp_dir: + tmp_path = os.path.join(tmp_dir, os.path.basename(s3_output_key)) + dataset.export_to_geotiff(tmp_path) + + logger.info("Writing raster to S3: %s", s3_output_key) + async with aiofiles.open(tmp_path, "rb") as f: + await s3_client.put_object(key=s3_output_key, body=await f.read()) + + if cog_key is not None: + generate_web_optimized_cog(input_path=tmp_path, output_path=cog_key) + + return PublishedRaster(output_key=s3_output_key, cog_key=cog_key) diff --git a/backend/packages/wps-sfms/src/wps_sfms/sfmsng_raster_addresser.py b/backend/packages/wps-sfms/src/wps_sfms/sfmsng_raster_addresser.py index 302eb3c78..db48cb88d 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/sfmsng_raster_addresser.py +++ b/backend/packages/wps-sfms/src/wps_sfms/sfmsng_raster_addresser.py @@ -7,6 +7,7 @@ from dataclasses import dataclass from datetime import datetime, timedelta +from typing import Mapping from wps_shared.run_type import RunType from wps_shared.sfms.raster_addresser import ( @@ -21,19 +22,15 @@ @dataclass(frozen=True) class FWIInputs: - """All S3 keys and metadata needed for a single FWI index calculation. + """All input key mappings and metadata needed for a single FWI calculation. - Input keys (temp_key, rh_key, precip_key, prev_fwi_key, cog_key) are - GDALPath values (/vsis3/...) for reading via GDAL. output_key is a plain - S3Key for writing via boto3. + `weather_keys` and `index_keys` values are GDALPath values (/vsis3/...) for + reading via GDAL. `output_key` is a plain S3Key for writing via boto3. """ - temp_key: GDALPath - rh_key: GDALPath - precip_key: GDALPath - prev_fwi_key: GDALPath + weather_keys: Mapping[SFMSInterpolatedWeatherParameter, GDALPath] + index_keys: Mapping[FWIParameter, GDALPath] output_key: S3Key - cog_key: GDALPath run_type: RunType @@ -61,7 +58,9 @@ def get_actual_weather_key( date = datetime_utc.date() param = weather_param.value date_str = date.isoformat().replace("-", "") - return S3Key(f"{self.root}/actual/{date.year:04d}/{date.month:02d}/{date.day:02d}/{param}_{date_str}.tif") + return S3Key( + f"{self.root}/actual/{date.year:04d}/{date.month:02d}/{date.day:02d}/{param}_{date_str}.tif" + ) def get_actual_index_key(self, datetime_utc: datetime, fwi_param: FWIParameter) -> S3Key: """ @@ -72,38 +71,74 @@ def get_actual_index_key(self, datetime_utc: datetime, fwi_param: FWIParameter) assert_all_utc(datetime_utc) date = datetime_utc.date() date_str = date.isoformat().replace("-", "") - return S3Key(f"{self.root}/actual/{date.year:04d}/{date.month:02d}/{date.day:02d}/{fwi_param.value}_{date_str}.tif") + return S3Key( + f"{self.root}/actual/{date.year:04d}/{date.month:02d}/{date.day:02d}/{fwi_param.value}_{date_str}.tif" + ) def get_actual_fwi_inputs( self, datetime_to_process: datetime, fwi_param: FWIParameter ) -> FWIInputs: """ - Build a FWIInputs for a station-interpolated actual run. + Build FWIInputs for one actual-run index calculation. - Uses yesterday's uploaded FWI value as seed and today's interpolated - weather rasters (temp, rh, precip) as inputs. + Dependency keys vary by requested index: + - FFMC, DMC, DC: yesterday's same index + today's weather + - ISI: today's FFMC + today's wind speed + - BUI: today's DMC + today's DC + - FWI: today's ISI + today's BUI :param datetime_to_process: UTC datetime being processed - :param fwi_param: Which FWI parameter to calculate (FFMC, DMC, or DC) + :param fwi_param: Which FWI parameter to calculate :return: FWIInputs ready for FWIProcessor """ assert_all_utc(datetime_to_process) yesterday = datetime_to_process - timedelta(days=1) - temp_key, rh_key, precip_key, prev_fwi_key = self.gdal_prefix_keys( - self.get_actual_weather_key(datetime_to_process, SFMSInterpolatedWeatherParameter.TEMP), - self.get_actual_weather_key(datetime_to_process, SFMSInterpolatedWeatherParameter.RH), - self.get_actual_weather_key( - datetime_to_process, SFMSInterpolatedWeatherParameter.PRECIP - ), - self.get_actual_index_key(yesterday, fwi_param), - ) + + weather_keys = { + param: self.gdal_path(self.get_actual_weather_key(datetime_to_process, param)) + for param in ( + SFMSInterpolatedWeatherParameter.TEMP, + SFMSInterpolatedWeatherParameter.RH, + SFMSInterpolatedWeatherParameter.PRECIP, + SFMSInterpolatedWeatherParameter.WIND_SPEED, + ) + } + + if fwi_param in (FWIParameter.FFMC, FWIParameter.DMC, FWIParameter.DC): + index_keys = { + fwi_param: self.gdal_path(self.get_actual_index_key(yesterday, fwi_param)) + } + elif fwi_param == FWIParameter.ISI: + index_keys = { + FWIParameter.FFMC: self.gdal_path( + self.get_actual_index_key(datetime_to_process, FWIParameter.FFMC) + ) + } + elif fwi_param == FWIParameter.BUI: + index_keys = { + FWIParameter.DMC: self.gdal_path( + self.get_actual_index_key(datetime_to_process, FWIParameter.DMC) + ), + FWIParameter.DC: self.gdal_path( + self.get_actual_index_key(datetime_to_process, FWIParameter.DC) + ), + } + elif fwi_param == FWIParameter.FWI: + index_keys = { + FWIParameter.ISI: self.gdal_path( + self.get_actual_index_key(datetime_to_process, FWIParameter.ISI) + ), + FWIParameter.BUI: self.gdal_path( + self.get_actual_index_key(datetime_to_process, FWIParameter.BUI) + ), + } + else: + raise ValueError(f"Unsupported FWI parameter: {fwi_param.value}") + output_key = self.get_actual_index_key(datetime_to_process, fwi_param) return FWIInputs( - temp_key=temp_key, - rh_key=rh_key, - precip_key=precip_key, - prev_fwi_key=prev_fwi_key, + weather_keys=weather_keys, + index_keys=index_keys, output_key=output_key, - cog_key=self.get_cog_key(output_key), run_type=RunType.ACTUAL, ) diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_field.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_field.py new file mode 100644 index 000000000..ffb6e5a2a --- /dev/null +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_field.py @@ -0,0 +1,66 @@ +import numpy as np +import pytest +from wps_shared.schemas.sfms import SFMSDailyActual + +from wps_sfms.interpolation.field import ( + build_attribute_field, + build_dewpoint_field, + build_temperature_field, + build_wind_vector_field, +) + + +class TestScalarFieldBuilders: + def test_build_attribute_field_filters_missing_values(self): + actuals = [ + SFMSDailyActual(code=1, lat=49.0, lon=-123.0, precipitation=1.5), + SFMSDailyActual(code=2, lat=49.1, lon=-123.1, precipitation=None), + ] + + field = build_attribute_field(actuals, "precipitation") + + np.testing.assert_allclose(field.lats, np.array([49.0], dtype=np.float32)) + np.testing.assert_allclose(field.lons, np.array([-123.0], dtype=np.float32)) + np.testing.assert_allclose(field.values, np.array([1.5], dtype=np.float32)) + + def test_build_temperature_field_applies_sea_level_adjustment(self): + actuals = [ + SFMSDailyActual( + code=1, lat=49.0, lon=-123.0, elevation=100.0, temperature=15.0 + ) + ] + + field = build_temperature_field(actuals) + + np.testing.assert_allclose(field.values, np.array([15.65], dtype=np.float32), atol=1e-4) + + def test_build_dewpoint_field_skips_missing_elevation_or_value(self): + actuals = [ + SFMSDailyActual(code=1, lat=49.0, lon=-123.0, elevation=None, dewpoint=10.0), + SFMSDailyActual(code=2, lat=49.1, lon=-123.1, elevation=100.0, dewpoint=None), + ] + + field = build_dewpoint_field(actuals) + + assert field.lats.size == 0 + assert field.lons.size == 0 + assert field.values.size == 0 + + +class TestWindVectorFieldBuilder: + def test_build_wind_vector_field_filters_unpaired_values(self): + actuals = [ + SFMSDailyActual(code=1, lat=49.0, lon=-123.0, wind_speed=10.0, wind_direction=90.0), + SFMSDailyActual(code=2, lat=49.1, lon=-123.1, wind_speed=8.0, wind_direction=None), + ] + + field = build_wind_vector_field(actuals) + + np.testing.assert_allclose(field.lats, np.array([49.0], dtype=np.float32)) + np.testing.assert_allclose(field.lons, np.array([-123.0], dtype=np.float32)) + np.testing.assert_allclose(field.u, np.array([-10.0], dtype=np.float32), atol=1e-5) + np.testing.assert_allclose(field.v, np.array([0.0], dtype=np.float32), atol=1e-5) + + def test_build_attribute_field_rejects_unknown_attribute(self): + with pytest.raises(ValueError, match="Unknown attribute"): + build_attribute_field([], "not_a_real_field") diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_fwi_processor.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_fwi_processor.py index 09057ffbb..0b98c84cf 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/tests/test_fwi_processor.py +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_fwi_processor.py @@ -7,20 +7,23 @@ from wps_shared.geospatial.wps_dataset import WPSDataset from wps_shared.run_type import RunType -from wps_shared.sfms.raster_addresser import FWIParameter +from wps_shared.sfms.raster_addresser import FWIParameter, SFMSInterpolatedWeatherParameter from wps_sfms.sfmsng_raster_addresser import FWIInputs from wps_shared.tests.geospatial.dataset_common import ( - create_mock_gdal_dataset, create_mock_input_dataset_context, create_test_dataset, ) from wps_shared.utils.s3_client import S3Client from wps_sfms.processors.fwi import ( + BUICalculator, DCCalculator, DMCCalculator, FFMCCalculator, FWICalculator, + FWIDatasets, + FWIFinalCalculator, FWIProcessor, + ISICalculator, ) TEST_DATETIME = datetime(2024, 10, 10, 20, tzinfo=timezone.utc) @@ -45,13 +48,39 @@ def make_fwi_inputs(fwi_param: FWIParameter, run_type: RunType = RunType.ACTUAL) prev_date_iso = "2024-10-09" s3_prefix = "/vsis3/test-bucket" output_key = f"sfms/calculated/{run_type.value}/{date_iso}/{param}{date_str}.tif" + + weather_keys = { + SFMSInterpolatedWeatherParameter.TEMP: f"{s3_prefix}/sfms/interpolated/temp/2024/10/10/temp_{date_str}.tif", + SFMSInterpolatedWeatherParameter.RH: f"{s3_prefix}/sfms/interpolated/rh/2024/10/10/rh_{date_str}.tif", + SFMSInterpolatedWeatherParameter.PRECIP: f"{s3_prefix}/sfms/interpolated/precip/2024/10/10/precip_{date_str}.tif", + SFMSInterpolatedWeatherParameter.WIND_SPEED: f"{s3_prefix}/sfms/interpolated/wind_speed/2024/10/10/wind_speed_{date_str}.tif", + } + + if fwi_param in (FWIParameter.FFMC, FWIParameter.DMC, FWIParameter.DC): + index_keys = { + fwi_param: f"{s3_prefix}/sfms/uploads/actual/{prev_date_iso}/{param}{prev_date_str}.tif" + } + elif fwi_param == FWIParameter.ISI: + index_keys = { + FWIParameter.FFMC: f"{s3_prefix}/sfms/calculated/actual/{date_iso}/ffmc{date_str}.tif" + } + elif fwi_param == FWIParameter.BUI: + index_keys = { + FWIParameter.DMC: f"{s3_prefix}/sfms/calculated/actual/{date_iso}/dmc{date_str}.tif", + FWIParameter.DC: f"{s3_prefix}/sfms/calculated/actual/{date_iso}/dc{date_str}.tif", + } + elif fwi_param == FWIParameter.FWI: + index_keys = { + FWIParameter.ISI: f"{s3_prefix}/sfms/calculated/actual/{date_iso}/isi{date_str}.tif", + FWIParameter.BUI: f"{s3_prefix}/sfms/calculated/actual/{date_iso}/bui{date_str}.tif", + } + else: + raise AssertionError(f"Unhandled FWI parameter for test: {fwi_param}") + return FWIInputs( - temp_key=f"{s3_prefix}/sfms/interpolated/temp/2024/10/10/temp_{date_str}.tif", - rh_key=f"{s3_prefix}/sfms/interpolated/rh/2024/10/10/rh_{date_str}.tif", - precip_key=f"{s3_prefix}/sfms/interpolated/precip/2024/10/10/precip_{date_str}.tif", - prev_fwi_key=f"{s3_prefix}/sfms/uploads/actual/{prev_date_iso}/{param}{prev_date_str}.tif", + weather_keys=weather_keys, + index_keys=index_keys, output_key=output_key, - cog_key=f"{s3_prefix}/{output_key.removesuffix('.tif')}_cog.tif", run_type=run_type, ) @@ -62,18 +91,17 @@ async def test_fwi_processor_ffmc(mocker: MockerFixture): processor = FWIProcessor(TEST_DATETIME) fwi_inputs = make_fwi_inputs(FWIParameter.FFMC) - _, mock_input_dataset_context = create_mock_input_dataset_context(4) + _, mock_input_dataset_context = create_mock_input_dataset_context(5) - mocker.patch("osgeo.gdal.Open", return_value=create_mock_gdal_dataset()) - generate_and_store_cog_spy = mocker.patch("wps_sfms.processors.fwi.generate_and_store_cog") + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) rasters_match_spy = mocker.patch("wps_sfms.processors.fwi.rasters_match", return_value=True) async with S3Client() as 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 processor.calculate_index( mock_s3_client, mock_input_dataset_context, FFMCCalculator(), fwi_inputs @@ -82,18 +110,11 @@ async def test_fwi_processor_ffmc(mocker: MockerFixture): # Verify weather + FWI keys were checked assert mock_all_objects_exist.call_count == 2 - # Verify all three weather rasters were checked against the FWI grid - assert rasters_match_spy.call_count == 3 + # Verify all required weather rasters were checked against the FWI grid + assert rasters_match_spy.call_count == 4 - # Verify output was persisted with the correct key - assert persist_raster_spy.call_count == 1 - assert persist_raster_spy.call_args[0][1] == fwi_inputs.output_key - - # Verify COG was generated with the correct key - assert generate_and_store_cog_spy.call_count == 1 - generate_and_store_cog_spy.assert_called_once_with( - src_ds=mocker.ANY, output_path=fwi_inputs.cog_key - ) + assert publish_spy.call_count == 1 + assert publish_spy.await_args.kwargs["output_key"] == fwi_inputs.output_key @pytest.mark.anyio @@ -104,14 +125,13 @@ async def test_fwi_processor_dmc(mocker: MockerFixture): _, mock_input_dataset_context = create_mock_input_dataset_context(4) - mocker.patch("osgeo.gdal.Open", return_value=create_mock_gdal_dataset()) - mocker.patch("wps_sfms.processors.fwi.generate_and_store_cog") + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) async with S3Client() as mock_s3_client: mocker.patch.object(mock_s3_client, "all_objects_exist", new=AsyncMock(return_value=True)) - persist_raster_spy = mocker.patch.object( - mock_s3_client, "persist_raster_data", return_value="test_key.tif" - ) await processor.calculate_index( mock_s3_client, @@ -120,8 +140,8 @@ async def test_fwi_processor_dmc(mocker: MockerFixture): fwi_inputs, ) - assert persist_raster_spy.call_count == 1 - assert persist_raster_spy.call_args[0][1] == fwi_inputs.output_key + assert publish_spy.call_count == 1 + assert publish_spy.await_args.kwargs["output_key"] == fwi_inputs.output_key @pytest.mark.anyio @@ -132,14 +152,13 @@ async def test_fwi_processor_dc(mocker: MockerFixture): _, mock_input_dataset_context = create_mock_input_dataset_context(4) - mocker.patch("osgeo.gdal.Open", return_value=create_mock_gdal_dataset()) - mocker.patch("wps_sfms.processors.fwi.generate_and_store_cog") + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) async with S3Client() as mock_s3_client: mocker.patch.object(mock_s3_client, "all_objects_exist", new=AsyncMock(return_value=True)) - persist_raster_spy = mocker.patch.object( - mock_s3_client, "persist_raster_data", return_value="test_key.tif" - ) await processor.calculate_index( mock_s3_client, @@ -148,8 +167,8 @@ async def test_fwi_processor_dc(mocker: MockerFixture): fwi_inputs, ) - assert persist_raster_spy.call_count == 1 - assert persist_raster_spy.call_args[0][1] == fwi_inputs.output_key + assert publish_spy.call_count == 1 + assert publish_spy.await_args.kwargs["output_key"] == fwi_inputs.output_key @pytest.mark.anyio @@ -162,19 +181,22 @@ async def test_fwi_processor_missing_weather_keys(mocker: MockerFixture): async with S3Client() as mock_s3_client: mocker.patch.object(mock_s3_client, "all_objects_exist", new=AsyncMock(return_value=False)) - persist_raster_spy = mocker.patch.object(mock_s3_client, "persist_raster_data") + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) - with pytest.raises(RuntimeError, match="Missing weather keys"): + with pytest.raises(RuntimeError, match="Missing weather dependency keys"): await processor.calculate_index( mock_s3_client, mock_input_dataset_context, FFMCCalculator(), fwi_inputs ) - persist_raster_spy.assert_not_called() + publish_spy.assert_not_called() @pytest.mark.anyio async def test_fwi_processor_missing_fwi_keys(mocker: MockerFixture): - """Test that processor bails when the previous day's FWI key is missing.""" + """Test that processor bails when required index dependency keys are missing.""" processor = FWIProcessor(TEST_DATETIME) fwi_inputs = make_fwi_inputs(FWIParameter.DMC) @@ -184,9 +206,12 @@ async def test_fwi_processor_missing_fwi_keys(mocker: MockerFixture): mocker.patch.object( mock_s3_client, "all_objects_exist", new=AsyncMock(side_effect=[True, False]) ) - persist_raster_spy = mocker.patch.object(mock_s3_client, "persist_raster_data") + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) - with pytest.raises(RuntimeError, match="Missing previous dmc raster for"): + with pytest.raises(RuntimeError, match="Missing index dependency keys"): await processor.calculate_index( mock_s3_client, mock_input_dataset_context, @@ -194,18 +219,19 @@ async def test_fwi_processor_missing_fwi_keys(mocker: MockerFixture): fwi_inputs, ) - persist_raster_spy.assert_not_called() + publish_spy.assert_not_called() @pytest.mark.anyio @pytest.mark.parametrize( "match_side_effects,expected_message", [ - ([False], "Temperature raster does not match FWI grid"), - ([True, False], "RH raster does not match FWI grid"), - ([True, True, False], "Precip raster does not match FWI grid"), + ([False], "temperature raster does not match FWI grid"), + ([True, False], "relative_humidity raster does not match FWI grid"), + ([True, True, False], "precipitation raster does not match FWI grid"), + ([True, True, True, False], "wind_speed raster does not match FWI grid"), ], - ids=["temp_mismatch", "rh_mismatch", "precip_mismatch"], + ids=["temp_mismatch", "rh_mismatch", "precip_mismatch", "wind_mismatch"], ) async def test_fwi_processor_raster_mismatch_raises( mocker: MockerFixture, match_side_effects, expected_message @@ -214,19 +240,22 @@ async def test_fwi_processor_raster_mismatch_raises( processor = FWIProcessor(TEST_DATETIME) fwi_inputs = make_fwi_inputs(FWIParameter.FFMC) - _, mock_input_dataset_context = create_mock_input_dataset_context(4) + _, mock_input_dataset_context = create_mock_input_dataset_context(5) mocker.patch("wps_sfms.processors.fwi.rasters_match", side_effect=match_side_effects) async with S3Client() as mock_s3_client: mocker.patch.object(mock_s3_client, "all_objects_exist", new=AsyncMock(return_value=True)) - persist_raster_spy = mocker.patch.object(mock_s3_client, "persist_raster_data") + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) with pytest.raises(ValueError, match=expected_message): await processor.calculate_index( mock_s3_client, mock_input_dataset_context, FFMCCalculator(), fwi_inputs ) - persist_raster_spy.assert_not_called() + publish_spy.assert_not_called() @pytest.mark.anyio @@ -242,25 +271,111 @@ async def test_fwi_processor_run_type_in_output_key(mocker: MockerFixture): assert "forecast" in forecast_inputs.output_key assert "actual" not in forecast_inputs.output_key - _, mock_input_dataset_context = create_mock_input_dataset_context(4) - mocker.patch("osgeo.gdal.Open", return_value=create_mock_gdal_dataset()) - mocker.patch("wps_sfms.processors.fwi.generate_and_store_cog") + _, mock_input_dataset_context = create_mock_input_dataset_context(5) + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) async with S3Client() as mock_s3_client: mocker.patch.object(mock_s3_client, "all_objects_exist", new=AsyncMock(return_value=True)) - persist_raster_spy = mocker.patch.object( - mock_s3_client, "persist_raster_data", return_value="test_key.tif" - ) await processor.calculate_index( mock_s3_client, mock_input_dataset_context, FFMCCalculator(), actual_inputs ) - output_key = persist_raster_spy.call_args[0][1] + output_key = publish_spy.await_args.kwargs["output_key"] assert "actual" in output_key assert "forecast" not in output_key +@pytest.mark.anyio +async def test_fwi_processor_isi(mocker: MockerFixture): + """Test that calculate_index with ISICalculator produces output with correct key.""" + processor = FWIProcessor(TEST_DATETIME) + fwi_inputs = make_fwi_inputs(FWIParameter.ISI) + + _, mock_input_dataset_context = create_mock_input_dataset_context(2) + + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) + rasters_match_spy = mocker.patch("wps_sfms.processors.fwi.rasters_match", return_value=True) + + async with S3Client() as mock_s3_client: + mocker.patch.object(mock_s3_client, "all_objects_exist", new=AsyncMock(return_value=True)) + + await processor.calculate_index( + mock_s3_client, + mock_input_dataset_context, + ISICalculator(), + fwi_inputs, + ) + + assert rasters_match_spy.call_count == 1 + assert publish_spy.call_count == 1 + assert publish_spy.await_args.kwargs["output_key"] == fwi_inputs.output_key + + +@pytest.mark.anyio +async def test_fwi_processor_bui(mocker: MockerFixture): + """Test that calculate_index with BUICalculator produces output with correct key.""" + processor = FWIProcessor(TEST_DATETIME) + fwi_inputs = make_fwi_inputs(FWIParameter.BUI) + + _, mock_input_dataset_context = create_mock_input_dataset_context(2) + + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) + rasters_match_spy = mocker.patch("wps_sfms.processors.fwi.rasters_match", return_value=True) + + async with S3Client() as mock_s3_client: + mocker.patch.object(mock_s3_client, "all_objects_exist", new=AsyncMock(return_value=True)) + + await processor.calculate_index( + mock_s3_client, + mock_input_dataset_context, + BUICalculator(), + fwi_inputs, + ) + + assert rasters_match_spy.call_count == 1 + assert publish_spy.call_count == 1 + assert publish_spy.await_args.kwargs["output_key"] == fwi_inputs.output_key + + +@pytest.mark.anyio +async def test_fwi_processor_fwi(mocker: MockerFixture): + """Test that calculate_index with FWIFinalCalculator produces output with correct key.""" + processor = FWIProcessor(TEST_DATETIME) + fwi_inputs = make_fwi_inputs(FWIParameter.FWI) + + _, mock_input_dataset_context = create_mock_input_dataset_context(2) + + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", + new=AsyncMock(), + ) + rasters_match_spy = mocker.patch("wps_sfms.processors.fwi.rasters_match", return_value=True) + + async with S3Client() as mock_s3_client: + mocker.patch.object(mock_s3_client, "all_objects_exist", new=AsyncMock(return_value=True)) + + await processor.calculate_index( + mock_s3_client, + mock_input_dataset_context, + FWIFinalCalculator(), + fwi_inputs, + ) + + assert rasters_match_spy.call_count == 1 + assert publish_spy.call_count == 1 + assert publish_spy.await_args.kwargs["output_key"] == fwi_inputs.output_key + + @pytest.mark.anyio async def test_fwi_processor_cog_failure_propagates(mocker: MockerFixture): """If COG generation fails after a successful persist, the error propagates. @@ -271,27 +386,22 @@ async def test_fwi_processor_cog_failure_propagates(mocker: MockerFixture): processor = FWIProcessor(TEST_DATETIME) fwi_inputs = make_fwi_inputs(FWIParameter.FFMC) - _, mock_input_dataset_context = create_mock_input_dataset_context(4) - mocker.patch("osgeo.gdal.Open", return_value=create_mock_gdal_dataset()) + _, mock_input_dataset_context = create_mock_input_dataset_context(5) mocker.patch("wps_sfms.processors.fwi.rasters_match", return_value=True) - mocker.patch( - "wps_sfms.processors.fwi.generate_and_store_cog", + publish_spy = mocker.patch( + "wps_sfms.processors.fwi.publish_dataset", side_effect=RuntimeError("COG generation failed"), ) async with S3Client() as mock_s3_client: mocker.patch.object(mock_s3_client, "all_objects_exist", new=AsyncMock(return_value=True)) - persist_raster_spy = mocker.patch.object( - mock_s3_client, "persist_raster_data", return_value="test_key.tif" - ) with pytest.raises(RuntimeError, match="COG generation failed"): await processor.calculate_index( mock_s3_client, mock_input_dataset_context, FFMCCalculator(), fwi_inputs ) - # Raster was persisted before COG generation was attempted - persist_raster_spy.assert_called_once() + publish_spy.assert_called_once() class TestFWINodeataPropagation: @@ -330,7 +440,18 @@ def test_nodata_propagates(self, calculator: FWICalculator, prev_value, nodata_i rh_ds = self.make_ds(50.0, nodata_at=nodata_pixel if nodata_input == "rh" else None) precip_ds = self.make_ds(0.0, nodata_at=nodata_pixel if nodata_input == "precip" else None) - result = calculator.calculate(prev_ds, temp_ds, rh_ds, precip_ds) + wind_ds = self.make_ds(10.0) + weather_datasets = { + SFMSInterpolatedWeatherParameter.TEMP: temp_ds, + SFMSInterpolatedWeatherParameter.RH: rh_ds, + SFMSInterpolatedWeatherParameter.PRECIP: precip_ds, + SFMSInterpolatedWeatherParameter.WIND_SPEED: wind_ds, + } + index_datasets = {calculator.reference_index_param: prev_ds} + + result = calculator.calculate( + FWIDatasets(index=index_datasets, weather=weather_datasets) + ) assert np.isnan(result.nodata_value) assert np.isnan(result.values[0, 0]), "nodata pixel must propagate as NaN" diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_grid_context.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_grid_context.py new file mode 100644 index 000000000..8692b78b7 --- /dev/null +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_grid_context.py @@ -0,0 +1,70 @@ +import uuid + +import numpy as np +import pytest +from osgeo import gdal + +from wps_sfms.interpolation.grid import build_grid_context +from wps_sfms.tests.conftest import create_test_raster + + +class TestBuildGridContext: + def test_loads_masked_grid_with_dem_and_temperature(self): + test_id = uuid.uuid4().hex + ref_path = f"/vsimem/reference_{test_id}.tif" + mask_path = f"/vsimem/mask_{test_id}.tif" + dem_path = f"/vsimem/dem_{test_id}.tif" + temp_path = f"/vsimem/temp_{test_id}.tif" + + try: + extent = (-123.1, -123.0, 49.0, 49.1) + create_test_raster(ref_path, 5, 5, extent, fill_value=1.0) + create_test_raster(dem_path, 5, 5, extent, fill_value=100.0) + create_test_raster(temp_path, 5, 5, extent, fill_value=15.0) + + mask_data = np.full((5, 5), 1.0, dtype=np.float32) + mask_data[2, 2] = 0.0 + create_test_raster(mask_path, 5, 5, extent, data=mask_data) + + grid = build_grid_context( + ref_path, + mask_path, + dem_path=dem_path, + temperature_raster_path=temp_path, + ) + + assert grid.x_size == 5 + assert grid.y_size == 5 + assert grid.total_pixels == 25 + assert grid.skipped_nodata_count == 1 + assert len(grid.valid_yi) == 24 + assert grid.valid_dem_values is not None + assert grid.temperature_data is not None + assert grid.valid_dem_values.shape == (24,) + assert grid.temperature_data.shape == (5, 5) + assert grid.valid_lats.dtype == np.float32 + assert grid.valid_lons.dtype == np.float32 + finally: + gdal.Unlink(ref_path) + gdal.Unlink(mask_path) + gdal.Unlink(dem_path) + gdal.Unlink(temp_path) + + def test_raises_when_dem_grid_does_not_match_reference(self): + test_id = uuid.uuid4().hex + ref_path = f"/vsimem/reference_{test_id}.tif" + mask_path = f"/vsimem/mask_{test_id}.tif" + dem_path = f"/vsimem/dem_{test_id}.tif" + + try: + extent = (-123.1, -123.0, 49.0, 49.1) + create_test_raster(ref_path, 5, 5, extent, fill_value=1.0) + create_test_raster(mask_path, 5, 5, extent, fill_value=1.0) + create_test_raster(dem_path, 4, 5, extent, fill_value=100.0) + + with pytest.raises(ValueError, match="DEM grid does not match reference raster"): + build_grid_context(ref_path, mask_path, dem_path=dem_path) + finally: + gdal.Unlink(ref_path) + gdal.Unlink(mask_path) + gdal.Unlink(dem_path) diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_idw_processor.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_idw_processor.py new file mode 100644 index 000000000..f84ef328a --- /dev/null +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_idw_processor.py @@ -0,0 +1,118 @@ +import numpy as np + +from wps_sfms.processors.idw import idw_on_valid_pixels + + +class TestIdwOnValidPixels: + def test_all_pixels_interpolate_successfully(self, monkeypatch): + # Verifies the full success path: all returned values are valid and + # are indexed back to the original valid pixel row/col arrays. + valid_lats = np.array([49.0, 49.1, 49.2], dtype=np.float32) + valid_lons = np.array([-123.0, -123.1, -123.2], dtype=np.float32) + valid_yi = np.array([1, 3, 4], dtype=np.intp) + valid_xi = np.array([2, 0, 5], dtype=np.intp) + station_lats = np.array([49.05, 49.15], dtype=np.float32) + station_lons = np.array([-123.05, -123.15], dtype=np.float32) + station_values = np.array([10.0, 20.0], dtype=np.float32) + + def fake_idw_interpolation(*args): + assert len(args) == 5 + return np.array([1.5, 2.5, 3.5], dtype=np.float64) + + monkeypatch.setattr("wps_sfms.processors.idw.idw_interpolation", fake_idw_interpolation) + + result = idw_on_valid_pixels( + valid_lats=valid_lats, + valid_lons=valid_lons, + valid_yi=valid_yi, + valid_xi=valid_xi, + station_lats=station_lats, + station_lons=station_lons, + station_values=station_values, + total_pixels=20, + label="test", + ) + + assert result.interpolated_values.dtype == np.float32 + np.testing.assert_allclose(result.interpolated_values, np.array([1.5, 2.5, 3.5])) + np.testing.assert_array_equal(result.succeeded_mask, np.array([True, True, True])) + np.testing.assert_array_equal(result.rows, valid_yi) + np.testing.assert_array_equal(result.cols, valid_xi) + np.testing.assert_allclose(result.values, np.array([1.5, 2.5, 3.5])) + assert result.interpolated_count == 3 + assert result.failed_interpolation_count == 0 + assert result.skipped_nodata_count == 17 + assert result.total_pixels == 20 + + def test_nan_interpolated_values_are_marked_failed(self, monkeypatch): + # Verifies NaN handling: failed pixels are excluded from rows/cols/values, + # while successful values retain their source indices. + valid_lats = np.array([49.0, 49.1, 49.2, 49.3], dtype=np.float32) + valid_lons = np.array([-123.0, -123.1, -123.2, -123.3], dtype=np.float32) + valid_yi = np.array([0, 1, 2, 3], dtype=np.intp) + valid_xi = np.array([9, 8, 7, 6], dtype=np.intp) + station_lats = np.array([49.05, 49.15], dtype=np.float32) + station_lons = np.array([-123.05, -123.15], dtype=np.float32) + station_values = np.array([10.0, 20.0], dtype=np.float32) + + monkeypatch.setattr( + "wps_sfms.processors.idw.idw_interpolation", + lambda *args: np.array([np.nan, 4.0, np.nan, -1.0], dtype=np.float64), + ) + + result = idw_on_valid_pixels( + valid_lats=valid_lats, + valid_lons=valid_lons, + valid_yi=valid_yi, + valid_xi=valid_xi, + station_lats=station_lats, + station_lons=station_lons, + station_values=station_values, + total_pixels=10, + label="test", + ) + + np.testing.assert_array_equal(result.succeeded_mask, np.array([False, True, False, True])) + np.testing.assert_array_equal(result.rows, np.array([1, 3], dtype=np.intp)) + np.testing.assert_array_equal(result.cols, np.array([8, 6], dtype=np.intp)) + np.testing.assert_allclose(result.values, np.array([4.0, -1.0], dtype=np.float32)) + assert result.interpolated_count == 2 + assert result.failed_interpolation_count == 2 + assert result.skipped_nodata_count == 6 + + def test_empty_valid_pixels_returns_empty_indexed_results(self, monkeypatch): + # Verifies edge case behavior when no valid raster pixels are provided: + # all outputs should be empty and skipped_nodata_count == total_pixels. + valid_lats = np.array([], dtype=np.float32) + valid_lons = np.array([], dtype=np.float32) + valid_yi = np.array([], dtype=np.intp) + valid_xi = np.array([], dtype=np.intp) + station_lats = np.array([49.05], dtype=np.float32) + station_lons = np.array([-123.05], dtype=np.float32) + station_values = np.array([10.0], dtype=np.float32) + + monkeypatch.setattr( + "wps_sfms.processors.idw.idw_interpolation", + lambda *args: np.array([], dtype=np.float64), + ) + + result = idw_on_valid_pixels( + valid_lats=valid_lats, + valid_lons=valid_lons, + valid_yi=valid_yi, + valid_xi=valid_xi, + station_lats=station_lats, + station_lons=station_lons, + station_values=station_values, + total_pixels=5, + label="test", + ) + + assert result.interpolated_values.size == 0 + assert result.succeeded_mask.size == 0 + assert result.rows.size == 0 + assert result.cols.size == 0 + assert result.values.size == 0 + assert result.interpolated_count == 0 + assert result.failed_interpolation_count == 0 + assert result.skipped_nodata_count == 5 diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_lapse_rate_adjusted_source.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_lapse_rate_adjusted_source.py index 8de616e87..a45fe1637 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/tests/test_lapse_rate_adjusted_source.py +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_lapse_rate_adjusted_source.py @@ -8,7 +8,7 @@ from hypothesis import given, strategies as st, settings import hypothesis.extra.numpy as hnp -from wps_sfms.interpolation.source import LAPSE_RATE, LapseRateAdjustedSource +from wps_sfms.interpolation.field import LAPSE_RATE, compute_adjusted_values, compute_sea_level_values finite_value_c = st.floats( @@ -52,8 +52,8 @@ def test_round_trip_property(values, elevs, lapse): values = values[:n].astype(np.float32, copy=False) elevs = elevs[:n].astype(np.float32, copy=False) - sea = LapseRateAdjustedSource.compute_sea_level_values(values, elevs, lapse) - back = LapseRateAdjustedSource.compute_adjusted_values(sea, elevs, lapse) + sea = compute_sea_level_values(values, elevs, lapse) + back = compute_adjusted_values(sea, elevs, lapse) assert_allclose(back, values, atol=1e-4) @@ -76,7 +76,7 @@ def test_monotone_cooling_with_positive_elevation(sea, elev, lapse): sea = sea[:n] elev = elev[:n] - out = LapseRateAdjustedSource.compute_adjusted_values(sea, elev, lapse) + out = compute_adjusted_values(sea, elev, lapse) assert np.all(out <= sea + 1e-6) @@ -99,5 +99,5 @@ def test_negative_elevation_warms(sea, elev, lapse): sea = sea[:n] elev = elev[:n] - out = LapseRateAdjustedSource.compute_adjusted_values(sea, elev, lapse) + out = compute_adjusted_values(sea, elev, lapse) assert np.all(out >= sea - 1e-6) diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_publish.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_publish.py new file mode 100644 index 000000000..58b2f19c8 --- /dev/null +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_publish.py @@ -0,0 +1,78 @@ +from unittest.mock import AsyncMock, MagicMock + +import pytest + +from wps_sfms.publish import publish_dataset +from wps_sfms.sfmsng_raster_addresser import SFMSNGRasterAddresser +from wps_shared.geospatial.wps_dataset import WPSDataset +from wps_shared.tests.geospatial.dataset_common import create_test_dataset + + +@pytest.mark.anyio +async def test_publish_dataset_uploads_raster_and_generates_cog_by_default(mocker): + output_key = "sfms_ng/actual/2024/07/04/temperature_20240704.tif" + expected_cog_key = SFMSNGRasterAddresser().get_cog_key(output_key) + generate_cog_spy = mocker.patch( + "wps_sfms.publish.generate_web_optimized_cog", + return_value=expected_cog_key, + ) + + mock_s3_client = MagicMock() + mock_s3_client.put_object = AsyncMock() + + gdal_ds = create_test_dataset( + "temperature_20240704.tif", + 1, + 1, + (-1.0, 1.0, -1.0, 1.0), + 4326, + fill_value=12.5, + no_data_value=-9999.0, + ) + + with WPSDataset(ds_path=None, ds=gdal_ds) as dataset: + published = await publish_dataset(mock_s3_client, dataset, output_key) + + mock_s3_client.put_object.assert_awaited_once() + assert mock_s3_client.put_object.await_args.kwargs["key"] == output_key + assert isinstance(mock_s3_client.put_object.await_args.kwargs["body"], bytes) + assert mock_s3_client.put_object.await_args.kwargs["body"] + + generate_cog_spy.assert_called_once() + assert generate_cog_spy.call_args.kwargs["output_path"] == expected_cog_key + assert generate_cog_spy.call_args.kwargs["input_path"].endswith("temperature_20240704.tif") + + assert published.output_key == output_key + assert published.cog_key == expected_cog_key + + +@pytest.mark.anyio +async def test_publish_dataset_can_skip_cog_generation(mocker): + output_key = "sfms_ng/actual/2024/07/04/temperature_20240704.tif" + generate_cog_spy = mocker.patch("wps_sfms.publish.generate_web_optimized_cog") + + mock_s3_client = MagicMock() + mock_s3_client.put_object = AsyncMock() + + gdal_ds = create_test_dataset( + "temperature_20240704.tif", + 1, + 1, + (-1.0, 1.0, -1.0, 1.0), + 4326, + fill_value=12.5, + no_data_value=-9999.0, + ) + + with WPSDataset(ds_path=None, ds=gdal_ds) as dataset: + published = await publish_dataset( + mock_s3_client, + dataset, + output_key, + generate_cog=False, + ) + + mock_s3_client.put_object.assert_awaited_once() + generate_cog_spy.assert_not_called() + assert published.output_key == output_key + assert published.cog_key is None diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_rh_interpolation.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_rh_interpolation.py index 6bdf2ad63..ab7d558ac 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/tests/test_rh_interpolation.py +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_rh_interpolation.py @@ -7,7 +7,7 @@ import pytest from osgeo import gdal from wps_shared.schemas.sfms import SFMSDailyActual -from wps_sfms.interpolation.source import StationDewPointSource +from wps_sfms.interpolation.field import build_dewpoint_field, compute_rh from wps_sfms.processors.relative_humidity import RHInterpolator from wps_sfms.tests.conftest import create_test_raster @@ -34,21 +34,21 @@ def test_dewpoint_equals_temp_gives_100_percent(self): """When dew point equals temperature, RH should be 100%.""" temp = np.array([20.0, 10.0, 0.0], dtype=np.float32) dewpoint = np.array([20.0, 10.0, 0.0], dtype=np.float32) - rh = StationDewPointSource.compute_rh(temp, dewpoint) + rh = compute_rh(temp, dewpoint) np.testing.assert_allclose(rh, 100.0, atol=0.01) def test_lower_dewpoint_gives_lower_rh(self): """Lower dew point relative to temperature should give lower RH.""" temp = np.array([20.0, 20.0, 20.0], dtype=np.float32) dewpoint = np.array([20.0, 15.0, 10.0], dtype=np.float32) - rh = StationDewPointSource.compute_rh(temp, dewpoint) + rh = compute_rh(temp, dewpoint) assert rh[0] > rh[1] > rh[2] def test_rh_clamped_to_0_100(self): """RH should be clamped between 0 and 100.""" temp = np.array([20.0], dtype=np.float32) dewpoint = np.array([-50.0], dtype=np.float32) - rh = StationDewPointSource.compute_rh(temp, dewpoint) + rh = compute_rh(temp, dewpoint) assert np.all(rh >= 0.0) assert np.all(rh <= 100.0) @@ -57,14 +57,14 @@ def test_known_values(self): # At 20°C with dewpoint of 10°C: e_s(10)/e_s(20) ≈ 0.5258 → ~52.58% temp = np.array([20.0], dtype=np.float32) dewpoint = np.array([10.0], dtype=np.float32) - rh = StationDewPointSource.compute_rh(temp, dewpoint) + rh = compute_rh(temp, dewpoint) np.testing.assert_allclose(rh, 52.58, atol=0.1) def test_output_dtype_is_float32(self): """Output should be float32.""" temp = np.array([20.0], dtype=np.float32) dewpoint = np.array([15.0], dtype=np.float32) - rh = StationDewPointSource.compute_rh(temp, dewpoint) + rh = compute_rh(temp, dewpoint) assert rh.dtype == np.float32 @@ -92,11 +92,14 @@ def test_interpolate_basic_success(self): dewpoints=[12.0, 11.0], elevations=[100.0, 200.0], ) - dewpoint_source = StationDewPointSource(actuals) + dewpoint_field = build_dewpoint_field(actuals) dataset = RHInterpolator( - mask_path=mask_path, dem_path=dem_path, temp_raster_path=temp_raster_path - ).interpolate(dewpoint_source, ref_path) + mask_path=mask_path, + dem_path=dem_path, + temp_raster_path=temp_raster_path, + field=dewpoint_field, + ).interpolate(ref_path) data = dataset.ds.GetRasterBand(1).ReadAsArray() nodata = dataset.ds.GetRasterBand(1).GetNoDataValue() @@ -135,11 +138,14 @@ def test_interpolate_skips_masked_cells(self): actuals = create_test_actuals( lats=[49.05], lons=[-123.05], dewpoints=[12.0], elevations=[100.0] ) - dewpoint_source = StationDewPointSource(actuals) + dewpoint_field = build_dewpoint_field(actuals) dataset = RHInterpolator( - mask_path=mask_path, dem_path=dem_path, temp_raster_path=temp_raster_path - ).interpolate(dewpoint_source, ref_path) + mask_path=mask_path, + dem_path=dem_path, + temp_raster_path=temp_raster_path, + field=dewpoint_field, + ).interpolate(ref_path) data = dataset.ds.GetRasterBand(1).ReadAsArray() nodata = dataset.ds.GetRasterBand(1).GetNoDataValue() @@ -171,11 +177,14 @@ def test_output_preserves_reference_properties(self): actuals = create_test_actuals( lats=[49.05], lons=[-123.05], dewpoints=[12.0], elevations=[100.0] ) - dewpoint_source = StationDewPointSource(actuals) + dewpoint_field = build_dewpoint_field(actuals) dataset = RHInterpolator( - mask_path=mask_path, dem_path=dem_path, temp_raster_path=temp_raster_path - ).interpolate(dewpoint_source, ref_path) + mask_path=mask_path, + dem_path=dem_path, + temp_raster_path=temp_raster_path, + field=dewpoint_field, + ).interpolate(ref_path) ref_ds = gdal.Open(ref_path) assert dataset.ds.RasterXSize == ref_ds.RasterXSize @@ -211,11 +220,14 @@ def test_rh_values_respond_to_elevation(self): actuals = create_test_actuals( lats=[49.05], lons=[-123.05], dewpoints=[12.0], elevations=[100.0] ) - dewpoint_source = StationDewPointSource(actuals) + dewpoint_field = build_dewpoint_field(actuals) dataset = RHInterpolator( - mask_path=mask_path, dem_path=dem_path, temp_raster_path=temp_raster_path - ).interpolate(dewpoint_source, ref_path) + mask_path=mask_path, + dem_path=dem_path, + temp_raster_path=temp_raster_path, + field=dewpoint_field, + ).interpolate(ref_path) data = dataset.ds.GetRasterBand(1).ReadAsArray() nodata = dataset.ds.GetRasterBand(1).GetNoDataValue() @@ -252,12 +264,15 @@ def test_interpolate_raises_when_no_valid_stations(self): actuals = [ SFMSDailyActual(code=1, lat=49.05, lon=-123.05, elevation=100.0, dewpoint=None) ] - dewpoint_source = StationDewPointSource(actuals) + dewpoint_field = build_dewpoint_field(actuals) with pytest.raises(RuntimeError, match="No pixels were successfully interpolated"): RHInterpolator( - mask_path=mask_path, dem_path=dem_path, temp_raster_path=temp_raster_path - ).interpolate(dewpoint_source, ref_path) + mask_path=mask_path, + dem_path=dem_path, + temp_raster_path=temp_raster_path, + field=dewpoint_field, + ).interpolate(ref_path) finally: gdal.Unlink(ref_path) gdal.Unlink(dem_path) diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_sfmsng_raster_addresser.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_sfmsng_raster_addresser.py index cccefa70f..0e15ea4d7 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/tests/test_sfmsng_raster_addresser.py +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_sfmsng_raster_addresser.py @@ -33,6 +33,10 @@ class TestGetActualWeatherKey: SFMSInterpolatedWeatherParameter.WIND_SPEED, "sfms_ng/actual/2024/04/15/wind_speed_20240415.tif", ), + ( + SFMSInterpolatedWeatherParameter.WIND_DIRECTION, + "sfms_ng/actual/2024/04/15/wind_direction_20240415.tif", + ), ( SFMSInterpolatedWeatherParameter.PRECIP, "sfms_ng/actual/2024/04/15/precipitation_20240415.tif", @@ -75,35 +79,76 @@ def test_non_utc_raises(self, addresser: SFMSNGRasterAddresser): class TestGetActualFwiInputs: - @pytest.mark.parametrize("fwi_param", [FWIParameter.FFMC, FWIParameter.DMC, FWIParameter.DC]) - def test_all_params(self, addresser: SFMSNGRasterAddresser, fwi_param: FWIParameter): + @pytest.mark.parametrize( + "fwi_param,index_key_checks", + [ + (FWIParameter.FFMC, {FWIParameter.FFMC: "2024/04/14/ffmc_20240414.tif"}), + (FWIParameter.DMC, {FWIParameter.DMC: "2024/04/14/dmc_20240414.tif"}), + (FWIParameter.DC, {FWIParameter.DC: "2024/04/14/dc_20240414.tif"}), + (FWIParameter.ISI, {FWIParameter.FFMC: "2024/04/15/ffmc_20240415.tif"}), + ( + FWIParameter.BUI, + { + FWIParameter.DMC: "2024/04/15/dmc_20240415.tif", + FWIParameter.DC: "2024/04/15/dc_20240415.tif", + }, + ), + ( + FWIParameter.FWI, + { + FWIParameter.ISI: "2024/04/15/isi_20240415.tif", + FWIParameter.BUI: "2024/04/15/bui_20240415.tif", + }, + ), + ], + ) + def test_all_params( + self, + addresser: SFMSNGRasterAddresser, + fwi_param: FWIParameter, + index_key_checks: dict[FWIParameter, str], + ): s3 = addresser.s3_prefix p = fwi_param.value result = addresser.get_actual_fwi_inputs(TEST_DATETIME, fwi_param) - assert result.temp_key == f"{s3}/sfms_ng/actual/2024/04/15/temperature_20240415.tif" - assert result.rh_key == f"{s3}/sfms_ng/actual/2024/04/15/relative_humidity_20240415.tif" - assert result.precip_key == f"{s3}/sfms_ng/actual/2024/04/15/precipitation_20240415.tif" - assert result.prev_fwi_key == f"{s3}/sfms_ng/actual/2024/04/14/{p}_20240414.tif" + assert ( + result.weather_keys[SFMSInterpolatedWeatherParameter.TEMP] + == f"{s3}/sfms_ng/actual/2024/04/15/temperature_20240415.tif" + ) + assert ( + result.weather_keys[SFMSInterpolatedWeatherParameter.RH] + == f"{s3}/sfms_ng/actual/2024/04/15/relative_humidity_20240415.tif" + ) + assert ( + result.weather_keys[SFMSInterpolatedWeatherParameter.PRECIP] + == f"{s3}/sfms_ng/actual/2024/04/15/precipitation_20240415.tif" + ) + assert ( + result.weather_keys[SFMSInterpolatedWeatherParameter.WIND_SPEED] + == f"{s3}/sfms_ng/actual/2024/04/15/wind_speed_20240415.tif" + ) + for index_param, expected_suffix in index_key_checks.items(): + assert result.index_keys[index_param] == f"{s3}/sfms_ng/actual/{expected_suffix}" assert result.output_key == f"sfms_ng/actual/2024/04/15/{p}_20240415.tif" - assert result.cog_key == f"{s3}/sfms_ng/actual/2024/04/15/{p}_20240415_cog.tif" assert result.run_type == RunType.ACTUAL def test_gdal_prefix_on_inputs_not_output(self, addresser: SFMSNGRasterAddresser): s3 = addresser.s3_prefix result = addresser.get_actual_fwi_inputs(TEST_DATETIME, FWIParameter.DMC) - assert result.temp_key.startswith(s3) - assert result.rh_key.startswith(s3) - assert result.precip_key.startswith(s3) - assert result.prev_fwi_key.startswith(s3) + for weather_key in result.weather_keys.values(): + assert weather_key.startswith(s3) + for index_key in result.index_keys.values(): + assert index_key.startswith(s3) assert not result.output_key.startswith(s3) def test_prev_fwi_key_uses_yesterday(self, addresser: SFMSNGRasterAddresser): result = addresser.get_actual_fwi_inputs(TEST_DATETIME, FWIParameter.DC) - assert "2024/04/14" in result.prev_fwi_key - assert "2024/04/15" not in result.prev_fwi_key + dc_dependency_key = result.index_keys[FWIParameter.DC] + assert "2024/04/14" in dc_dependency_key + assert "2024/04/15" not in dc_dependency_key def test_output_key_uses_actual_run_type(self, addresser: SFMSNGRasterAddresser): result = addresser.get_actual_fwi_inputs(TEST_DATETIME, FWIParameter.DMC) diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_source.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_source.py deleted file mode 100644 index 230d7ae55..000000000 --- a/backend/packages/wps-sfms/src/wps_sfms/tests/test_source.py +++ /dev/null @@ -1,32 +0,0 @@ -import pytest -from wps_shared.schemas.sfms import SFMSDailyActual -from wps_sfms.interpolation.source import StationActualSource - - -class TestStationActualSource: - """Tests for StationActualSource attribute validation.""" - - def test_valid_attribute_constructs_successfully(self): - """Test that a valid attribute name constructs without error.""" - actuals = [SFMSDailyActual(code=1, lat=49.0, lon=-123.0, precipitation=5.0)] - source = StationActualSource("precipitation", actuals) - lats, _, values = source.get_interpolation_data() - assert len(lats) == 1 - assert values[0] == pytest.approx(5.0) - - def test_invalid_attribute_raises_value_error(self): - """Test that an unknown attribute name raises ValueError at construction.""" - actuals = [SFMSDailyActual(code=1, lat=49.0, lon=-123.0)] - with pytest.raises(ValueError, match="Unknown attribute"): - StationActualSource("not_a_real_field", actuals) - - def test_none_values_are_excluded(self): - """Test that stations with None for the target attribute are excluded.""" - actuals = [ - SFMSDailyActual(code=1, lat=49.0, lon=-123.0, precipitation=3.0), - SFMSDailyActual(code=2, lat=49.1, lon=-123.1, precipitation=None), - ] - source = StationActualSource("precipitation", actuals) - lats, _, values = source.get_interpolation_data() - assert len(lats) == 1 - assert values[0] == pytest.approx(3.0) diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_station_dewpoint_source.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_station_dewpoint_source.py deleted file mode 100644 index 56c4df68d..000000000 --- a/backend/packages/wps-sfms/src/wps_sfms/tests/test_station_dewpoint_source.py +++ /dev/null @@ -1,121 +0,0 @@ -""" -Unit tests for StationDewPointSource. -""" - -from numpy.testing import assert_allclose - -from wps_shared.schemas.sfms import SFMSDailyActual -from wps_sfms.interpolation.source import ( - DEW_POINT_LAPSE_RATE, - LapseRateAdjustedSource, - StationDewPointSource, -) - - -def _make_actual(code, lat, lon, elevation, dewpoint): - return SFMSDailyActual( - code=code, - lat=lat, - lon=lon, - elevation=elevation, - dewpoint=dewpoint, - ) - - -def test_dewpoint_values_read_directly(): - """Dew point values should be read directly from actuals.""" - actuals = [ - _make_actual(1, 49.0, -123.0, 100.0, 12.0), - _make_actual(2, 49.1, -123.1, 200.0, 11.0), - ] - source = StationDewPointSource(actuals) - _, _, _, dewpoints = source.get_station_arrays(only_valid=True) - - assert_allclose(dewpoints[0], 12.0, atol=1e-4) - assert_allclose(dewpoints[1], 11.0, atol=1e-4) - - -def test_missing_dewpoint_excluded(): - """Stations with missing dewpoint should be filtered out.""" - actuals = [ - _make_actual(1, 49.0, -123.0, 100.0, 12.0), - _make_actual(2, 49.1, -123.1, 200.0, None), # Missing dewpoint - ] - source = StationDewPointSource(actuals) - lats, _, _, dewpoints = source.get_station_arrays(only_valid=True) - assert len(lats) == 1 - assert len(dewpoints) == 1 - - -def test_missing_elevation_excluded(): - """Stations with missing elevation should be filtered out.""" - actuals = [ - _make_actual(1, 49.0, -123.0, None, 12.0), # Missing elevation - _make_actual(2, 49.1, -123.1, 200.0, 11.0), - ] - source = StationDewPointSource(actuals) - lats, _, _, _ = source.get_station_arrays(only_valid=True) - assert len(lats) == 1 - - -def test_get_interpolation_data_returns_sea_level_dewpoints(): - """get_interpolation_data should return sea-level adjusted dew points.""" - actuals = [ - _make_actual(1, 49.0, -123.0, 500.0, 12.0), - ] - source = StationDewPointSource(actuals) - _, _, sea_level_td = source.get_interpolation_data(lapse_rate=DEW_POINT_LAPSE_RATE) - - _, _, _, raw_dewpoints = source.get_station_arrays(only_valid=True) - _, _, elevs, _ = source.get_station_arrays(only_valid=True) - - # Sea level dew point should be warmer than actual (positive elevation) - expected_sea = LapseRateAdjustedSource.compute_sea_level_values( - raw_dewpoints, elevs, DEW_POINT_LAPSE_RATE - ) - assert_allclose(sea_level_td, expected_sea, atol=1e-5) - assert sea_level_td[0] > raw_dewpoints[0] - - -def test_empty_actuals(): - """Empty actuals should return empty arrays.""" - source = StationDewPointSource([]) - lats, lons, sea_td = source.get_interpolation_data() - assert len(lats) == 0 - assert len(lons) == 0 - assert len(sea_td) == 0 - - -def test_all_invalid_returns_empty(): - """If all stations have missing data, should return empty arrays.""" - actuals = [ - _make_actual(1, 49.0, -123.0, None, None), - _make_actual(2, 49.1, -123.1, None, None), - ] - source = StationDewPointSource(actuals) - lats, _, _ = source.get_interpolation_data() - assert len(lats) == 0 - - -def test_station_count(): - """get_station_count should return total stations, not just valid ones.""" - actuals = [ - _make_actual(1, 49.0, -123.0, 100.0, 12.0), - _make_actual(2, 49.1, -123.1, None, None), # Invalid - ] - source = StationDewPointSource(actuals) - assert source.get_station_count() == 2 - - -def test_round_trip_lapse_rate(): - """Sea-level adjusted dew point adjusted back to station elevation should match original.""" - actuals = [ - _make_actual(1, 49.0, -123.0, 500.0, 19.0), - _make_actual(2, 49.1, -123.1, 1000.0, 0.0), - ] - source = StationDewPointSource(actuals) - _, _, elevs, original_td = source.get_station_arrays(only_valid=True) - _, _, sea_td = source.get_interpolation_data(lapse_rate=DEW_POINT_LAPSE_RATE) - - back = LapseRateAdjustedSource.compute_adjusted_values(sea_td, elevs, DEW_POINT_LAPSE_RATE) - assert_allclose(back, original_td, atol=1e-4) diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_station_temperature_source.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_station_temperature_source.py index f4d32f9c2..69e457352 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/tests/test_station_temperature_source.py +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_station_temperature_source.py @@ -1,31 +1,31 @@ import numpy as np from numpy.testing import assert_allclose -from wps_sfms.interpolation.source import LAPSE_RATE, LapseRateAdjustedSource +from wps_sfms.interpolation.field import LAPSE_RATE, compute_adjusted_values, compute_sea_level_values def test_zero_elevation_identity(): temps = np.array([-5.0, 0.0, 12.3], dtype=np.float32) elevs = np.zeros_like(temps) - out = LapseRateAdjustedSource.compute_sea_level_values(temps, elevs, LAPSE_RATE) + out = compute_sea_level_values(temps, elevs, LAPSE_RATE) assert_allclose(out, temps, atol=0) # round-trip with actual temps: - rt = LapseRateAdjustedSource.compute_adjusted_values(out, elevs, LAPSE_RATE) + rt = compute_adjusted_values(out, elevs, LAPSE_RATE) assert_allclose(rt, temps, atol=0) def test_positive_elevation_cools_when_applying_to_terrain(): sea = np.array([20.0, 15.0, 10.0], dtype=np.float32) elev = np.array([1000.0, 1500.0, 2000.0], dtype=np.float32) - out = LapseRateAdjustedSource.compute_adjusted_values(sea, elev, LAPSE_RATE) + out = compute_adjusted_values(sea, elev, LAPSE_RATE) assert np.all(out < sea) def test_negative_elevation_warms_when_applying_to_terrain(): sea = np.array([10.0], dtype=np.float32) elev = np.array([-100.0], dtype=np.float32) - out = LapseRateAdjustedSource.compute_adjusted_values(sea, elev, LAPSE_RATE) + out = compute_adjusted_values(sea, elev, LAPSE_RATE) assert out[0] > sea[0] @@ -33,7 +33,7 @@ def test_broadcasting_shapes(): # sea: (2,1), elev: (3,) => result (2,3) temps = np.array([[20.0], [15.0]], dtype=np.float32) elevs = np.array([0.0, 1000.0, 2000.0], dtype=np.float32) - out = LapseRateAdjustedSource.compute_adjusted_values(temps, elevs, LAPSE_RATE) + out = compute_adjusted_values(temps, elevs, LAPSE_RATE) expected = temps - elevs * np.float32(LAPSE_RATE) assert_allclose(out, expected, atol=1e-6) assert out.shape == (2, 3) @@ -44,10 +44,10 @@ def test_dtype_is_float32_even_if_inputs_float64(): temps = np.array([20.0, 10.0], dtype=np.float64) elevs = np.array([0.0, 1000.0], dtype=np.float64) - sea = LapseRateAdjustedSource.compute_sea_level_values(temps, elevs, LAPSE_RATE) + sea = compute_sea_level_values(temps, elevs, LAPSE_RATE) assert sea.dtype == np.float32 - act = LapseRateAdjustedSource.compute_adjusted_values(sea, elevs.astype(np.float32), LAPSE_RATE) + act = compute_adjusted_values(sea, elevs.astype(np.float32), LAPSE_RATE) assert act.dtype == np.float32 @@ -55,6 +55,6 @@ def test_round_trip_sea_then_actual_returns_original(): # compute_actual(compute_sea(temps, elev), elev) == temps temps = np.array([5.0, 10.0, -3.0, 0.0], dtype=np.float32) elevs = np.array([0.0, 100.0, 500.0, 2000.0], dtype=np.float32) - sea = LapseRateAdjustedSource.compute_sea_level_values(temps, elevs, LAPSE_RATE) - back = LapseRateAdjustedSource.compute_adjusted_values(sea, elevs, LAPSE_RATE) + sea = compute_sea_level_values(temps, elevs, LAPSE_RATE) + back = compute_adjusted_values(sea, elevs, LAPSE_RATE) assert_allclose(back, temps, rtol=0, atol=1e-6) diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_temperature_interpolation.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_temperature_interpolation.py index b6a30b246..74bfcd510 100644 --- a/backend/packages/wps-sfms/src/wps_sfms/tests/test_temperature_interpolation.py +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_temperature_interpolation.py @@ -7,7 +7,7 @@ import pytest from osgeo import gdal from wps_shared.schemas.sfms import SFMSDailyActual -from wps_sfms.interpolation.source import StationTemperatureSource +from wps_sfms.interpolation.field import build_temperature_field from wps_sfms.processors.temperature import TemperatureInterpolator from wps_sfms.tests.conftest import create_test_raster @@ -52,11 +52,11 @@ def test_interpolate_basic_success(self): temps=[15.0, 12.0], elevations=[100.0, 200.0], ) - temperature_source = StationTemperatureSource(actuals) + temperature_field = build_temperature_field(actuals) dataset = TemperatureInterpolator( - mask_path=mask_path, dem_path=dem_path - ).interpolate(temperature_source, ref_path) + mask_path=mask_path, dem_path=dem_path, field=temperature_field + ).interpolate(ref_path) data = dataset.ds.GetRasterBand(1).ReadAsArray() nodata = dataset.ds.GetRasterBand(1).GetNoDataValue() @@ -88,11 +88,11 @@ def test_interpolate_skips_masked_cells(self): actuals = create_test_actuals( lats=[49.05], lons=[-123.05], temps=[15.0], elevations=[100.0] ) - temperature_source = StationTemperatureSource(actuals) + temperature_field = build_temperature_field(actuals) dataset = TemperatureInterpolator( - mask_path=mask_path, dem_path=dem_path - ).interpolate(temperature_source, ref_path) + mask_path=mask_path, dem_path=dem_path, field=temperature_field + ).interpolate(ref_path) data = dataset.ds.GetRasterBand(1).ReadAsArray() nodata = dataset.ds.GetRasterBand(1).GetNoDataValue() @@ -128,11 +128,11 @@ def test_interpolate_with_elevation_adjustment(self): actuals = create_test_actuals( lats=[49.05], lons=[-123.05], temps=[15.0], elevations=[100.0] ) - temperature_source = StationTemperatureSource(actuals) + temperature_field = build_temperature_field(actuals) dataset = TemperatureInterpolator( - mask_path=mask_path, dem_path=dem_path - ).interpolate(temperature_source, ref_path) + mask_path=mask_path, dem_path=dem_path, field=temperature_field + ).interpolate(ref_path) data = dataset.ds.GetRasterBand(1).ReadAsArray() @@ -167,11 +167,11 @@ def test_output_preserves_reference_properties(self): actuals = create_test_actuals( lats=[49.05], lons=[-123.05], temps=[15.0], elevations=[100.0] ) - temperature_source = StationTemperatureSource(actuals) + temperature_field = build_temperature_field(actuals) dataset = TemperatureInterpolator( - mask_path=mask_path, dem_path=dem_path - ).interpolate(temperature_source, ref_path) + mask_path=mask_path, dem_path=dem_path, field=temperature_field + ).interpolate(ref_path) ref_ds = gdal.Open(ref_path) assert dataset.ds.RasterXSize == ref_ds.RasterXSize @@ -200,12 +200,12 @@ def test_interpolate_raises_when_no_valid_stations(self): actuals = [ SFMSDailyActual(code=1, lat=49.05, lon=-123.05, elevation=100.0, temperature=None) ] - temperature_source = StationTemperatureSource(actuals) + temperature_field = build_temperature_field(actuals) with pytest.raises(RuntimeError, match="No pixels were successfully interpolated"): TemperatureInterpolator( - mask_path=mask_path, dem_path=dem_path - ).interpolate(temperature_source, ref_path) + mask_path=mask_path, dem_path=dem_path, field=temperature_field + ).interpolate(ref_path) finally: gdal.Unlink(ref_path) gdal.Unlink(dem_path) diff --git a/backend/packages/wps-sfms/src/wps_sfms/tests/test_wind_interpolation.py b/backend/packages/wps-sfms/src/wps_sfms/tests/test_wind_interpolation.py new file mode 100644 index 000000000..2ba184a3c --- /dev/null +++ b/backend/packages/wps-sfms/src/wps_sfms/tests/test_wind_interpolation.py @@ -0,0 +1,162 @@ +import uuid + +import numpy as np +import pytest +from osgeo import gdal + +import wps_sfms.processors.wind as wind_module +from wps_shared.schemas.sfms import SFMSDailyActual +from wps_sfms.interpolation.field import build_wind_vector_field +from wps_sfms.processors.idw import ValidPixelIDWResult +from wps_sfms.processors.wind import WindDirectionInterpolator +from wps_sfms.tests.conftest import create_test_raster + + +class TestWindDirectionInterpolator: + def test_compute_direction_from_uv_matches_legacy_special_cases(self): + u = np.array([-1.0, 1.0, 0.0], dtype=np.float32) + v = np.array([0.0, 0.0, 0.0], dtype=np.float32) + + direction = WindDirectionInterpolator.compute_direction_from_uv(u, v) + + np.testing.assert_allclose(direction, np.array([90.0, 270.0, 0.0], dtype=np.float32)) + + def test_compute_direction_from_uv_uses_arctan_branch_when_v_nonzero(self): + # Four quadrants with v != 0 should follow: deg(atan2(u, v)) + 180 + u = np.array([1.0, -1.0, 1.0, -1.0], dtype=np.float32) + v = np.array([1.0, 1.0, -1.0, -1.0], dtype=np.float32) + + direction = WindDirectionInterpolator.compute_direction_from_uv(u, v) + + np.testing.assert_allclose( + direction, + np.array([225.0, 135.0, 315.0, 45.0], dtype=np.float32), + atol=1e-5, + ) + + def test_compute_direction_from_uv_mixed_zero_and_nonzero_v(self): + # Mixed inputs should apply atan2 where v != 0 and legacy overrides where v == 0. + u = np.array([1.0, -1.0, 0.0], dtype=np.float32) + v = np.array([1.0, 0.0, 0.5], dtype=np.float32) + + direction = WindDirectionInterpolator.compute_direction_from_uv(u, v) + + np.testing.assert_allclose( + direction, + np.array([225.0, 90.0, 180.0], dtype=np.float32), + atol=1e-5, + ) + + def test_interpolate_basic_success(self): + test_id = uuid.uuid4().hex + ref_path = f"/vsimem/reference_{test_id}.tif" + mask_path = f"/vsimem/mask_{test_id}.tif" + + try: + extent = (-123.1, -123.0, 49.0, 49.1) + create_test_raster(ref_path, 10, 10, extent, fill_value=1.0) + create_test_raster(mask_path, 10, 10, extent, fill_value=1.0) + + actuals = [ + SFMSDailyActual( + code=100, lat=49.05, lon=-123.05, wind_speed=10.0, wind_direction=90.0 + ), + SFMSDailyActual( + code=101, lat=49.08, lon=-123.02, wind_speed=8.0, wind_direction=180.0 + ), + ] + field = build_wind_vector_field(actuals) + + dataset = WindDirectionInterpolator(mask_path=mask_path, field=field).interpolate( + ref_path + ) + data = dataset.ds.GetRasterBand(1).ReadAsArray() + nodata = dataset.ds.GetRasterBand(1).GetNoDataValue() + valid = data[data != nodata] + + assert valid.size > 0 + assert np.all(valid >= 0.0) + assert np.all(valid <= 360.0) + finally: + gdal.Unlink(ref_path) + gdal.Unlink(mask_path) + + def test_interpolate_raises_without_paired_stations(self): + test_id = uuid.uuid4().hex + ref_path = f"/vsimem/reference_{test_id}.tif" + mask_path = f"/vsimem/mask_{test_id}.tif" + + try: + extent = (-123.1, -123.0, 49.0, 49.1) + create_test_raster(ref_path, 5, 5, extent, fill_value=1.0) + create_test_raster(mask_path, 5, 5, extent, fill_value=1.0) + + actuals = [ + SFMSDailyActual( + code=100, lat=49.05, lon=-123.05, wind_speed=10.0, wind_direction=None + ) + ] + field = build_wind_vector_field(actuals) + + with pytest.raises(RuntimeError, match="No pixels were successfully interpolated"): + WindDirectionInterpolator(mask_path=mask_path, field=field).interpolate(ref_path) + finally: + gdal.Unlink(ref_path) + gdal.Unlink(mask_path) + + @pytest.mark.parametrize("successful_component_label", ["wind-u component", "wind-v component"]) + def test_interpolate_raises_when_only_one_component_succeeds( + self, monkeypatch, successful_component_label + ): + test_id = uuid.uuid4().hex + ref_path = f"/vsimem/reference_{test_id}.tif" + mask_path = f"/vsimem/mask_{test_id}.tif" + + try: + extent = (-123.1, -123.0, 49.0, 49.1) + create_test_raster(ref_path, 5, 5, extent, fill_value=1.0) + create_test_raster(mask_path, 5, 5, extent, fill_value=1.0) + + actuals = [ + SFMSDailyActual( + code=100, lat=49.05, lon=-123.05, wind_speed=10.0, wind_direction=90.0 + ), + SFMSDailyActual( + code=101, lat=49.08, lon=-123.02, wind_speed=8.0, wind_direction=180.0 + ), + ] + field = build_wind_vector_field(actuals) + + def _fake_idw_on_valid_pixels(**kwargs): + valid_yi = kwargs["valid_yi"] + valid_xi = kwargs["valid_xi"] + total_pixels = kwargs["total_pixels"] + label = kwargs["label"] + n = len(valid_yi) + + if label == successful_component_label: + interpolated_values = np.full(n, 4.0, dtype=np.float32) + succeeded_mask = np.ones(n, dtype=bool) + else: + interpolated_values = np.full(n, np.nan, dtype=np.float32) + succeeded_mask = np.zeros(n, dtype=bool) + + return ValidPixelIDWResult( + interpolated_values=interpolated_values, + succeeded_mask=succeeded_mask, + rows=valid_yi[succeeded_mask], + cols=valid_xi[succeeded_mask], + values=interpolated_values[succeeded_mask].astype(np.float32, copy=False), + total_pixels=total_pixels, + interpolated_count=int(np.sum(succeeded_mask)), + failed_interpolation_count=n - int(np.sum(succeeded_mask)), + skipped_nodata_count=total_pixels - n, + ) + + monkeypatch.setattr(wind_module, "idw_on_valid_pixels", _fake_idw_on_valid_pixels) + + with pytest.raises(RuntimeError, match="No pixels were successfully interpolated"): + WindDirectionInterpolator(mask_path=mask_path, field=field).interpolate(ref_path) + finally: + gdal.Unlink(ref_path) + gdal.Unlink(mask_path) diff --git a/backend/packages/wps-shared/src/wps_shared/db/models/sfms_run.py b/backend/packages/wps-shared/src/wps_shared/db/models/sfms_run.py index 1375c14e9..0b0350405 100644 --- a/backend/packages/wps-shared/src/wps_shared/db/models/sfms_run.py +++ b/backend/packages/wps-shared/src/wps_shared/db/models/sfms_run.py @@ -15,6 +15,8 @@ class SFMSRunLogJobName(str, enum.Enum): TEMPERATURE_INTERPOLATION = "temperature_interpolation" PRECIPITATION_INTERPOLATION = "precipitation_interpolation" RH_INTERPOLATION = "rh_interpolation" + WIND_SPEED_INTERPOLATION = "wind_speed_interpolation" + WIND_DIRECTION_INTERPOLATION = "wind_direction_interpolation" FFMC_INTERPOLATION = "ffmc_interpolation" DMC_INTERPOLATION = "dmc_interpolation" DC_INTERPOLATION = "dc_interpolation" diff --git a/backend/packages/wps-shared/src/wps_shared/geospatial/wps_dataset.py b/backend/packages/wps-shared/src/wps_shared/geospatial/wps_dataset.py index 0b06e4346..e3e0eedf3 100644 --- a/backend/packages/wps-shared/src/wps_shared/geospatial/wps_dataset.py +++ b/backend/packages/wps-shared/src/wps_shared/geospatial/wps_dataset.py @@ -5,7 +5,7 @@ import numpy as np import io -from wps_shared.geospatial.geospatial import GDALResamplingMethod +from wps_shared.geospatial.geospatial import GDALResamplingMethod, rasters_match gdal.UseExceptions() @@ -390,24 +390,8 @@ def apply_mask(self, mask_ds: "WPSDataset") -> np.ndarray: :return: Boolean array where True = valid, False = masked :raises ValueError: If the mask grid does not match this dataset's grid """ - mismatches = [] - if self.ds.RasterXSize != mask_ds.ds.RasterXSize or self.ds.RasterYSize != mask_ds.ds.RasterYSize: - mismatches.append( - f"size: reference=({self.ds.RasterXSize}x{self.ds.RasterYSize}), " - f"mask=({mask_ds.ds.RasterXSize}x{mask_ds.ds.RasterYSize})" - ) - if self.ds.GetGeoTransform() != mask_ds.ds.GetGeoTransform(): - mismatches.append( - f"geotransform: reference={self.ds.GetGeoTransform()}, " - f"mask={mask_ds.ds.GetGeoTransform()}" - ) - if self.ds.GetProjection() != mask_ds.ds.GetProjection(): - mismatches.append("projection differs") - - if mismatches: - raise ValueError( - "Mask grid does not match reference grid: " + "; ".join(mismatches) - ) + if not rasters_match(self.ds, mask_ds.ds): + raise ValueError("Mask grid does not match reference grid") mask_band: gdal.Band = mask_ds.ds.GetRasterBand(1) mask_data = mask_band.ReadAsArray() diff --git a/backend/packages/wps-shared/src/wps_shared/logging.json b/backend/packages/wps-shared/src/wps_shared/logging.json index 7709b5ee0..376cb7887 100644 --- a/backend/packages/wps-shared/src/wps_shared/logging.json +++ b/backend/packages/wps-shared/src/wps_shared/logging.json @@ -42,6 +42,13 @@ ], "propagate": false }, + "wps_sfms": { + "level": "INFO", + "handlers": [ + "console" + ], + "propagate": false + }, "__main__": { "level": "INFO", "handlers": [ diff --git a/backend/packages/wps-shared/src/wps_shared/sfms/raster_addresser.py b/backend/packages/wps-shared/src/wps_shared/sfms/raster_addresser.py index c3ca4e868..033d1548f 100644 --- a/backend/packages/wps-shared/src/wps_shared/sfms/raster_addresser.py +++ b/backend/packages/wps-shared/src/wps_shared/sfms/raster_addresser.py @@ -17,6 +17,7 @@ class SFMSInterpolatedWeatherParameter(enum.Enum): TEMP = "temperature" RH = "relative_humidity" WIND_SPEED = "wind_speed" + WIND_DIRECTION = "wind_direction" PRECIP = "precipitation" diff --git a/backend/packages/wps-shared/src/wps_shared/tests/geospatial/dataset_common.py b/backend/packages/wps-shared/src/wps_shared/tests/geospatial/dataset_common.py index 2192b353f..291c8bbe5 100644 --- a/backend/packages/wps-shared/src/wps_shared/tests/geospatial/dataset_common.py +++ b/backend/packages/wps-shared/src/wps_shared/tests/geospatial/dataset_common.py @@ -6,7 +6,16 @@ from wps_shared.geospatial.wps_dataset import WPSDataset -def create_test_dataset(filename, width, height, extent, projection, data_type=gdal.GDT_Float32, fill_value=None, no_data_value=None) -> gdal.Dataset: +def create_test_dataset( + filename, + width, + height, + extent, + projection, + data_type=gdal.GDT_Float32, + fill_value=None, + no_data_value=None, +) -> gdal.Dataset: """ Create a test GDAL dataset. """ @@ -42,7 +51,9 @@ def create_test_dataset(filename, width, height, extent, projection, data_type=g def create_mock_gdal_dataset(): extent = (-1, 1, -1, 1) # xmin, xmax, ymin, ymax - return create_test_dataset(f"{str(uuid.uuid4())}.tif", 1, 1, extent, 4326, data_type=gdal.GDT_Byte, fill_value=1) + return create_test_dataset( + f"{str(uuid.uuid4())}.tif", 1, 1, extent, 4326, data_type=gdal.GDT_Byte, fill_value=1 + ) # Create a mock for the WPSDataset class @@ -59,7 +70,10 @@ def create_mock_input_dataset_context(num: int): input_datasets = create_mock_wps_datasets(num) @contextmanager - def mock_input_dataset_context(_: List[str]): + def mock_input_dataset_context(dataset_keys: List[str]): + for ds, key in zip(input_datasets, dataset_keys, strict=True): + ds.ds_path = key + try: # Enter each dataset's context and yield the list of instances with ExitStack() as stack: diff --git a/backend/packages/wps-tools/src/wps_tools/interpolation_plotter.py b/backend/packages/wps-tools/src/wps_tools/interpolation_plotter.py index 941819daf..a022dfcac 100644 --- a/backend/packages/wps-tools/src/wps_tools/interpolation_plotter.py +++ b/backend/packages/wps-tools/src/wps_tools/interpolation_plotter.py @@ -1,18 +1,19 @@ import argparse import asyncio -import sys from datetime import datetime, timedelta, timezone from pathlib import Path from typing import List, Optional import numpy as np from aiohttp import ClientSession -from wps_sfms.interpolation.source import ( +from wps_sfms.interpolation.field import ( DEW_POINT_LAPSE_RATE, LAPSE_RATE, - StationActualSource, - StationDewPointSource, - StationTemperatureSource, + build_attribute_field, + build_wind_vector_field, + compute_adjusted_values, + compute_rh, + compute_sea_level_values, ) from wps_shared.geospatial.spatial_interpolation import idw_interpolation from wps_shared.schemas.sfms import SFMSDailyActual @@ -22,53 +23,113 @@ from wps_tools.sfms_scatterplot import create_sfms_scatterplot -def interpolate_temp(sfms_actuals: List[SFMSDailyActual]): - temp_source = StationTemperatureSource(sfms_actuals) - lats, lons, elevs, values = temp_source.get_station_arrays(only_valid=True) - sea = StationTemperatureSource.compute_sea_level_values(values, elevs, LAPSE_RATE) - assert len(lats) == len(lons) == len(elevs) == len(values) == len(sea) +MIN_WIND_SPEED_FOR_DIRECTION_KMH = 1.0 + + +def leave_one_out_idw(lats: np.ndarray, lons: np.ndarray, values: np.ndarray) -> np.ndarray: + """Leave-one-out IDW interpolation for station arrays.""" n = len(lats) - # Mask used for 'leave one out' analysis mask = np.ones(n, dtype=bool) - sea_interpolated = np.full(n, -sys.float_info.max) + interpolated_values = np.full(n, np.nan, dtype=np.float32) + if n < 2: + return interpolated_values for i in range(n): mask[i] = False - lat = lats[i] - lon = lons[i] - sea_interpolated[i] = idw_interpolation(lat, lon, lats[mask], lons[mask], sea[mask]) + interpolated_values[i] = idw_interpolation( + lats[i], lons[i], lats[mask], lons[mask], values[mask] + ) mask[i] = True - interpolated_values = StationTemperatureSource.compute_adjusted_values( - sea_interpolated, elevs, LAPSE_RATE - ) + return interpolated_values + + +def interpolate_temp(sfms_actuals: List[SFMSDailyActual]): + valid = [s for s in sfms_actuals if s.temperature is not None and s.elevation is not None] + lats = np.array([s.lat for s in valid], dtype=np.float32) + lons = np.array([s.lon for s in valid], dtype=np.float32) + elevs = np.array([s.elevation for s in valid], dtype=np.float32) + values = np.array([s.temperature for s in valid], dtype=np.float32) + sea = compute_sea_level_values(values, elevs, LAPSE_RATE) + assert len(lats) == len(lons) == len(elevs) == len(values) == len(sea) + sea_interpolated = leave_one_out_idw(lats, lons, sea) + + interpolated_values = compute_adjusted_values(sea_interpolated, elevs, LAPSE_RATE) return (elevs, values, interpolated_values) def interpolate_dewpoint_temp(sfms_actuals: List[SFMSDailyActual]): - dewpoint_source = StationDewPointSource(sfms_actuals) - lats, lons, elevs, values = dewpoint_source.get_station_arrays(only_valid=True) - sea = StationDewPointSource.compute_sea_level_values(values, elevs, DEW_POINT_LAPSE_RATE) + valid = [s for s in sfms_actuals if s.dewpoint is not None and s.elevation is not None] + lats = np.array([s.lat for s in valid], dtype=np.float32) + lons = np.array([s.lon for s in valid], dtype=np.float32) + elevs = np.array([s.elevation for s in valid], dtype=np.float32) + values = np.array([s.dewpoint for s in valid], dtype=np.float32) + sea = compute_sea_level_values(values, elevs, DEW_POINT_LAPSE_RATE) assert len(lats) == len(lons) == len(elevs) == len(values) == len(sea) - n = len(lats) - mask = np.ones(n, dtype=bool) - sea_interpolated = np.ones(n, dtype=np.float32) + sea_interpolated = leave_one_out_idw(lats, lons, sea) - for i in range(n): - mask[i] = False - lat = lats[i] - lon = lons[i] - sea_interpolated[i] = idw_interpolation(lat, lon, lats[mask], lons[mask], sea[mask]) - mask[i] = True - - interpolated_values = StationDewPointSource.compute_adjusted_values( + interpolated_values = compute_adjusted_values( sea_interpolated, elevs, DEW_POINT_LAPSE_RATE ) return (elevs, values, interpolated_values) -async def main(start: datetime, end: datetime, out_dir: Path): +def interpolate_wind_speed(sfms_actuals: List[SFMSDailyActual]): + valid = [s for s in sfms_actuals if s.wind_speed is not None] + lats = np.array([s.lat for s in valid], dtype=np.float32) + lons = np.array([s.lon for s in valid], dtype=np.float32) + observed_wind_speed = np.array([s.wind_speed for s in valid], dtype=np.float32) + interpolated_wind_speed = leave_one_out_idw(lats, lons, observed_wind_speed) + return (observed_wind_speed, interpolated_wind_speed) + + +def circular_difference_degrees(generated: np.ndarray, reference: np.ndarray) -> np.ndarray: + """Smallest signed angular difference in degrees in [-180, 180).""" + return (generated - reference + 180.0) % 360.0 - 180.0 + + +def interpolate_wind_direction(sfms_actuals: List[SFMSDailyActual]): + valid = [ + s + for s in sfms_actuals + if s.wind_speed is not None + and s.wind_direction is not None + and s.lat is not None + and s.lon is not None + ] + observed_wind_speed = np.array([s.wind_speed for s in valid], dtype=np.float32) + observed_wind_direction = np.array([s.wind_direction for s in valid], dtype=np.float32) + + wind_vector_field = build_wind_vector_field(valid) + lats, lons, station_u, station_v = ( + wind_vector_field.lats, + wind_vector_field.lons, + wind_vector_field.u, + wind_vector_field.v, + ) + interpolated_u = leave_one_out_idw(lats, lons, station_u) + interpolated_v = leave_one_out_idw(lats, lons, station_v) + + interpolated_wind_direction = np.zeros(len(interpolated_u), dtype=np.float32) + eps = np.float32(1e-6) + zero_v = np.abs(interpolated_v) < eps + nonzero_v = ~zero_v + interpolated_wind_direction[nonzero_v] = ( + np.degrees(np.arctan2(interpolated_u[nonzero_v], interpolated_v[nonzero_v])) + 180.0 + ) + zero_u = np.abs(interpolated_u) < eps + interpolated_wind_direction[zero_v & (interpolated_u < 0.0)] = 90.0 + interpolated_wind_direction[zero_v & (interpolated_u > 0.0)] = 270.0 + interpolated_wind_direction[zero_v & zero_u] = 0.0 + return ( + observed_wind_speed, + observed_wind_direction, + interpolated_wind_direction, + ) + + +async def main(start: datetime, end: datetime, out_dir: Optional[Path]): # Assemble the data needed for comparison async with ClientSession() as client_session: wfwx_api = WfwxApi(client_session) @@ -80,11 +141,12 @@ async def main(start: datetime, end: datetime, out_dir: Path): elevs, observed_temps, interpolated_temps = interpolate_temp(sfms_actuals) _, _, interpolated_dewpoint_temps = interpolate_dewpoint_temp(sfms_actuals) - interpolated_rh_from_observed_temps = StationDewPointSource.compute_rh(observed_temps, interpolated_dewpoint_temps) + interpolated_rh_from_observed_temps = compute_rh( + observed_temps, interpolated_dewpoint_temps + ) # Get observed rh values for comparison - rh_source = StationActualSource("relative_humidity", sfms_actuals) - _, _, observed_rh = rh_source.get_interpolation_data() + observed_rh = build_attribute_field(sfms_actuals, "relative_humidity").values observed_rh_np = np.array(observed_rh) # Weather stations show a rh of 0.0 when no observation is present observed_rh_masked = observed_rh_np[observed_rh_np > 0.0] @@ -133,6 +195,46 @@ async def main(start: datetime, end: datetime, out_dir: Path): f"{out_dir}/rh_diff_histogram_{datetime_to_process.date()}.png", ) + observed_wind_speed, interpolated_wind_speed = interpolate_wind_speed(sfms_actuals) + wind_speed_difference = observed_wind_speed - interpolated_wind_speed + + create_sfms_histogram( + wind_speed_difference, + "Difference: Observed - Interpolated", + "", + f"Wind Speed Difference: Observed - Interpolated - {datetime_to_process.date()}", + f"{out_dir}/wind_speed_diff_histogram_{datetime_to_process.date()}.png", + ) + + ( + wind_direction_observed_speed, + observed_wind_direction, + interpolated_wind_direction, + ) = interpolate_wind_direction(sfms_actuals) + wind_direction_difference = circular_difference_degrees( + observed_wind_direction, interpolated_wind_direction + ) + + wind_direction_mask = ( + np.isfinite(wind_direction_difference) + & np.isfinite(wind_direction_observed_speed) + & (wind_direction_observed_speed >= MIN_WIND_SPEED_FOR_DIRECTION_KMH) + ) + + masked_direction_error = wind_direction_difference[wind_direction_mask] + + create_sfms_histogram( + masked_direction_error, + "Difference: Observed - Interpolated (deg)", + "", + ( + "Wind Direction Difference: Observed - Interpolated " + f"(wind >= {MIN_WIND_SPEED_FOR_DIRECTION_KMH:.1f} km/h) - " + f"{datetime_to_process.date()}" + ), + f"{out_dir}/wind_direction_diff_histogram_{datetime_to_process.date()}.png", + ) + if __name__ == "__main__": parser = argparse.ArgumentParser( @@ -147,8 +249,7 @@ async def main(start: datetime, end: datetime, out_dir: Path): args = parser.parse_args() - out_dir: Optional[Path] = None start = datetime.strptime(args.start, "%Y-%m-%d") end = datetime.strptime(args.end, "%Y-%m-%d") - out_dir = args.out_dir + out_dir: Optional[Path] = Path(args.out_dir) if args.out_dir is not None else None asyncio.run(main(start, end, out_dir)) diff --git a/backend/packages/wps-tools/src/wps_tools/sfms_diff.py b/backend/packages/wps-tools/src/wps_tools/sfms_diff.py index 4f55f569e..b5c190219 100644 --- a/backend/packages/wps-tools/src/wps_tools/sfms_diff.py +++ b/backend/packages/wps-tools/src/wps_tools/sfms_diff.py @@ -1,5 +1,6 @@ import numpy as np import matplotlib +from pathlib import Path from wps_tools.sfms_raster_input_parser import parse_input, create_parser matplotlib.use("Agg") # Non-interactive backend @@ -10,6 +11,8 @@ def main(): data1, data2, valid1, valid2, diff, args = parse_input( create_parser(description="Compare two SFMS raster files and generate a difference image") ) + generated_name = Path(args.generated).name + reference_name = Path(args.reference).name # Compute value ranges from valid data valid_data1 = data1[valid1] @@ -26,19 +29,19 @@ def main(): # Generated ax1 = axes[0] im1 = ax1.imshow(np.where(valid1, data1, np.nan), cmap="coolwarm", vmin=data_min, vmax=data_max) - ax1.set_title("Generated") + ax1.set_title(f"Generated\n{generated_name}") plt.colorbar(im1, ax=ax1) # Reference ax2 = axes[1] im2 = ax2.imshow(np.where(valid2, data2, np.nan), cmap="coolwarm", vmin=data_min, vmax=data_max) - ax2.set_title("Reference") + ax2.set_title(f"Reference\n{reference_name}") plt.colorbar(im2, ax=ax2) # Difference ax3 = axes[2] im3 = ax3.imshow(diff, cmap="RdBu_r", vmin=-diff_abs_max, vmax=diff_abs_max) - ax3.set_title("Difference (Gen - Ref)") + ax3.set_title(f"Difference (Gen - Ref)\n{generated_name} - {reference_name}") plt.colorbar(im3, ax=ax3) plt.tight_layout() diff --git a/backend/packages/wps-tools/src/wps_tools/sfms_histogram.py b/backend/packages/wps-tools/src/wps_tools/sfms_histogram.py index 37981015b..d1e1db4bc 100644 --- a/backend/packages/wps-tools/src/wps_tools/sfms_histogram.py +++ b/backend/packages/wps-tools/src/wps_tools/sfms_histogram.py @@ -1,5 +1,6 @@ import matplotlib import numpy as np +from pathlib import Path from wps_tools.sfms_raster_input_parser import create_parser, parse_input @@ -44,9 +45,11 @@ def main(): _, _, _, _, diff, args = parse_input( create_parser(description="Generate histogram of differences between two SFMS raster files") ) + generated_name = Path(args.generated).name + reference_name = Path(args.reference).name x_label = "Difference (Generated - Reference)" y_label = "Pixel Count" - title = "Histogram of Raster Differences" + title = f"Histogram of Raster Differences\n{generated_name} - {reference_name}" create_sfms_histogram(diff, x_label, y_label, title, args.output) diff --git a/backend/uv.lock b/backend/uv.lock index 16ae9523c..70a61ff32 100644 --- a/backend/uv.lock +++ b/backend/uv.lock @@ -4659,6 +4659,8 @@ version = "0.1.0" source = { editable = "packages/wps-tools" } dependencies = [ { name = "aiobotocore" }, + { name = "aiohttp" }, + { name = "matplotlib" }, { name = "pandas" }, { name = "python-decouple" }, { name = "requests" }, @@ -4668,6 +4670,8 @@ dependencies = [ [package.metadata] requires-dist = [ { name = "aiobotocore", specifier = ">=2.4.1,<3" }, + { name = "aiohttp", specifier = ">=3.13.2" }, + { name = "matplotlib", specifier = ">=3.10.8" }, { name = "pandas", specifier = ">=1.5.2,<3" }, { name = "python-decouple", specifier = ">=3.6,<4" }, { name = "requests", specifier = ">=2.31.0,<3" },