Merge pull request #145 from Deltares/feat/DEI-263-monthly-average

Feat/DEI-263 multi year monthly average

Author: wschoonveld

decoimpact/business/entities/rules/time_aggregation_rule.py

@@ -11,11 +11,10 @@
     TimeAggregationRule
 """
 
-from typing import List
+from typing import List, Optional
 
 import numpy as _np
 import xarray as _xr
-from xarray.core.resample import DataArrayResample
 
 from decoimpact.business.entities.rules.i_array_based_rule import IArrayBasedRule
 from decoimpact.business.entities.rules.rule_base import RuleBase
@@ -36,18 +35,44 @@ def __init__(
         name: str,
         input_variable_names: List[str],
         operation_type: TimeOperationType,
+        multi_year_start: Optional[int] = None,
+        multi_year_end: Optional[int] = None,
     ):
+        # pylint: disable=too-many-arguments
+        # pylint: disable=too-many-positional-arguments
         super().__init__(name, input_variable_names)
         self._settings = TimeOperationSettings({"month": "ME", "year": "YE"})
         self._settings.percentile_value = 0
         self._settings.operation_type = operation_type
         self._settings.time_scale = "year"
+        self._multi_year_start = multi_year_start
+        self._multi_year_end = multi_year_end
 
     @property
     def settings(self):
         """Time operation settings"""
         return self._settings
 
+    @property
+    def multi_year_start(self) -> Optional[int]:
+        """Start year for the aggregation (inclusive)"""
+        return self._multi_year_start
+
+    @multi_year_start.setter
+    def multi_year_start(self, value: Optional[int]) -> None:
+        # optional validation here
+        self._multi_year_start = value
+
+    @property
+    def multi_year_end(self) -> Optional[int]:
+        """End year for the aggregation (inclusive)"""
+        return self._multi_year_end
+
+    @multi_year_end.setter
+    def multi_year_end(self, value: Optional[int]) -> None:
+        # optional validation here
+        self._multi_year_end = value
+
     def validate(self, logger: ILogger) -> bool:
         """Validates if the rule is valid
 
@@ -56,6 +81,21 @@ def validate(self, logger: ILogger) -> bool:
         """
         return self.settings.validate(self.name, logger)
 
+    def filter_years(
+        self, time_dim_name: str, value_array: _xr.DataArray
+    ) -> _xr.DataArray:
+        """Filters the values for the specified start and end year
+
+        Args:
+            value_array (DataArray): value to filter
+        """
+        start = (
+            str(self._multi_year_start) if self._multi_year_start is not None else None
+        )
+        end = str(self._multi_year_end) if self._multi_year_end is not None else None
+        slice_obj = slice(start, end)
+        return value_array.sel({time_dim_name: slice_obj})
+
     def execute(self, value_array: _xr.DataArray, logger: ILogger) -> _xr.DataArray:
         """Aggregates the values for the specified start and end date
 
@@ -81,14 +121,19 @@ def execute(self, value_array: _xr.DataArray, logger: ILogger) -> _xr.DataArray:
 
         dim_name = get_dict_element(settings.time_scale, settings.time_scale_mapping)
 
+        # get the name of the time dimension
         time_dim_name = get_time_dimension_name(value_array, logger)
-        aggregated_values = value_array.resample({time_dim_name: dim_name}, skipna=True)
-
-        result = self._perform_operation(aggregated_values)
-        # create a new aggregated time dimension based on original time dimension
 
-        result_time_dim_name = f"{time_dim_name}_{settings.time_scale}"
-        result = result.rename({time_dim_name: result_time_dim_name})
+        # perform aggregations in case of multi-year monthly average
+        if TimeOperationType.MULTI_YEAR_MONTHLY_AVERAGE == settings.operation_type:
+            result, result_time_dim_name = self._handle_grouping_operation(
+                time_dim_name, value_array
+            )
+        # perform the operation in case of other operation types
+        else:
+            result, result_time_dim_name = self._handle_aggregation_operation(
+                settings, dim_name, time_dim_name, value_array
+            )
 
         for key, value in value_array[time_dim_name].attrs.items():
             if value:
@@ -102,12 +147,73 @@ def execute(self, value_array: _xr.DataArray, logger: ILogger) -> _xr.DataArray:
 
         return result
 
