Average time

HISTORY.rst

@@ -7,6 +7,7 @@ v0.8.0 (unreleased)
 New Features
 ^^^^^^^^^^^^
 * ``clisops.core.average.average_shape`` copies the global and variable attributes from the input data to the results.
+* ``clisops.ops.average.average_time`` and ``clisops.core.average.average_time`` added. Allowing averaging over time frequencies of day, month and year.
 
 Bug fixes
 ^^^^^^^^^

clisops/core/average.py

@@ -1,15 +1,24 @@
 """Average module."""
-from pathlib import Path
-from typing import Tuple, Union, Sequence
 import warnings
-import numpy as np
+from pathlib import Path
+from typing import Sequence, Tuple, Union
+
 import geopandas as gpd
+import numpy as np
 import xarray as xr
 from roocs_utils.exceptions import InvalidParameterValue
-from roocs_utils.xarray_utils.xarray_utils import get_coord_type, known_coord_types
+from roocs_utils.xarray_utils.xarray_utils import (
+    get_coord_by_type,
+    get_coord_type,
+    known_coord_types,
+)
+
 from .subset import shape_bbox_indexer
 
-__all__ = ["average_over_dims", "average_shape"]
+__all__ = ["average_over_dims", "average_shape", "average_time"]
+
+# see https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#dateoffset-objects
+freqs = {"day": "1D", "month": "1MS", "year": "1AS"}
 
 
 def average_shape(
@@ -216,3 +225,55 @@ def average_over_dims(
     if isinstance(ds, xr.Dataset):
         return xr.merge((ds_averaged_over_dims, untouched_ds))
     return ds_averaged_over_dims
+
+
+def average_time(
+    ds: Union[xr.DataArray, xr.Dataset],
+    freq: str,
+) -> Union[xr.DataArray, xr.Dataset]:
+    """
+    Average a DataArray or Dataset over the time frequency specified.
+
+    Parameters
+    ----------
+    ds : Union[xr.DataArray, xr.Dataset]
+      Input values.
+    freq: str
+      The frequency to average over. One of "month", "year".
+
+    Returns
+    -------
+    Union[xr.DataArray, xr.Dataset]
+      New Dataset or DataArray object averaged over the indicated time frequency.
+
+    Examples
+    --------
+    >>> from clisops.core.average import average_time  # doctest: +SKIP
+    >>> pr = xr.open_dataset(path_to_pr_file).pr  # doctest: +SKIP
+    ...
+    # Average data array over each month
+    >>> prAvg = average_time(pr, freq='month')  # doctest: +SKIP
+    """
+
+    if not freq:
+        raise InvalidParameterValue(
+            "At least one frequency for averaging must be provided"
+        )
+
+    if freq not in list(freqs.keys()):
+        raise InvalidParameterValue(
+            f"Time frequency for averaging must be one of {list(freqs.keys())}."
+        )
+
+    # check time coordinate exists and get name
+    t = get_coord_by_type(ds, "time", ignore_aux_coords=False)
+
+    if t is None:
+        raise Exception("Time dimension could not be found")
+
+    # resample and average over time
+    ds_avg_over_time = ds.resample(indexer={t.name: freqs[freq]}).mean(
+        dim=t.name, skipna=True, keep_attrs=True
+    )
+
+    return ds_avg_over_time

clisops/ops/average.py

@@ -2,6 +2,7 @@
 from typing import List, Optional, Tuple, Union
 
 import xarray as xr
+from roocs_utils.exceptions import InvalidParameterValue
 from roocs_utils.parameter.dimension_parameter import DimensionParameter
 from roocs_utils.xarray_utils.xarray_utils import (
     convert_coord_to_axis,
@@ -17,9 +18,7 @@
 from clisops.utils.file_namers import get_file_namer
 from clisops.utils.output_utils import get_output, get_time_slices
 
-__all__ = [
-    "average_over_dims",
-]
+__all__ = ["average_over_dims", "average_time"]
 
 LOGGER = logging.getLogger(__file__)
 
@@ -97,3 +96,75 @@ def average_over_dims(
     """
     op = Average(**locals())
     return op.process()
+
+
+class AverageTime(Operation):
+    def _resolve_params(self, **params):
+        freq = params.get("freq", None)
+
+        if not freq:
+            raise InvalidParameterValue(
+                "At least one frequency for averaging must be provided"
+            )
+
+        if freq not in list(average.freqs.keys()):
+            raise InvalidParameterValue(
+                f"Time frequency for averaging must be one of {list(average.freqs.keys())}."
+            )
+
+        self.params = {"freq": freq}
+
+    def _get_file_namer(self):
+        extra = f"_avg-{self.params.get('freq')}"
+        namer = get_file_namer(self._file_namer)(extra=extra)
+
+        return namer
+
+    def _calculate(self):
+        avg_ds = average.average_time(
+            self.ds,
+            self.params.get("freq", None),
+        )
+
+        return avg_ds
+
+
+def average_time(
+    ds,
+    freq: str,
+    output_dir: Optional[Union[str, Path]] = None,
+    output_type="netcdf",
+    split_method="time:auto",
+    file_namer="standard",
+) -> List[Union[xr.Dataset, str]]:
+    """
+
+    Parameters
+    ----------
+    ds: Union[xr.Dataset, str]
+    freq: str
+      The frequency to average over. One of "month", "year".
+    output_dir: Optional[Union[str, Path]] = None
+    output_type: {"netcdf", "nc", "zarr", "xarray"}
+    split_method: {"time:auto"}
+    file_namer: {"standard", "simple"}
+
+    Returns
+    -------
+    List[Union[xr.Dataset, str]]
+    A list of the outputs in the format selected; str corresponds to file paths if the
+    output format selected is a file.
+
+    Examples
+    --------
+    | ds: xarray Dataset or "cmip5.output1.MOHC.HadGEM2-ES.rcp85.mon.atmos.Amon.r1i1p1.latest.tas"
+    | dims: ['latitude', 'longitude']
+    | ignore_undetected_dims: False
+    | output_dir: "/cache/wps/procs/req0111"
+    | output_type: "netcdf"
+    | split_method: "time:auto"
+    | file_namer: "standard"
+
+    """
+    op = AverageTime(**locals())
+    return op.process()

docs/api.rst

@@ -29,8 +29,8 @@ Subset operation
    :show-inheritance:
 
 
-Average operation
-=================
+Average operations
+==================
 
 .. automodule:: clisops.ops.average
    :noindex:

tests/core/test_average.py

@@ -6,10 +6,12 @@
 import xarray as xr
 from pkg_resources import parse_version
 from roocs_utils.exceptions import InvalidParameterValue
+from roocs_utils.xarray_utils import xarray_utils as xu
 
 from clisops.core import average
 from clisops.utils import get_file
 
+from .._common import CMIP5_RH, CMIP6_SICONC_DAY
 from .._common import XCLIM_TESTS_DATA as TESTS_DATA
 
 try:
@@ -173,3 +175,54 @@ def test_average_wrong_format(self):
             str(exc.value)
             == "Dimensions for averaging must be one of ['time', 'level', 'latitude', 'longitude']"
         )
+
+
+class TestAverageTime:
+    month_ds = CMIP5_RH
+    day_ds = CMIP6_SICONC_DAY
+
+    def test_average_month(self):
+        ds = xu.open_xr_dataset(self.day_ds)
+        assert ds.time.shape == (60225,)
+
+        avg_ds = average.average_time(ds, freq="month")
+        assert avg_ds.time.shape == (1980,)
+
+    def test_average_year(self):
+        ds = xu.open_xr_dataset(self.month_ds)
+        assert ds.time.shape == (1752,)
+
+        avg_ds = average.average_time(ds, freq="year")
+        assert avg_ds.time.shape == (147,)
+
+    def test_no_freq(self):
+        ds = xu.open_xr_dataset(self.month_ds)
+
+        with pytest.raises(InvalidParameterValue) as exc:
+            average.average_time(ds, freq=None)
+        assert str(exc.value) == "At least one frequency for averaging must be provided"
+
+    def test_incorrect_freq(self):
+        ds = xu.open_xr_dataset(self.month_ds)
+
+        with pytest.raises(InvalidParameterValue) as exc:
+            average.average_time(ds, freq="wrong")
+        assert (
+            str(exc.value)
+            == "Time frequency for averaging must be one of ['day', 'month', 'year']."
+        )
+
+    def test_freq_wrong_format(self):
+        ds = xu.open_xr_dataset(self.month_ds)
+
+        with pytest.raises(InvalidParameterValue) as exc:
+            average.average_time(ds, freq=0)
+        assert str(exc.value) == "At least one frequency for averaging must be provided"
+
+    def test_no_time(self):
+        ds = xu.open_xr_dataset(self.month_ds)
+        ds = ds.drop_dims("time")
+
+        with pytest.raises(Exception) as exc:
+            average.average_time(ds, freq="year")
+        assert str(exc.value) == "Time dimension could not be found"

tests/ops/test_average.py

@@ -8,17 +8,17 @@
 
 import clisops
 from clisops import CONFIG
-from clisops.ops.average import average_over_dims
+from clisops.ops.average import average_over_dims, average_time
 
-from .._common import CMIP5_TAS
+from .._common import C3S_CORDEX_EUR_ZG500, CMIP5_TAS, CMIP6_SICONC_DAY
 
 
 def _check_output_nc(result, fname="output_001.nc"):
     assert fname in [os.path.basename(_) for _ in result]
 
 
 def _load_ds(fpath):
-    return xr.open_mfdataset(fpath)
+    return xr.open_mfdataset(fpath, use_cftime=True)
 
 
 def test_average_basic_data_array(cmip5_tas_file):
@@ -217,3 +217,110 @@ def test_aux_variables():
     )
 
     assert "do_i_get_written" in result[0].variables
+
+
+def test_average_over_years():
+    ds = _load_ds(CMIP5_TAS)  # monthly dataset
+
+    # check initial dataset
+    assert ds.time.shape == (3530,)
+    assert ds.time.values[0].isoformat() == "2005-12-16T00:00:00"
+    assert ds.time.values[-1].isoformat() == "2299-12-16T00:00:00"
+
+    result = average_time(
+        CMIP5_TAS,
+        freq="year",
+        output_type="xarray",
+    )
+
+    time_length = ds.time.values[-1].year - ds.time.values[0].year + 1
+
+    assert result[0].time.shape == (time_length,)  # get number of years
+    assert result[0].time.values[0].isoformat() == "2005-01-01T00:00:00"
+    assert result[0].time.values[-1].isoformat() == "2299-01-01T00:00:00"
+
+
+def test_average_over_months():
+    ds = _load_ds(CMIP6_SICONC_DAY)  # monthly dataset
+
+    # check initial dataset
+    assert ds.time.shape == (60225,)
+    assert ds.time.values[0].isoformat() == "1850-01-01T12:00:00"
+    assert ds.time.values[-1].isoformat() == "2014-12-31T12:00:00"
+
+    # average over time
+    result = average_time(
+        CMIP6_SICONC_DAY,
+        freq="month",
+        output_type="xarray",
+    )
+
+    time_length = (
+        ds.time.values[-1].year - ds.time.values[0].year + 1
+    ) * 12  # get number of months
+
+    assert result[0].time.shape == (time_length,)
+    assert result[0].time.values[0].isoformat() == "1850-01-01T00:00:00"
+    assert result[0].time.values[-1].isoformat() == "2014-12-01T00:00:00"
+
+
+def test_average_time_no_freq():
+    with pytest.raises(InvalidParameterValue) as exc:
+        # average over time
+        average_time(
+            CMIP6_SICONC_DAY,
+            freq=None,
+            output_type="xarray",
+        )
+    assert str(exc.value) == "At least one frequency for averaging must be provided"
+
+
+def test_average_time_incorrect_freq():
+    with pytest.raises(InvalidParameterValue) as exc:
+        # average over time
+        average_time(
+            CMIP6_SICONC_DAY,
+            freq="week",
+            output_type="xarray",
+        )
+    assert (
+        str(exc.value)
+        == "Time frequency for averaging must be one of ['day', 'month', 'year']."
+    )
+
+
+def test_average_time_file_name(tmpdir):
+    result = average_time(
+        CMIP5_TAS,
+        freq="year",
+        output_type="nc",
+        output_dir=tmpdir,
+    )
+
+    _check_output_nc(
+        result, fname="tas_mon_HadGEM2-ES_rcp85_r1i1p1_20050101-22990101_avg-year.nc"
+    )
+
+
+def test_average_time_cordex():
+    ds = _load_ds(C3S_CORDEX_EUR_ZG500)
+
+    # check initial dataset
+    assert ds.time.shape == (3653,)
+    assert ds.time.values[0].isoformat() == "2071-01-01T12:00:00"
+    assert ds.time.values[-1].isoformat() == "2080-12-31T12:00:00"
+
+    # average over time
+    result = average_time(
+        C3S_CORDEX_EUR_ZG500,
+        freq="month",
+        output_type="xarray",
+    )
+
+    time_length = (
+        ds.time.values[-1].year - ds.time.values[0].year + 1
+    ) * 12  # get number of months
+
+    assert result[0].time.shape == (time_length,)
+    assert result[0].time.values[0].isoformat() == "2071-01-01T00:00:00"
+    assert result[0].time.values[-1].isoformat() == "2080-12-01T00:00:00"