Repository: Ouranosinc/xclim

Issue:
Improve missing data check

### Description
Currently missing data checks in indicator classes are very strict (single missing day or single nan value) will result in nan for a given resample time step


Options include:

1.  Percent missing : https://github.com/Ouranosinc/xclim/issues/364#issuecomment-585804628
2. WMO data completeness https://github.com/Ouranosinc/xclim/issues/364#issuecomment-585798512


For option 1 need to decide on :

- How % missing data is determined on a annual frequency : see https://github.com/Ouranosinc/xclim/issues/364#issuecomment-585828201

- Whether missing data is interpolated or not before averaging / summing : https://github.com/Ouranosinc/xclim/issues/364#issuecomment-585838629


Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index f0f5444e..6883462f 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -15,6 +15,8 @@ History
 * Ensembles : Allow alignment of datasets with same frequency but different offsets.
 * Bug fixes in run_length for run-with-dates methods when the date is not found in the run.
 * Remove deepcopy from subset.subset_shape to improve memory usage.
+* Add `missing_wmo` function, identifying null calculations based on criteria from WMO.
+* Add `missing_pct` function, identifying null calculations based on percentage of missing values.
 
 0.15.x (2020-03-12)
 -------------------
diff --git a/xclim/core/checks.py b/xclim/core/checks.py
index 50805aaa..fb7eef91 100644
--- a/xclim/core/checks.py
+++ b/xclim/core/checks.py
@@ -161,8 +161,110 @@ def check_is_dataarray(comp):
 
 
 # This function can probably be made simpler once CFPeriodIndex is implemented.
+class MissingBase:
+    def __init__(self, da, freq, **indexer):
+        self.null, self.count = self.prepare(da, freq, **indexer)
+
+    @staticmethod
+    def split_freq(freq):
+        if "-" in freq:
+            pfreq, anchor = freq.split("-")
+        else:
+            pfreq, anchor = freq, None
+
+        return pfreq, anchor
+
+    @staticmethod
+    def is_null(da, freq, **indexer):
+        # Compute the number of days in the time series during each period at the given frequency.
+        selected = generic.select_time(da, **indexer)
+        if selected.time.size == 0:
+            raise ValueError("No data for selected period.")
+
+        return selected.isnull().resample(time=freq)
+
+    def prepare(self, da, freq, **indexer):
+        pfreq, anchor = self.split_freq(freq)
+        null = self.is_null(da, freq, **indexer)
+
+        c = null.sum(dim="time")
+
+        # Otherwise simply use the start and end dates to find the expected number of days.
+        if pfreq.endswith("S"):
+            start_time = c.indexes["time"]
+            end_time = start_time.shift(1, freq=freq)
+        else:
+            end_time = c.indexes["time"]
+            start_time = end_time.shift(-1, freq=freq)
+
+        if indexer:
+            # Create a full synthetic time series and compare the number of days with the original series.
+            t0 = str(start_time[0].date())
+            t1 = str(end_time[-1].date())
+            if isinstance(c.indexes["time"], xr.CFTimeIndex):
+                cal = da.time.encoding.get("calendar")
+                t = xr.cftime_range(t0, t1, freq="D", calendar=cal)
+            else:
+                t = pd.date_range(t0, t1, freq="D")
+
+            sda = xr.DataArray(data=np.ones(len(t)), coords={"time": t}, dims=("time",))
+            st = generic.select_time(sda, **indexer)
+            count = st.notnull().resample(time=freq).sum(dim="time")
+
+        else:
+            n = (end_time - start_time).days
+            count = xr.DataArray(n.values, coords={"time": c.time}, dims="time")
+
+        return null, count
+
+    def is_missing(self, null, count, **kwargs):
+        """Return whether or not the values within each period should be considered missing or not."""
+        raise NotImplementedError
+
+    def __call__(self, **kwargs):
+        return self.is_missing(self.null, self.count, **kwargs)
+
+
+class MissingAny(MissingBase):
+    def is_missing(self, null, count, **kwargs):
+        cond0 = null.count(dim="time") != count  # Check total number of days
+        cond1 = null.sum(dim="time") > 0  # Check if any is missing
+        return cond0 | cond1
+
+
+class MissingWMO(MissingAny):
+    def __init__(self, da, freq, **indexer):
+        # Force computation on monthly frequency
+        if not freq.startswith("M"):
+            raise ValueError
+        super().__init__(da, freq, **indexer)
+
+    def is_missing(self, null, count, nm=11, nc=5):
+        import xclim.indices.run_length as rl
+
+        # Check total number of days
+        cond0 = null.count(dim="time") != count
+
+        # Check if more than threshold is missing
+        cond1 = null.sum(dim="time") >= nm
+
+        # Check for consecutive missing values
+        cond2 = null.map(rl.longest_run, dim="time") >= nc
+
+        return cond0 | cond1 | cond2
+
+
+class MissingPct(MissingBase):
+    def is_missing(self, null, count, tolerance=0.1):
+        if tolerance < 0 or tolerance > 1:
+            raise ValueError("tolerance should be between 0 and 1.")
+
+        n = count - null.count(dim="time") + null.sum(dim="time")
+        return n / count >= tolerance
+
+
 def missing_any(da, freq, **indexer):
