From 57c26220d62f5825f836bb177c80053c029eec56 Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Mon, 2 Oct 2023 12:35:01 -0400
Subject: [PATCH 1/8] Resampling from coarser than daily

---
 xscen/extract.py | 61 +++++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 52 insertions(+), 9 deletions(-)

diff --git a/xscen/extract.py b/xscen/extract.py
index 7575a526..222475ba 100644
--- a/xscen/extract.py
+++ b/xscen/extract.py
@@ -397,15 +397,30 @@ def resample(
     initial_frequency_td = pd.Timedelta(
         CV.xrfreq_to_timedelta(xr.infer_freq(da.time.dt.round("T")), None)
     )
-    if initial_frequency_td == pd.Timedelta("1D"):
-        logger.warning(
-            "You appear to be resampling daily data using extract_dataset. "
-            "It is advised to use compute_indicators instead, as it is far more robust."
+    days_per_step = None
+    if initial_frequency != "undetected" and initial_frequency_td > pd.Timedelta(
+        7, "days"
+    ):
+        # More than a week -> non-uniform sampling length!
+        t = xr.date_range(
+            da.indexes["time"][0],
+            periods=da.time.size + 1,
+            freq=initial_frequency,
+            calendar=da.time.dt.calendar,
+        )
+        # This is the number of days in each sampling period
+        days_per_step = (
+            xr.DataArray(t, dims=("time",), coords={"time": t})
+            .diff("time", label="lower")
+            .dt.days
         )
-    elif initial_frequency_td > pd.Timedelta("1D"):
-        logger.warning(
-            "You appear to be resampling data that is coarser than daily. "
-            "Be aware that this is not currently explicitely supported by xscen and might result in erroneous manipulations."
+        days_per_period = (
+            days_per_step.resample(time=target_frequency)
+            .sum()  # Total number of days per period
+            .sel(time=days_per_step.time, method="ffill")  # Upsample to initial freq
+            .assign_coords(
+                time=days_per_step.time
+            )  # Not sure why we need this, but time coord is from the resample even after sel
         )
 
     if method is None:
@@ -454,7 +469,15 @@ def resample(
             )
 
         # Resample first to find the average wind speed and components
-        ds = ds.resample(time=target_frequency).mean(dim="time", keep_attrs=True)
+        if days_per_step is not None:
+            with xr.set_options(keep_attrs=True):
+                ds = (
+                    (ds * (days_per_step / days_per_period))
+                    .reample(time=target_frequency)
+                    .sum(time="time")
+                )
+        else:
+            ds = ds.resample(time=target_frequency).mean(dim="time", keep_attrs=True)
 
         # Based on Vector Magnitude and Direction equations
         # For example: https://www.khanacademy.org/math/precalculus/x9e81a4f98389efdf:vectors/x9e81a4f98389efdf:component-form/a/vector-magnitude-and-direction-review
@@ -475,6 +498,26 @@ def resample(
         else:
             out = ds[var_name]
 
+    elif days_per_step is not None and method in ["mean", "median", "std", "var"]:
+        if method == "mean":
+            # Avoiding resample().map() is much more performant
+            with xr.set_options(keep_attrs=True):
+                out = (
+                    (da * (days_per_step / days_per_period))
+                    .reample(time=target_frequency)
+                    .sum(time="time")
+                )
+        elif hasattr(xr.core.weighted.DataArrayWeighted, method):
+            ds = xr.merge([da, days_per_step.rename("weights")])
+            da = ds.resample(time="QS-DEC").map(
+                lambda grp: getattr(
+                    grp.drop_vars("weights").weighted(grp.weights), method
+                )(dim="time")
+            )[da.name]
+        else:
+            raise NotImplementedError(
+                f"Weighted resampling not implemented for method {method}."
+            )
     else:
         out = getattr(da.resample(time=target_frequency), method)(
             dim="time", keep_attrs=True

From ce99daa81120e92fc7c6030f1250f1edab9628ab Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Mon, 2 Oct 2023 16:49:47 -0400
Subject: [PATCH 2/8] Add test - missing - doc

---
 docs/goodtoknow.rst   |  8 ++---
 tests/test_extract.py | 75 +++++++++++++++++++++++++++++++++++++++++++
 xscen/extract.py      | 61 +++++++++++++++++++++++++++--------
 3 files changed, 126 insertions(+), 18 deletions(-)

diff --git a/docs/goodtoknow.rst b/docs/goodtoknow.rst
index 1b79028f..6314b133 100644
--- a/docs/goodtoknow.rst
+++ b/docs/goodtoknow.rst
@@ -37,12 +37,12 @@ There are many ways to open data in xscen workflows. The list below tries to mak
 Which function to use when resampling data
 ------------------------------------------
 
-:extract_dataset: :py:func:`~xscen.extract.extract_dataset` has resampling capabilities to provide daily data from finer sources.
+:extract_dataset: :py:func:`~xscen.extract.extract_dataset`'s resampling capabilities are meant to provide daily data from finer sources.
 
-:xclim indicators: Through :py:func:`~xscen.indicators.compute_indicators`, xscen workflows can easily use `xclim indicators <https://xclim.readthedocs.io/en/stable/indicators.html>`_
-    to go from daily data to coarser (monthly, seasonal, annual).
+:resample: :py:func`xscen.extract.resample` extends xarray's `resample` methods with support for weighted resampling when starting from data coarser than daily and for handling of missing timesteps or values.
 
-What is currently not covered by either `xscen` or `xclim` is a method to resample from data coarser than daily, where the base period is non-uniform (ex: resampling from monthly to annual data, taking into account the number of days per month).
+:xclim indicators: Through :py:func:`~xscen.indicators.compute_indicators`, xscen workflows can easily use `xclim indicators <https://xclim.readthedocs.io/en/stable/indicators.html>`_
+    to go from daily data to coarser (monthly, seasonal, annual), with missing values handling. This option will add more metadata than the two firsts.
 
 
 Metadata translation
diff --git a/tests/test_extract.py b/tests/test_extract.py
index 983f6454..bc511621 100644
--- a/tests/test_extract.py
+++ b/tests/test_extract.py
@@ -3,6 +3,7 @@
 import numpy as np
 import pandas as pd
 import pytest
+import xarray as xr
 from conftest import notebooks
 from xclim.testing.helpers import test_timeseries as timeseries
 
@@ -413,3 +414,77 @@ def test_outofrange(self):
 
     def test_none(self):
         assert xs.subset_warming_level(TestSubsetWarmingLevel.ds, wl=20) is None
+
+
+class TestResample:
+    @pytest.mark.parametrize(
+        "infrq,meth,outfrq,exp",
+        [
+            ["MS", "mean", "QS-DEC", [0.47457627, 3, 6.01086957, 9, 11.96666667]],
+            [
+                "QS-DEC",
+                "mean",
+                "YS",
+                [1.49041096, 5.49041096, 9.49453552, 13.49041096, 17.49041096],
+            ],
+            ["MS", "std", "2YS", [6.92009239, 6.91557206]],
+        ],
+    )
+    def test_mean_from_monthly(self, infrq, meth, outfrq, exp):
+        da = timeseries(
+            np.arange(48),
+            variable="tas",
+            start="2001-01-01",
+            freq=infrq,
+        )
+        out = xs.extract.resample(da, outfrq, method=meth)
+        np.testing.assert_allclose(out.isel(time=slice(0, 5)), exp)
+
+    def test_wind_from_monthly(self):
+        uas = timeseries(
+            np.arange(48),
+            variable="uas",
+            start="2001-01-01",
+            freq="MS",
+        )
+        vas = timeseries(
+            np.arange(48),
+            variable="vas",
+            start="2001-01-01",
+            freq="MS",
+        )
+        ds = xr.merge([uas, vas])
+        out = xs.extract.resample(ds.uas, "YS", method="wind_direction", ds=ds)
+        np.testing.assert_allclose(out, [5.5260274, 17.5260274, 29.5260274, 41.5136612])
+
+    def test_missing(self):
+        da = timeseries(
+            np.arange(48),
+            variable="tas",
+            start="2001-01-01",
+            freq="MS",
+        )
+        out = xs.extract.resample(da, "QS-DEC", method="mean", missing="drop")
+        assert out.size == 15
+
+        out = xs.extract.resample(da, "QS-DEC", method="mean", missing="mask")
+        assert out.isel(time=0).isnull().all()
+
+    def test_missing_xclim(self):
+        arr = np.arange(48).astype(float)
+        arr[0] = np.nan
+        arr[40:] = np.nan
+        da = timeseries(
+            arr,
+            variable="tas",
+            start="2001-01-01",
+            freq="MS",
+        )
+        out = xs.extract.resample(da, "YS", method="mean", missing={"method": "any"})
+        assert out.isel(time=0).isnull().all()
+
+        out = xs.extract.resample(
+            da, "YS", method="mean", missing={"method": "pct", "tolerance": 0.6}
+        )
+        assert out.isel(time=0).notnull().all()
+        assert out.isel(time=-1).isnull().all()
diff --git a/xscen/extract.py b/xscen/extract.py
index 222475ba..cd18467c 100644
--- a/xscen/extract.py
+++ b/xscen/extract.py
@@ -368,28 +368,37 @@ def resample(
     *,
     ds: Optional[xr.Dataset] = None,
     method: Optional[str] = None,
+    missing: Union[str, dict] = None,
 ) -> xr.DataArray:
     """Aggregate variable to the target frequency.
 
+    If the input frequency is greater than a week, the resampling operation is weighted by
+    the number of days in each sampling period.
+
     Parameters
     ----------
     da  : xr.DataArray
         DataArray of the variable to resample, must have a "time" dimension and be of a
-        finer temporal resolution than "target_timestep".
+        finer temporal resolution than "target_frequency".
     target_frequency : str
-        The target frequency/freq str, must be one of the frequency supported by pandas.
+        The target frequency/freq str, must be one of the frequency supported by xarray.
     ds : xr.Dataset, optional
         The "wind_direction" resampling method needs extra variables, which can be given here.
     method : {'mean', 'min', 'max', 'sum', 'wind_direction'}, optional
         The resampling method. If None (default), it is guessed from the variable name and frequency,
         using the mapping in CVs/resampling_methods.json. If the variable is not found there,
         "mean" is used by default.
+    missing: {'mask', 'drop'} or dict, optional
+        If 'mask' or 'drop', target periods with fewer steps than expected are masked or dropped.
+        For example, for daily data beginning in January with a `target_frequency` of "QS-DEC". the first season is missing one month.
+        If a dict, it points to a xclim check missing method which will mask periods according to their number of NaN values.
+        The dict must contain a "method" field corresponding to the xclim method name and may contain
+        any other args to pass. Options are documented in :py:mod:`xclim.core.missing`.
 
     Returns
     -------
     xr.DataArray
         Resampled variable
-
     """
     var_name = da.name
 
@@ -397,7 +406,8 @@ def resample(
     initial_frequency_td = pd.Timedelta(
         CV.xrfreq_to_timedelta(xr.infer_freq(da.time.dt.round("T")), None)
     )
-    days_per_step = None
+
+    weights = None
     if initial_frequency != "undetected" and initial_frequency_td > pd.Timedelta(
         7, "days"
     ):
@@ -422,6 +432,7 @@ def resample(
                 time=days_per_step.time
             )  # Not sure why we need this, but time coord is from the resample even after sel
         )
+        weights = days_per_step / days_per_period
 
     if method is None:
         if (
@@ -471,11 +482,7 @@ def resample(
         # Resample first to find the average wind speed and components
         if days_per_step is not None:
             with xr.set_options(keep_attrs=True):
-                ds = (
-                    (ds * (days_per_step / days_per_period))
-                    .reample(time=target_frequency)
-                    .sum(time="time")
-                )
+                ds = (ds * weights).resample(time=target_frequency).sum(dim="time")
         else:
             ds = ds.resample(time=target_frequency).mean(dim="time", keep_attrs=True)
 
@@ -503,13 +510,14 @@ def resample(
             # Avoiding resample().map() is much more performant
             with xr.set_options(keep_attrs=True):
                 out = (
-                    (da * (days_per_step / days_per_period))
-                    .reample(time=target_frequency)
-                    .sum(time="time")
+                    (da * weights)
+                    .resample(time=target_frequency)
+                    .sum(dim="time")
+                    .rename(da.name)
                 )
         elif hasattr(xr.core.weighted.DataArrayWeighted, method):
-            ds = xr.merge([da, days_per_step.rename("weights")])
-            da = ds.resample(time="QS-DEC").map(
+            ds = xr.merge([da, weights.rename("weights")])
+            out = ds.resample(time=target_frequency).map(
                 lambda grp: getattr(
                     grp.drop_vars("weights").weighted(grp.weights), method
                 )(dim="time")
@@ -523,6 +531,31 @@ def resample(
             dim="time", keep_attrs=True
         )
 
+    if missing in ["mask", "drop"]:
+        steps_per_period = (
+            xr.ones_like(da.time, dtype="int").resample(time=target_frequency).sum()
+        )
+        t = xr.date_range(
+            steps_per_period.indexes["time"][0],
+            periods=steps_per_period.time.size + 1,
+            freq=target_frequency,
+        )
+        expected = (
+            xr.DataArray(t, dims=("time",), coords={"time": t}).diff(
+                "time", label="lower"
+            )
+            / initial_frequency_td
+        )
+        complete = (steps_per_period / expected) > 0.95
+    elif isinstance(missing, dict):
+        missmeth = missing.pop("method")
+        complete = ~xc.core.missing.MISSING_METHODS[missmeth](
+            da, target_frequency, initial_frequency
+        )(**missing)
+        missing = "mask"
+    if missing in {"mask", "drop"}:
+        out = out.where(complete, drop=(missing == "drop"))
+
     new_history = (
         f"[{datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')}] {method} "
         f"resample from {initial_frequency} to {target_frequency} - xarray v{xr.__version__}"

From 98c37e01e8ea81bd76915b154c63e4ceeb4833db Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Mon, 2 Oct 2023 16:52:57 -0400
Subject: [PATCH 3/8] upd hist

---
 HISTORY.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/HISTORY.rst b/HISTORY.rst
index c23cba17..a8bae594 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -16,6 +16,7 @@ New features and enhancements
 * ``xs.spatial_mean`` now accepts the ``region="global"`` keyword to perform a global average (:issue:`94`, :pull:`260`).
 * ``xs.spatial_mean`` with ``method='xESMF'`` will also automatically segmentize polygons (down to a 1° resolution) to ensure a correct average (:pull:`260`).
 * Added documentation for `require_all_on` in `search_data_catalogs`. (:pull:`263`).
+* Better ``xs.extract.resample`` : support for weighted resampling operations when starting with frequencies coarser than daily and missing timesteps/values handling. (:issue:`80`, :issue:`93`, :pull:`265`).
 
 Breaking changes
 ^^^^^^^^^^^^^^^^

From 5ae762dc01e1b08e48648c043df9453eef874b55 Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Mon, 2 Oct 2023 17:05:56 -0400
Subject: [PATCH 4/8] implement weighted median

---
 xscen/extract.py | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/xscen/extract.py b/xscen/extract.py
index cd18467c..a2534e52 100644
--- a/xscen/extract.py
+++ b/xscen/extract.py
@@ -515,17 +515,15 @@ def resample(
                     .sum(dim="time")
                     .rename(da.name)
                 )
-        elif hasattr(xr.core.weighted.DataArrayWeighted, method):
+        else:
+            kws = {"q": 0.5} if method == "median" else {}
             ds = xr.merge([da, weights.rename("weights")])
             out = ds.resample(time=target_frequency).map(
                 lambda grp: getattr(
-                    grp.drop_vars("weights").weighted(grp.weights), method
-                )(dim="time")
+                    grp.drop_vars("weights").weighted(grp.weights),
+                    method if method != "median" else "quantile",
+                )(dim="time", **kws)
             )[da.name]
-        else:
-            raise NotImplementedError(
-                f"Weighted resampling not implemented for method {method}."
-            )
     else:
         out = getattr(da.resample(time=target_frequency), method)(
             dim="time", keep_attrs=True

From 01cf1539fac60619c408b5940b9eb22fe8fd315d Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Mon, 2 Oct 2023 17:26:05 -0400
Subject: [PATCH 5/8] add metadata

---
 tests/test_extract.py | 10 ++++++++--
 xscen/extract.py      | 21 ++++++++++++++++-----
 2 files changed, 24 insertions(+), 7 deletions(-)

diff --git a/tests/test_extract.py b/tests/test_extract.py
index bc511621..cd36f34a 100644
--- a/tests/test_extract.py
+++ b/tests/test_extract.py
@@ -428,9 +428,15 @@ class TestResample:
                 [1.49041096, 5.49041096, 9.49453552, 13.49041096, 17.49041096],
             ],
             ["MS", "std", "2YS", [6.92009239, 6.91557206]],
+            [
+                "QS",
+                "median",
+                "YS",
+                [1.516437, 5.516437, 9.516437, 13.51092864, 17.516437],
+            ],
         ],
     )
-    def test_mean_from_monthly(self, infrq, meth, outfrq, exp):
+    def test_weighted(self, infrq, meth, outfrq, exp):
         da = timeseries(
             np.arange(48),
             variable="tas",
@@ -440,7 +446,7 @@ def test_mean_from_monthly(self, infrq, meth, outfrq, exp):
         out = xs.extract.resample(da, outfrq, method=meth)
         np.testing.assert_allclose(out.isel(time=slice(0, 5)), exp)
 
-    def test_wind_from_monthly(self):
+    def test_weighted_wind(self):
         uas = timeseries(
             np.arange(48),
             variable="uas",
diff --git a/xscen/extract.py b/xscen/extract.py
index a2534e52..4c55f931 100644
--- a/xscen/extract.py
+++ b/xscen/extract.py
@@ -408,8 +408,10 @@ def resample(
     )
 
     weights = None
-    if initial_frequency != "undetected" and initial_frequency_td > pd.Timedelta(
-        7, "days"
+    if (
+        initial_frequency != "undetected"
+        and initial_frequency_td > pd.Timedelta(7, "days")
+        and method in ["mean", "median", "std", "var"]
     ):
         # More than a week -> non-uniform sampling length!
         t = xr.date_range(
@@ -505,7 +507,7 @@ def resample(
         else:
             out = ds[var_name]
 
-    elif days_per_step is not None and method in ["mean", "median", "std", "var"]:
+    elif weights is not None:
         if method == "mean":
             # Avoiding resample().map() is much more performant
             with xr.set_options(keep_attrs=True):
@@ -529,6 +531,7 @@ def resample(
             dim="time", keep_attrs=True
         )
 
+    missing_note = " "
     if missing in ["mask", "drop"]:
         steps_per_period = (
             xr.ones_like(da.time, dtype="int").resample(time=target_frequency).sum()
@@ -545,18 +548,26 @@ def resample(
             / initial_frequency_td
         )
         complete = (steps_per_period / expected) > 0.95
+        action = "masking" if missing == "mask" else "dropping"
+        missing_note = f", {action} incomplete periods "
     elif isinstance(missing, dict):
         missmeth = missing.pop("method")
         complete = ~xc.core.missing.MISSING_METHODS[missmeth](
             da, target_frequency, initial_frequency
         )(**missing)
+        funcstr = xc.core.formatting.gen_call_string(
+            f"xclim.core.missing_{missmeth}", **missing
+        )
         missing = "mask"
+        missing_note = f", masking incomplete periods according to {funcstr} "
     if missing in {"mask", "drop"}:
         out = out.where(complete, drop=(missing == "drop"))
 
     new_history = (
-        f"[{datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')}] {method} "
-        f"resample from {initial_frequency} to {target_frequency} - xarray v{xr.__version__}"
+        f"[{datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')}] "
+        f"{'weighted' if weights is not None else ''} {method} "
+        f"resample from {initial_frequency} to {target_frequency}"
+        f"{missing_note}- xarray v{xr.__version__}"
     )
     history = (
         new_history + " \n " + out.attrs["history"]

From f5d338fe0bfc1d26e3db03b9ec0b9fe582750a88 Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Thu, 5 Oct 2023 16:45:16 -0400
Subject: [PATCH 6/8] Update xscen/extract.py

Co-authored-by: RondeauG <38501935+RondeauG@users.noreply.github.com>
---
 xscen/extract.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/xscen/extract.py b/xscen/extract.py
index 4c55f931..da2c4c32 100644
--- a/xscen/extract.py
+++ b/xscen/extract.py
@@ -389,8 +389,8 @@ def resample(
         using the mapping in CVs/resampling_methods.json. If the variable is not found there,
         "mean" is used by default.
     missing: {'mask', 'drop'} or dict, optional
-        If 'mask' or 'drop', target periods with fewer steps than expected are masked or dropped.
-        For example, for daily data beginning in January with a `target_frequency` of "QS-DEC". the first season is missing one month.
+        If 'mask' or 'drop', target periods that would have been computed from fewer timesteps than expected are masked or dropped, using a threshold of 5% of missing data.
+        For example, the first season of a `target_frequency` of "QS-DEC" will be masked or dropped if data only starts in January.
         If a dict, it points to a xclim check missing method which will mask periods according to their number of NaN values.
         The dict must contain a "method" field corresponding to the xclim method name and may contain
         any other args to pass. Options are documented in :py:mod:`xclim.core.missing`.

From 8e8d49a1517422e68bb9e72367382186985c797b Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Thu, 5 Oct 2023 17:00:25 -0400
Subject: [PATCH 7/8] Better xrfreq2timedelta - fix tests - warn

---
 xscen/extract.py | 62 ++++++++++++++++++++++++++----------------------
 xscen/utils.py   |  6 +++++
 2 files changed, 39 insertions(+), 29 deletions(-)

diff --git a/xscen/extract.py b/xscen/extract.py
index da2c4c32..b2ee1382 100644
--- a/xscen/extract.py
+++ b/xscen/extract.py
@@ -14,6 +14,7 @@
 import xarray as xr
 import xclim as xc
 from intake_esm.derived import DerivedVariableRegistry
+from xclim.core.calendar import compare_offsets
 
 from .catalog import DataCatalog  # noqa
 from .catalog import ID_COLUMNS, concat_data_catalogs, generate_id, subset_file_coverage
@@ -23,7 +24,7 @@
 from .spatial import subset
 from .utils import CV
 from .utils import ensure_correct_time as _ensure_correct_time
-from .utils import get_cat_attrs, natural_sort, standardize_periods
+from .utils import get_cat_attrs, natural_sort, standardize_periods, xrfreq_to_timedelta
 
 logger = logging.getLogger(__name__)
 
@@ -403,15 +404,36 @@ def resample(
     var_name = da.name
 
     initial_frequency = xr.infer_freq(da.time.dt.round("T")) or "undetected"
-    initial_frequency_td = pd.Timedelta(
-        CV.xrfreq_to_timedelta(xr.infer_freq(da.time.dt.round("T")), None)
-    )
+    if initial_frequency == "undetected":
+        warnings.warn(
+            "Could not infer the frequency of the dataset. Be aware that this might result in erroneous manipulations."
+        )
+
+    if method is None:
+        if (
+            target_frequency in CV.resampling_methods.dict
+            and var_name in CV.resampling_methods.dict[target_frequency]
+        ):
+            method = CV.resampling_methods(target_frequency)[var_name]
+            logger.info(
+                f"Resampling method for {var_name}: '{method}', based on variable name and frequency."
+            )
+
+        elif var_name in CV.resampling_methods.dict["any"]:
+            method = CV.resampling_methods("any")[var_name]
+            logger.info(
+                f"Resampling method for {var_name}: '{method}', based on variable name."
+            )
+
+        else:
+            method = "mean"
+            logger.info(f"Resampling method for {var_name} defaulted to: 'mean'.")
 
     weights = None
     if (
         initial_frequency != "undetected"
-        and initial_frequency_td > pd.Timedelta(7, "days")
-        and method in ["mean", "median", "std", "var"]
+        and compare_offsets(initial_frequency, ">", "W")
+        and method in ["mean", "median", "std", "var", "wind_direction"]
     ):
         # More than a week -> non-uniform sampling length!
         t = xr.date_range(
@@ -436,26 +458,6 @@ def resample(
         )
         weights = days_per_step / days_per_period
 
-    if method is None:
-        if (
-            target_frequency in CV.resampling_methods.dict
-            and var_name in CV.resampling_methods.dict[target_frequency]
-        ):
-            method = CV.resampling_methods(target_frequency)[var_name]
-            logger.info(
-                f"Resampling method for {var_name}: '{method}', based on variable name and frequency."
-            )
-
-        elif var_name in CV.resampling_methods.dict["any"]:
-            method = CV.resampling_methods("any")[var_name]
-            logger.info(
-                f"Resampling method for {var_name}: '{method}', based on variable name."
-            )
-
-        else:
-            method = "mean"
-            logger.info(f"Resampling method for {var_name} defaulted to: 'mean'.")
-
     # TODO : Support non-surface wind?
     if method == "wind_direction":
         ds[var_name] = da
@@ -482,7 +484,7 @@ def resample(
             )
 
         # Resample first to find the average wind speed and components
-        if days_per_step is not None:
+        if weights is not None:
             with xr.set_options(keep_attrs=True):
                 ds = (ds * weights).resample(time=target_frequency).sum(dim="time")
         else:
@@ -532,7 +534,8 @@ def resample(
         )
 
     missing_note = " "
-    if missing in ["mask", "drop"]:
+    initial_td = xrfreq_to_timedelta(initial_frequency)
+    if missing in ["mask", "drop"] and not pd.isnull(initial_td):
         steps_per_period = (
             xr.ones_like(da.time, dtype="int").resample(time=target_frequency).sum()
         )
@@ -541,11 +544,12 @@ def resample(
             periods=steps_per_period.time.size + 1,
             freq=target_frequency,
         )
+
         expected = (
             xr.DataArray(t, dims=("time",), coords={"time": t}).diff(
                 "time", label="lower"
             )
-            / initial_frequency_td
+            / initial_td
         )
         complete = (steps_per_period / expected) > 0.95
         action = "masking" if missing == "mask" else "dropping"
diff --git a/xscen/utils.py b/xscen/utils.py
index 0ad5a790..c7d931a8 100644
--- a/xscen/utils.py
+++ b/xscen/utils.py
@@ -1274,3 +1274,9 @@ def season_sort_key(idx: pd.Index, name: str = None):
         # TypeError if season element was not a string.
         pass
     return idx
+
+
+def xrfreq_to_timedelta(freq):
+    """Approximate the length of a period based on its frequency offset."""
+    N, B, _, _ = parse_offset(freq)
+    return N * pd.Timedelta(CV.xrfreq_to_timedelta(B, "NaT"))

From 1f3de06bbe1118b810542763051a4edac2d15193 Mon Sep 17 00:00:00 2001
From: "bumpversion[bot]" <bumpversion[bot]@ouranos.ca>
Date: Mon, 9 Oct 2023 14:37:52 +0000
Subject: [PATCH 8/8] =?UTF-8?q?Bump=20version:=200.7.12-beta=20=E2=86=92?=
 =?UTF-8?q?=200.7.13-beta?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 .cruft.json         | 2 +-
 setup.cfg           | 2 +-
 setup.py            | 2 +-
 tests/test_xscen.py | 2 +-
 xscen/__init__.py   | 2 +-
 5 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/.cruft.json b/.cruft.json
index 0cdc4aa8..19851606 100644
--- a/.cruft.json
+++ b/.cruft.json
@@ -11,7 +11,7 @@
       "project_slug": "xscen",
       "project_short_description": "A climate change scenario-building analysis framework, built with xclim/xarray.",
       "pypi_username": "RondeauG",
-      "version": "0.7.12-beta",
+      "version": "0.7.13-beta",
       "use_pytest": "y",
       "use_black": "y",
       "add_pyup_badge": "n",
diff --git a/setup.cfg b/setup.cfg
index c7973609..585a0776 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.7.12-beta
+current_version = 0.7.13-beta
 commit = True
 tag = False
 parse = (?P<major>\d+)\.(?P<minor>\d+).(?P<patch>\d+)(\-(?P<release>[a-z]+))?
diff --git a/setup.py b/setup.py
index 145a0812..61d18c4a 100644
--- a/setup.py
+++ b/setup.py
@@ -102,6 +102,6 @@ def run(self):
     test_suite="tests",
     extras_require={"dev": dev_requirements},
     url="https://github.com/Ouranosinc/xscen",
-    version="0.7.12-beta",
+    version="0.7.13-beta",
     zip_safe=False,
 )
diff --git a/tests/test_xscen.py b/tests/test_xscen.py
index 9a3b5caf..4b35f02a 100644
--- a/tests/test_xscen.py
+++ b/tests/test_xscen.py
@@ -28,4 +28,4 @@ def test_package_metadata(self):
         contents = f.read()
         assert """Gabriel Rondeau-Genesse""" in contents
         assert '__email__ = "rondeau-genesse.gabriel@ouranos.ca"' in contents
-        assert '__version__ = "0.7.12-beta"' in contents
+        assert '__version__ = "0.7.13-beta"' in contents
diff --git a/xscen/__init__.py b/xscen/__init__.py
index d6476feb..896a42db 100644
--- a/xscen/__init__.py
+++ b/xscen/__init__.py
@@ -52,7 +52,7 @@
 
 __author__ = """Gabriel Rondeau-Genesse"""
 __email__ = "rondeau-genesse.gabriel@ouranos.ca"
-__version__ = "0.7.12-beta"
+__version__ = "0.7.13-beta"
 
 
 # monkeypatch so that warnings.warn() doesn't mention itself