-    def _perform_operation(self, aggregated_values: DataArrayResample) -> _xr.DataArray:
+    def _handle_grouping_operation(
+        self,
+        time_dim_name: str,
+        value_array: _xr.DataArray,
+    ) -> tuple[_xr.DataArray, str]:
+        """Handles the grouping operation for multi-year monthly average"""
+        filtered_values = self.filter_years(time_dim_name, value_array)
+        grouped_values = filtered_values.groupby(f"{time_dim_name}.month")
+        grouping_result = self._perform_grouping_operation(grouped_values)
+        # create a new aggregated time dimension based on original time dimension
+        _result_time_dim_name = f"{time_dim_name}_monthly"
+        grouping_result = grouping_result.rename({"month": _result_time_dim_name})
+        return grouping_result, _result_time_dim_name
+
+    def _handle_aggregation_operation(
+        self,
+        settings: TimeOperationSettings,
+        dim_name: str,
+        time_dim_name: str,
+        value_array: _xr.DataArray,
+    ) -> tuple[_xr.DataArray, str]:
+        """Handles the aggregation operation for other operation types"""
+        aggregated_values = value_array.resample({time_dim_name: dim_name}, skipna=True)
+        # create a new aggregated time dimension based on original time dimension
+        aggregation_result = self._perform_operation(
+            aggregated_values, settings.operation_type
+        )
+        _result_time_dim_name = f"{time_dim_name}_{settings.time_scale}"
+        aggregation_result = aggregation_result.rename(
+            {time_dim_name: _result_time_dim_name}
+        )
+        return aggregation_result, _result_time_dim_name
+
+    def _perform_grouping_operation(
+        self,
+        grouped_values,
+    ) -> _xr.DataArray:
+        """Returns the values based on the grouping operation type
+
+        Args:
+            aggregated_values (DataArrayGroupBy): aggregate values
+            operation_type (TimeOperationType): the operation type
+            time_dim_name (str): name of the time dimension (for MONTHLY_AVERAGE)
+        Raises:
+            NotImplementedError: If operation type is not supported
+
+        Returns:
+            DataArray: Values of operation type
+        """
+        # Compute mean across years for each calendar month
+        monthly = grouped_values.mean(skipna=True)
+        # Ensure all 12 months are present (1..12), insert NaNs using reindex
+        months = _np.arange(1, 13)
+        result = monthly.reindex({"month": months})
+
+        return _xr.DataArray(result)
+
+    def _perform_operation(
+        self,
+        aggregated_values,
+        operation_type: TimeOperationType,
+    ) -> _xr.DataArray:
         """Returns the values based on the operation type
 
         Args:
             aggregated_values (DataArrayResample): aggregate values
-
+            operation_type (TimeOperationType): the operation type
         Raises:
             NotImplementedError: If operation type is not supported
 
@@ -120,8 +226,6 @@ def _perform_operation(self, aggregated_values: DataArrayResample) -> _xr.DataAr
             TimeOperationType.AVG_DURATION_PERIODS,
         ]
 
-        operation_type = self.settings.operation_type
-
         if operation_type is TimeOperationType.ADD:
             result = aggregated_values.sum()
 

decoimpact/business/workflow/model_builder.py

@@ -147,6 +147,8 @@ def _create_rule(rule_data: IRuleData) -> IRule:
             )
             rule.settings.percentile_value = rule_data.percentile_value
             rule.settings.time_scale = rule_data.time_scale
+            rule.multi_year_start = rule_data.multi_year_start
+            rule.multi_year_end = rule_data.multi_year_end
         elif isinstance(rule_data, IRollingStatisticsRuleData):
             rule = RollingStatisticsRule(
                 rule_data.name, [rule_data.input_variable], rule_data.operation
@@ -159,7 +161,7 @@ def _create_rule(rule_data: IRuleData) -> IRule:
                 rule_data.name,
                 rule_data.input_variable_names,
                 MultiArrayOperationType[rule_data.operation_type],
-                rule_data.ignore_nan
+                rule_data.ignore_nan,
             )
         elif isinstance(rule_data, IResponseCurveRuleData):
             rule = ResponseCurveRule(

decoimpact/data/api/i_time_aggregation_rule_data.py

@@ -14,6 +14,7 @@
 
 
 from abc import ABC, abstractmethod
+from typing import Optional
 
 from decoimpact.data.api.i_rule_data import IRuleData
 from decoimpact.data.api.time_operation_type import TimeOperationType
@@ -41,3 +42,13 @@ def percentile_value(self) -> float:
     @abstractmethod
     def time_scale(self) -> str:
         """Time scale"""
+
+    @property
+    @abstractmethod
+    def multi_year_start(self) -> Optional[int]:
+        """Start year for the aggregation (inclusive)"""
+
+    @property
+    @abstractmethod
+    def multi_year_end(self) -> Optional[int]:
+        """End year for the aggregation (inclusive)"""

decoimpact/data/api/time_operation_type.py

@@ -26,3 +26,4 @@ class TimeOperationType(IntEnum):
     AVG_DURATION_PERIODS = 8
     STDEV = 9
     PERCENTILE = 10
+    MULTI_YEAR_MONTHLY_AVERAGE = 11

decoimpact/data/entities/time_aggregation_rule_data.py

@@ -12,6 +12,8 @@
 
 """
 
+from typing import Optional
+
 from decoimpact.data.api.i_time_aggregation_rule_data import ITimeAggregationRuleData
 from decoimpact.data.api.time_operation_type import TimeOperationType
 from decoimpact.data.entities.time_operation_rule_data import TimeOperationRuleData
@@ -20,11 +22,42 @@
 class TimeAggregationRuleData(TimeOperationRuleData, ITimeAggregationRuleData):
     """Class for storing data related to time_aggregation rule"""
 
