New multirun plotter

validation/multirun/profiles_plotter/config.yaml

@@ -0,0 +1,35 @@
+sources:
+  enea_raw:
+    path: /media/internal/disk2TiB/data/FORCINGS_ENEA_RAW_CONVERTED/AVESCAN_OUTPUT/STAT_PROFILES_COMPACT
+    meshmask: /media/internal/disk2TiB/data/FORCINGS_ENEA_RAW_CONVERTED/meshmask_MITgcm.nc
+  enea_interpolated:
+    path: /somewhere/else
+    meshmask: /media/internal/disk2TiB/data/FORCINGS_ENEA_RAW_CONVERTED/meshmask_MITgcm.nc
+
+variable_selections:
+  v1: [P_l, N3n, N1p, N4n, N5s, P_c, O2o]
+  v2: [N1p, N3n, N4n, O2o, P_l, P_c, DIC, ppn, ALK, DIC, pH, pCO2, O3c, O3h]
+  var_enea: [THETA, SALT]
+
+levels: ["0", "0 - 20", "100", "150"] # in meters
+
+plots:
+  enea_original:
+    source: enea_raw
+    variables: var_enea
+    color: red
+    alpha_for_time_series: 0.2
+    zorder: 1
+    legend: enea
+
+  enea_trend:
+    source: enea_raw
+    variables: var_enea
+    active: true
+    color: blue
+    alpha: 1
+    zorder: 1
+    filter: MovingAverage(10y)
+    draw_depth_profile: false
+
+output_dir: /dev/shm

validation/multirun/profiles_plotter/filters/filter.py

@@ -0,0 +1,24 @@
+from abc import ABC, abstractmethod
+from typing import Union
+
+from tools.data_object import DataObject
+
+
+class InvalidFilterDescription(Exception):
+    pass
+
+
+class FilteredObject(DataObject, ABC):
+    def __init__(self, data_object: DataObject):
+        self._original_data = data_object
+
+
+class Filter(ABC):
+    @abstractmethod
+    def get_filtered_object(self, data_object) -> FilteredObject:
+        raise NotImplementedError
+
+    @staticmethod
+    @abstractmethod
+    def initialize_from_string(args_str: Union[str, None]):
+        raise NotImplementedError

validation/multirun/profiles_plotter/filters/moving_average_filter.py

