Repository: Open-EO/FuseTS (language: python)

Issue:
SRD-05 Multi-temporal outlier filter
Implement an outlier filter based on single pixel timeseries.
Should we also consider elimination of false positives, like Sentinel-2 flagging water as cloud?


Discussion / hints:
I had a bit of a look at this. I feel the scope of this issue is quite undefined:
- what kind of algorithm do we wish to apply?
- what kind of data are we working with?
- what kind of time-series are we working with (which use case)?

Already the 3 questions above make the scope of this issue explode. 

There was also one article posted next to this issue which might come in handy: https://www.sciencedirect.com/science/article/pii/S0034425714002259, but I don't believe that one to be useful, as it is targeted to Landsat, which we have never used so far in the project. 

I would suggest that if a user has a specific problem, the user implements their own outlier filter as a UDF, suitable to their problem, as opposed to us loosing time by providing different methods of outlier filtering as a preventive measure, not really being sure it would be used or not.

So i would propose going with something simple and generic, such as the z-score outlier filtering, which is basically quantitatively evaluating which points around the trend line are potential outliers, based on the size of the standard deviation from the mean. This is something intuitive and can be applied to take care of the occasional sharp spike.
As discussed during the meeting, the proposed approach is the preferred way to go.
Adding some more information: 
- Example solution: https://stackoverflow.com/a/75938737
- will follow the peak-valley example: https://github.com/Open-EO/FuseTS/blob/main/src/fusets/peakvalley.py



Fix patch:
diff --git a/src/fusets/temporal_outliers.py b/src/fusets/temporal_outliers.py
new file mode 100644
index 0000000..3bf9c7a
--- /dev/null
+++ b/src/fusets/temporal_outliers.py
@@ -0,0 +1,70 @@
+from datetime import datetime
+from typing import List, Sequence, Union
+
+import numpy as np
+import pandas as pd
+import xarray
+from xarray import DataArray
+
+from fusets._xarray_utils import _extract_dates, _time_dimension
+
+
+def temporal_outliers(
+    array: DataArray, window: Union[int, str], threshold: float, variables: List[str] = None
+) -> DataArray:
+    """
+    Algorithm for a z-score-based filtering of time series outliers
+
+    Args:
+        array: input data array
+        window: pandas-based window-size, can be integer or string, e.g., '20D' for a 20 days window
+        threshold: the threshold to be applied on the z-scores for filtering outlier (amount of st. dev.)
+        variables: The list of variable names that should be affected, or None to use all variables
+
+    Returns:
+        data array with outliers filtered out and with the rolling mean of time series
+    """
+
+    dates = np.array(_extract_dates(array))
+    time_dimension = _time_dimension(array, None)
+
+    if variables is not None:
+        array = array.drop_vars([var for var in list(array.data_vars) if var not in variables])
+
+    def callback(timeseries):
+        return temporal_outliers_f(dates, timeseries, window, threshold)
+
+    result = xarray.apply_ufunc(
+        callback,
+        array,
+        input_core_dims=[[time_dimension]],
+        output_core_dims=[[time_dimension]],
+        vectorize=True,
+    )
+
+    return result
+
+
+def temporal_outliers_f(x: Sequence[datetime], y: np.ndarray, window: Union[int, str], threshold: float) -> np.ndarray:
+    """
+    Algorithm for a z-score-based filtering of time series outliers
+
+    Args:
+        x: array of timestamps
+        y: array of input feature values
+        window: pandas-based window-size, can be integer or string, e.g., '20D' for a 20 days window
+        threshold: the threshold to be applied on the z-scores for filtering outlier (amount of st. dev.)
+
+    Returns:
+        array with outliers filtered out and replaced with the rolling mean of time series
+    """
+
+    timeseries = pd.Series(y, index=x)
+    rolling = timeseries.rolling(window=window, center=True, closed="both")
+    ts_mean = rolling.mean()
+    ts_std = rolling.std(ddof=1)
+
+    ts_zscore = timeseries.sub(ts_mean).div(ts_std)
+    ts_mask = ts_zscore.between(-threshold, threshold)
+
+    return timeseries.where(ts_mask, ts_mean).to_numpy(dtype='float32')


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 665fd5c..d1e0b95 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,9 +1,11 @@
 from datetime import datetime, timedelta
 
 import numpy as np
+import pandas as pd
 import pytest
 import xarray
 
+
 @pytest.fixture
 def connection():
     """Connection fixture to a backend of given version with some image collections."""
@@ -61,6 +63,15 @@ def generate_data(xs: np.array):
                      dims=["time"],
                      coords=dict(time=ts))
 
+@pytest.fixture
+def outlier_timeseries():
+    rng = np.random.default_rng(42)
+    dates = pd.date_range('2019-01-01', '2019-12-31', periods=300)
+    values = np.sin(np.linspace(0, 4*np.pi, len(dates))) + rng.random(len(dates))*0.2
+    values[rng.choice(range(len(dates)), 4).astype(int)] += rng.choice([-1, 1], 4) * 5
+
+    return xarray.DataArray(data=values, dims=["time"], coords=dict(time=dates))
+
 @pytest.fixture
 def areas():
     return {
diff --git a/tests/test_temporal_outliers.py b/tests/test_temporal_outliers.py
new file mode 100644
index 0000000..22bf6bd
--- /dev/null
+++ b/tests/test_temporal_outliers.py
@@ -0,0 +1,16 @@
+import numpy as np
+
+from fusets._xarray_utils import _extract_dates
+from fusets.temporal_outliers import temporal_outliers_f
+from numpy.testing import assert_almost_equal
+
+
+def test_temporal_outlier_filtering(outlier_timeseries):
+    dates = np.array(_extract_dates(outlier_timeseries))
+    vals = outlier_timeseries.values
+
+    filtered_vals = temporal_outliers_f(dates, vals, window='20D', threshold=3)
+
+    assert_almost_equal(filtered_vals.mean(), 0.09904716, decimal=6)
+    assert_almost_equal(filtered_vals.std(), 0.71552783, decimal=6)
+