-    r"""Return a boolean DataArray indicating whether there are missing days in the resampled array.
+    r"""Return whether there are missing days in the array.
 
     Parameters
     ----------
@@ -178,44 +280,66 @@ def missing_any(da, freq, **indexer):
     Returns
     -------
     out : DataArray
-      A boolean array set to True if any month or year has missing values.
+      A boolean array set to True if period has missing values.
     """
-    if "-" in freq:
-        pfreq, anchor = freq.split("-")
-    else:
-        pfreq = freq
-
-    # Compute the number of days in the time series during each period at the given frequency.
-    selected = generic.select_time(da, **indexer)
-    if selected.time.size == 0:
-        raise ValueError("No data for selected period.")
-
-    c = selected.notnull().resample(time=freq).sum(dim="time")
-
-    # Otherwise simply use the start and end dates to find the expected number of days.
-    if pfreq.endswith("S"):
-        start_time = c.indexes["time"]
-        end_time = start_time.shift(1, freq=freq)
-    else:
-        end_time = c.indexes["time"]
-        start_time = end_time.shift(-1, freq=freq)
-
-    if indexer:
-        # Create a full synthetic time series and compare the number of days with the original series.
-        t0 = str(start_time[0].date())
-        t1 = str(end_time[-1].date())
-        if isinstance(c.indexes["time"], xr.CFTimeIndex):
-            cal = da.time.encoding.get("calendar")
-            t = xr.cftime_range(t0, t1, freq="D", calendar=cal)
-        else:
-            t = pd.date_range(t0, t1, freq="D")
+    return MissingAny(da, freq, **indexer)()
 
-        sda = xr.DataArray(data=np.ones(len(t)), coords={"time": t}, dims=("time",))
-        st = generic.select_time(sda, **indexer)
-        sn = st.notnull().resample(time=freq).sum(dim="time")
-        miss = sn != c
-        return miss
 
-    n = (end_time - start_time).days
-    nda = xr.DataArray(n.values, coords={"time": c.time}, dims="time")
-    return c != nda
+def missing_wmo(da, freq, nm=11, nc=5, **indexer):
+    r"""Return whether a series fails WMO criteria for missing days.
+
+    The World Meteorological Organisation recommends that where monthly means are computed from daily values,
+    it should considered missing if either of these two criteria are met:
+      – observations are missing for 11 or more days during the month;
+      – observations are missing for a period of 5 or more consecutive days during the month.
+
+    Stricter criteria are sometimes used in practice, with a tolerance of 5 missing values or 3 consecutive missing
+    values.
+
+    Parameters
+    ----------
+    da : DataArray
+      Input array at daily frequency.
+    freq : str
+      Resampling frequency.
+    nm : int
+      Number of missing values per month that should not be exceeded.
+    nc : int
+      Number of consecutive missing values per month that should not be exceeded.
+    **indexer : {dim: indexer, }, optional
+      Time attribute and values over which to subset the array. For example, use season='DJF' to select winter values,
+      month=1 to select January, or month=[6,7,8] to select summer months. If not indexer is given, all values are
+      considered.
+
+    Returns
+    -------
+    out : DataArray
+      A boolean array set to True if period has missing values.
+    """
+    missing = MissingWMO(da, "M", **indexer)(nm=nm, nc=nc)
+    return missing.resample(time=freq).any()
+
+
+def missing_pct(da, freq, tolerance, **indexer):
+    r"""Return whether there are more missing days in the array than a given percentage.
+
+    Parameters
+    ----------
+    da : DataArray
+      Input array at daily frequency.
+    freq : str
+      Resampling frequency.
+    tolerance : float
+      Fraction of missing values that is tolerated.
+    **indexer : {dim: indexer, }, optional
+      Time attribute and values over which to subset the array. For example, use season='DJF' to select winter
+      values,
+      month=1 to select January, or month=[6,7,8] to select summer months. If not indexer is given, all values are
+      considered.
+
+    Returns
+    -------
+    out : DataArray
+      A boolean array set to True if period has missing values.
+    """
+    return MissingPct(da, freq, **indexer)(tolerance=tolerance)
diff --git a/xclim/indices/run_length.py b/xclim/indices/run_length.py
index 3e0866da..0c50b132 100644
--- a/xclim/indices/run_length.py
+++ b/xclim/indices/run_length.py
@@ -95,7 +95,7 @@ def longest_run(
     ufunc_1dim : Union[str, bool]
       Use the 1d 'ufunc' version of this function : default (auto) will attempt to select optimal
       usage based on number of data points.  Using 1D_ufunc=True is typically more efficient
-      for dataarray with a small number of gridpoints.
+      for DataArray with a small number of grid points.
     Returns
     -------
     N-dimensional array (int)


Test patch:
diff --git a/tests/test_checks.py b/tests/test_checks.py
index c00e22ff..3cae7161 100644
--- a/tests/test_checks.py
+++ b/tests/test_checks.py
@@ -130,3 +130,34 @@ class TestMissingAnyFills:
 
         miss = checks.missing_any(ds.q_sim, freq="YS", season="JJA")
         np.testing.assert_array_equal(miss, False)
+
+
+class TestMissingWMO:
+    def test_missing_days(self, tas_series):
+        a = np.arange(360.0)
+        a[5:7] = np.nan  # Number of missing values under acceptable limit in a month
+        a[40:45] = np.nan  # Too many consecutive missing values
+        a[70:92:2] = np.nan  # Too many non-consecutive missing values
+        ts = tas_series(a)
+        out = checks.missing_wmo(ts, freq="MS")
+        assert not out[0]
+        assert out[1]
+        assert out[2]
+
+    def test_missing_days_in_quarter(self, tas_series):
+        a = np.arange(360.0)
+        a[5:16] = np.nan  # Number of missing values under acceptable limit in a month
+        ts = tas_series(a)
+        out = checks.missing_wmo(ts, freq="QS-JAN")
+        np.testing.assert_array_equal(out, [True, False, False, True])
+
+
+class TestMissingPct:
+    def test_missing_days(self, tas_series):
+        a = np.arange(360.0)
+        a[5:7] = np.nan  # Number of missing values under acceptable limit in a month
+        a[40:45] = np.nan  # Too many missing values
+        ts = tas_series(a)
+        out = checks.missing_pct(ts, freq="MS", tolerance=0.1)
+        assert not out[0]
+        assert out[1]