@@ -0,0 +1,186 @@
+from datetime import timedelta
+import re
+
+import numpy as np
+
+from filters.filter import Filter, FilteredObject
+from tools.data_object import DataObject
+
+
+class MovingAverageFilteredObject(FilteredObject):
+    def __init__(self, data_object: DataObject, window_size: int):
+        super().__init__(data_object)
+
+        self._window_size = int(window_size)
+
+        # How many points to we lose because of the average (the ones on the
+        # boundary of the time series)
+        lost_points = self._window_size - 1
+
+        original_times = self._original_data.get_time_steps()
+
+        if self._window_size % 2 == 1:
+            lost_per_side = lost_points // 2
+            self._time_steps = tuple(
+                original_times[lost_per_side: - lost_per_side]
+            )
+        else:
+            lost_per_side = (lost_points - 1) // 2
+            new_time_steps = []
+            averaged_times = original_times[lost_per_side:- lost_per_side]
+            for i, t1, in enumerate(averaged_times[:-1]):
+                t2 = averaged_times[i + 1]
+                delta_time = t2 - t1
+                new_time_steps.append(t1 + delta_time // 2)
+            self._time_steps = tuple(new_time_steps)
+            assert len(self._time_steps) == len(original_times) - lost_points
+
+    def get_time_steps(self):
+        return self._time_steps
+
+    def get_values(self, time_steps, basin, level_index, indicator, coasts=1):
+        fixed_axis = set()
+        if not isinstance(basin, slice):
+            fixed_axis.add('basin')
+        if not isinstance(level_index, slice):
+            fixed_axis.add('level')
+        if not isinstance(coasts, slice):
+            fixed_axis.add('coasts')
+        if not isinstance(indicator, slice):
+            fixed_axis.add('indicator')
+
+        time_axis = self.get_axis('time', fixed_axis)
+
+        if isinstance(time_steps, slice):
+            if time_steps.stop is None:
+                new_slice_stop = None
+            else:
+                if time_steps.stop < 0:
+                    new_slice_stop = time_steps.stop
+                else:
+                    new_slice_stop = time_steps.stop + self._window_size - 1
+            time_slice = slice(time_steps.start, new_slice_stop)
+        else:
+            # Here we can suppose that time_steps is an integer
+            if time_steps >= 0:
+                time_slice_start = time_steps
+                time_slice_stop = time_steps + self._window_size
+            else:
+                time_slice_start = time_steps - self._window_size
+                time_slice_stop = time_steps
+            time_slice = slice(time_slice_start, time_slice_stop)
+
+        with self._original_data:
+            original_data = self._original_data.get_values(
+                time_steps=time_slice,
+                basin=basin,
+                level_index=level_index,
+                indicator=indicator,
+                shores=coasts
+            )
+
+        if isinstance(time_steps, int):
+            if original_data.shape[time_axis] != self._window_size:
+                raise IndexError(
+                    'Index {} outside valid range'.format(time_steps)
+                )
+
+            return np.mean(original_data, axis=time_axis)
+
+        all_slice_none = [
+            slice(None) for _ in range(len(original_data.shape))
+        ]
+
+        time_index_geq_0 = all_slice_none.copy()
+        time_index_geq_0[time_axis] = slice(1, None)
+        time_index_geq_0 = tuple(time_index_geq_0)
+
+        cut_first_time_indices = all_slice_none.copy()
+        cut_first_time_indices[time_axis] = slice(self._window_size - 1, None)
+        cut_first_time_indices = tuple(cut_first_time_indices)
+
+        cut_last_time_indices = all_slice_none.copy()
+        cut_last_time_indices[time_axis] = slice(None, - self._window_size)
+        cut_last_time_indices = tuple(cut_last_time_indices)
+
+        v1 = np.cumsum(
+            original_data,
+            axis=time_axis
+        )[cut_first_time_indices]
+
+        v2 = np.zeros_like(v1)
+        v2[time_index_geq_0] = np.cumsum(
+            original_data,
+            axis=time_axis
+        )[cut_last_time_indices]
+
+        moving_average = (v1 - v2) / self._window_size
+
+        return moving_average
+
+
+class MovingAverageFilter(Filter):
+    WINDOW_DESCRIPTION_MASK = re.compile(
+        r'^(?P<size>\d+)\s{0,2}(?P<unit>(d|D|m|M|y|Y)?)\s*$'
+    )
+
+    def __init__(self, window_description: str):
+        mask_match = self.WINDOW_DESCRIPTION_MASK.match(window_description)
+        if mask_match is None:
+            raise ValueError(
+                'The string that describes the window size of the moving '
+                'average filter must be an integer, optionally followed by one '
+                'of the following letters: d, m or y'
+            )
+        self._window_size = int(mask_match.group('size'))
+        if mask_match.group('unit') is not None:
+            self._unit = mask_match.group('unit').lower()
+        else:
+            self._unit = None
+
+    def get_filtered_object(self, data_object):
+        # If there are no units, we can simply return the filtered object as is
+        if self._unit is None:
+            return MovingAverageFilteredObject(
+                data_object,
+                window_size=self._windows_size
+            )
+
+        # Otherwise, we need to estimate the new size of the window
+        if self._unit == 'd':
+            unit_time = timedelta(days=1)
+        elif self._unit == 'm':
+            unit_time = timedelta(days=30)
+        else:
+            unit_time = timedelta(days=365)
+
+        time_steps = data_object.get_time_steps()
+        time_index = 0
+        for time_index, current_time in enumerate(time_steps):
+            if current_time - time_steps[0] >= unit_time:
+                break
+
+        if time_index == 0:
+            raise ValueError('No time steps found in the data_object')
+
+        window_size = self._window_size * time_index
+        return MovingAverageFilteredObject(data_object, window_size=window_size)
+
+    @staticmethod
+    def initialize_from_string(args_str):
+        if args_str is None:
+            raise ValueError(
+                'A moving average filter requires a mandatory argument that '
+                'must be the size of the average window'
+            )
+
+        args_split = args_str.split(',')
+        if len(args_split) > 1:
+            raise ValueError(
+                'A moving average filter requires only one argument. '
+                'Received: {}'.format(
+                    ', '.join([k.strip() for k in args_split])
+                )
+            )
+        window_desc = args_split[0].strip()
+        return MovingAverageFilter(window_desc)

validation/multirun/profiles_plotter/filters/read_filter_description.py

@@ -0,0 +1,45 @@
+import re
+from filters.filter import InvalidFilterDescription
+from filters.moving_average_filter import MovingAverageFilter
+
+
+FILTER_DESCRIPTION_MASK = re.compile(
+    r'^\s*(?P<filter_name>[a-zA-Z][a-zA-Z0-9_]*)((?P<args>\(.*\))?)\s*$'
+)
+
+
+# Add here the name of the filter with its association every time you introduce
+# a new filter
+_FILTER_ASSOCIATION = {
+    'movingaverage': MovingAverageFilter
+}
+
+
+def read_filter_description(filter_description: str):
+    mask_match = FILTER_DESCRIPTION_MASK.match(filter_description)
+    if mask_match is None:
+        raise InvalidFilterDescription(
+            'Invalid filter description! A filter must satisfy the following '
+            'regular expression: {}\nReceived: "{}"'.format(
+                FILTER_DESCRIPTION_MASK,
+                filter_description
+            )
+        )
+    filter_name = mask_match.group('filter_name')
+    filter_args = mask_match.group('args')
+
+    filter_name_lower = filter_name.lower()
+    if filter_name_lower in _FILTER_ASSOCIATION:
+        # Remove the parenthesis
+        if filter_args is not None:
+            filter_args = filter_args[1:-1]
+        return _FILTER_ASSOCIATION[filter_name_lower].initialize_from_string(
+            filter_args
+        )
+
+    raise InvalidFilterDescription(
+        'Invalid filter specified: "{}"\nAllowed values are: {}'.format(
+            filter_name,
+            ', '.join(sorted(list(_FILTER_ASSOCIATION.keys())))
+        )
+    )