-    def __init__(self, name: str, operation: TimeOperationType, input_variable: str):
+    def __init__(
+        self,
+        name: str,
+        operation: TimeOperationType,
+        input_variable: str,
+        multi_year_start: Optional[int] = None,
+        multi_year_end: Optional[int] = None,
+    ):
+        # pylint: disable=too-many-arguments
+        # pylint: disable=too-many-positional-arguments
         super().__init__(name, operation)
         self._input_variable = input_variable
+        self._multi_year_start = multi_year_start
+        self._multi_year_end = multi_year_end
 
     @property
     def input_variable(self) -> str:
         """Name of the input variable"""
         return self._input_variable
+
+    @property
+    def multi_year_start(self) -> Optional[int]:
+        """Start year for aggregation"""
+        return self._multi_year_start
+
+    @multi_year_start.setter
+    def multi_year_start(self, value: Optional[int]) -> None:
+        # optional validation here
+        self._multi_year_start = value
+
+    @property
+    def multi_year_end(self) -> Optional[int]:
+        """End year for aggregation"""
+        return self._multi_year_end
+
+    @multi_year_end.setter
+    def multi_year_end(self, value: Optional[int]) -> None:
+        # optional validation here
+        self._multi_year_end = value

decoimpact/data/parsers/parser_time_aggregation_rule.py

@@ -42,6 +42,8 @@ def parse_dict(self, dictionary: Dict[str, Any], logger: ILogger) -> IRuleData:
         input_variable_name: str = get_dict_element("input_variable", dictionary)
         operation: str = get_dict_element("operation", dictionary)
         time_scale: str = get_dict_element("time_scale", dictionary)
+        multi_year_start: int = get_dict_element("multi_year_start", dictionary, False)
+        multi_year_end: int = get_dict_element("multi_year_end", dictionary, False)
         output_variable_name: str = get_dict_element("output_variable", dictionary)
 
         operation_value, percentile_value = parse_operation_values(operation)
@@ -50,6 +52,8 @@ def parse_dict(self, dictionary: Dict[str, Any], logger: ILogger) -> IRuleData:
 
         rule_data.percentile_value = percentile_value
         rule_data.time_scale = time_scale
+        rule_data.multi_year_start = multi_year_start
+        rule_data.multi_year_end = multi_year_end
         rule_data.output_variable = output_variable_name
         rule_data.description = description
 

docs/manual/rules/time_aggregation_rule.md

@@ -9,12 +9,14 @@ FORMAT
       time_scale : <time_aggregation_applied>
       input_variable: <one_input_variable_name>
       output_variable: <one_output_variable_name>
+      multi_year_start: <start_year>
+      multi_year_end: <end_year>
 ```
 
 The time aggregation rule allows for calculating a statistical summary over the time axes of 3D and 2D variables. This could be used for calculating the maximum value over a year (e.g., for water level) or the minimum value over a month (e.g., oxygen concentration). The rule operates both on 3D variables and 2D variables as long as they have a time axis and returns a 3D or 2D result depending on input with the statistic calculated for a new time axis (e.g., year or month).
-Operations available: Add, Average, Median, Min, Max, period statistics, Stdev and Percentile(n). When using percentile, add a number for the nth percentile with brackets like this: percentile(10). Stdev calculates the standard- deviation over the time period. Under period statistics are explained further in the text.
+Operations available: Add, Average, Multi_year_monthly_average, Median, Min, Max, period statistics, Stdev and Percentile(n). When using percentile, add a number for the nth percentile with brackets like this: percentile(10). Stdev calculates the standard- deviation over the time period. Under period statistics are explained further in the text.
 
-Time aggregation available: Year, Month
+Time aggregation available: Year, Month, Multi_year_monthly_average
 
 The rule needs to be applied to an existing 2D/3D variable with time axis. A new 2D/3D variable with new time axis is created when the rule is executed. With a year timestep the result is written to the last day of the year, with a month timestep the result is written to the last day of the month per year.
 
@@ -29,6 +31,10 @@ The rule needs to be applied to an existing 2D/3D variable with time axis. A new
       output_variable: MAX_water_level_year
 ```
 
+To aggregate an average per month, use the combination of operation AVERAGE with time_scale MONTH.
+To aggregate a multi-year average over all months, use the combination of operation MULTI_YEAR_MONTHLY_AVERAGE with time_scale MONTH. Then, for each month, the average over the time frame is calculated. The output variable is related to an extra time dimension time_monthly with all 12 months (in numbers).
+In combination with the operation MULTI_YEAR_MONTHLY_AVERAGE it is possible do this is for a specific range of years by adding the optional parameters MULTI_YEAR_START and MULTI_YEAR_END.
+
 ![Result Time aggregation rule](../../assets/images/3_result_time_aggregation.png "Water level (in m NAP, left-hand side) with a timestep every 10 days has been summarized to the maximum for each year (right-hand side) while maintaining the face dimension (layer dimension is not present in this example, but would be maintained).")
 
 Period statistics: Time aggregation rule with COUNT_PERIODS, AVG_DURATION_PERIODS, MIN_DURATION_PERIODS and MAX_DURATION_PERIODS

template_input.yaml

@@ -29,6 +29,8 @@ rules:
       time_scale:
       input_variable:
       output_variable:
+      multi_year_start:
+      multi_year_end:
 
     - step_function_rule:
       name: