Repository: solo-spice/sospice (language: python)

Issue:
Merge develop



Discussion / hints:
## [Codecov](https://app.codecov.io/gh/solo-spice/sospice/pull/34?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None) Report
Patch coverage is **`100.00%`** of modified lines.

| [Files Changed](https://app.codecov.io/gh/solo-spice/sospice/pull/34?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None) | Coverage |
|---|---|
| [sospice/calibrate/uncertainties.py](https://app.codecov.io/gh/solo-spice/sospice/pull/34?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None#diff-c29zcGljZS9jYWxpYnJhdGUvdW5jZXJ0YWludGllcy5weQ==) | `100.00%` |
| [sospice/catalog/catalog.py](https://app.codecov.io/gh/solo-spice/sospice/pull/34?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None#diff-c29zcGljZS9jYXRhbG9nL2NhdGFsb2cucHk=) | `100.00%` |
| [sospice/catalog/file\_metadata.py](https://app.codecov.io/gh/solo-spice/sospice/pull/34?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None#diff-c29zcGljZS9jYXRhbG9nL2ZpbGVfbWV0YWRhdGEucHk=) | `100.00%` |
| [sospice/instrument\_modelling/observation.py](https://app.codecov.io/gh/solo-spice/sospice/pull/34?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None#diff-c29zcGljZS9pbnN0cnVtZW50X21vZGVsbGluZy9vYnNlcnZhdGlvbi5weQ==) | `100.00%` |


:loudspeaker: Thoughts on this report? [Let us know!](https://about.codecov.io/pull-request-comment-report/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None).


Fix patch:
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index d859114..d7c3ee1 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -37,7 +37,9 @@ jobs:
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest
       run: |
-        python -m setuptools_scm
+        python -c "import setuptools_scm; version = setuptools_scm.get_version(); setuptools_scm.dump_version('sospice', version, '_version.py')"
+        ls -l _version.py || echo "_version.py not found"
+        ls -l sospice/_version.py ||  echo "sospice/_version.py not found"
         pytest --cov --cov-report=xml
     - name: Upload coverage reports to Codecov with GitHub Action
       uses: codecov/codecov-action@v3
diff --git a/changelog/35.feature.rst b/changelog/35.feature.rst
new file mode 100644
index 0000000..d2bff04
--- /dev/null
+++ b/changelog/35.feature.rst
@@ -0,0 +1,1 @@
+Add local sigma clipping
diff --git a/changelog/README.rst b/changelog/README.rst
index d3a5eef..d0dd4e9 100644
--- a/changelog/README.rst
+++ b/changelog/README.rst
@@ -12,7 +12,7 @@ The ``CHANGELOG`` will be read by users, so this description should be aimed at
 
 Make sure to use full sentences with correct case and punctuation, for example::
 
-    Add support for PSF correction`.
+    Add support for PSF correction.
 
 Please try to use Sphinx intersphinx using backticks.
 
diff --git a/docs/source/index.rst b/docs/source/index.rst
index d5d06a9..9c77879 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -31,6 +31,7 @@ Documentation for this package is available on `Read the Docs <https://sospice.r
    getting_started
    calibration
    catalog
+   utilities
    citation
    api
    develop
diff --git a/docs/source/utilities.rst b/docs/source/utilities.rst
new file mode 100644
index 0000000..fd2addb
--- /dev/null
+++ b/docs/source/utilities.rst
@@ -0,0 +1,9 @@
+Utilities
+=========
+
+Local sigma clipping
+--------------------
+
+This provides a version of
+`astropy.stats.sigma_clip <https://docs.astropy.org/en/stable/api/astropy.stats.sigma_clip.html>`__
+working on a local neighbourhood.
diff --git a/sospice/__init__.py b/sospice/__init__.py
index 9f17d9c..5209c60 100644
--- a/sospice/__init__.py
+++ b/sospice/__init__.py
@@ -3,3 +3,4 @@
 from .catalog.release import Release
 from .catalog.file_metadata import FileMetadata
 from .calibrate.uncertainties import spice_error
+from .util.sigma_clipping import sigma_clip
diff --git a/sospice/util/sigma_clipping.py b/sospice/util/sigma_clipping.py
new file mode 100644
index 0000000..45fab5e
--- /dev/null
+++ b/sospice/util/sigma_clipping.py
@@ -0,0 +1,87 @@
+import numpy as np
+from scipy.ndimage import generic_filter
+from numpy import ma
+
+
+def _get_numpy_function(data, name):
+    """
+    Get NaN-aware version of numpy array function if there are any NaNs in data
+
+    Parameters
+    ----------
+    data: numpy.array
+        Data
+    name: str
+        Numpy function name, in its not NaN-aware version
+
+    Return
+    ------
+    function
+        Numpy function
+    """
+    if np.isnan(data).any():
+        name = "nan" + name
+    return getattr(np, name)
+
+
+def sigma_clip(
+    data,
+    size,
+    sigma=3,
+    sigma_lower=None,
+    sigma_upper=None,
+    maxiters=5,
+    centerfunc="median",
+    masked=True,
+):
+    """
+     Performs sigma-clipping of the input array.
+
+     Parameters
+     ----------
+    data: numpy.ndarray
+        Input array
+    size: int or tuple[int]
+        Size of the kernel used to compute the running median (or mean) and standard deviation
+    sigma: float
+        The number of standard deviations to use for both the lower and upper clipping limit.
+        This is overriden by `sigma_lower` and `sigma_upper`
+    sigma_lower: float
+        Low threshold, in units of the standard deviation of the local intensity distribution
+    sigmer_upper: float
+        High threshold, in units of the standard deviation of the local intensity distribution
+    maxiters: int
+        Maximum number of iterations to perform
+    centerfunc: str
+        Method used to estimate the center of the local intensity distribution ("median" (default) or "mean")
+    masked: bool
+        Return a `numpy.ma.MaskedArray` (default) instead of an `numpy.array`
+
+    Returns
+    -------
+    numpy.ndarray
+        Filtered array, with clipped pixels replaced by the estimated value of the center of the
+        local intensity distribution (either median or mean).
+    """
+    output = np.copy(data)
+    if type(size) is int:
+        size = (size,) * data.ndim
+    sigma_lower = sigma_lower or sigma
+    sigma_upper = sigma_upper or sigma
+    maxiters = maxiters or np.inf
+    nchanged = 1
+    iteration = 0
+    while nchanged != 0 and (iteration < maxiters):
+        iteration += 1
+        center = generic_filter(output, _get_numpy_function(output, centerfunc), size)
+        stddev = generic_filter(output, _get_numpy_function(output, "std"), size)
+        diff = output - center
+        new_mask = (diff > sigma_upper * stddev) | (diff < -sigma_lower * stddev)
+        output[new_mask] = np.nan
+        nchanged = np.count_nonzero(new_mask)
+    nan = np.isnan(output)
+    output[nan] = center[nan]
+    if masked:
+        return ma.masked_array(output, mask=nan)
+    else:
+        return output


Test patch:
diff --git a/sospice/util/tests/test_sigma_clipping.py b/sospice/util/tests/test_sigma_clipping.py
new file mode 100644
index 0000000..d00e443
--- /dev/null
+++ b/sospice/util/tests/test_sigma_clipping.py
@@ -0,0 +1,42 @@
+import pytest
+import numpy as np
+
+import astropy.units as u
+
+from ..sigma_clipping import _get_numpy_function, sigma_clip
+
+
+@pytest.fixture
+def image_1d():
+    return np.array([2, 3, 2, 4, 25, 5, 1, 2, 3, 2], dtype=float)
+
+
+@pytest.fixture
+def image_2d(image_1d):
+    image_2d = np.ones((10, 10))
+    image_2d[3, 4] = 100
+    return image_2d
+
+
+class TestSigmaClip:
+    def test_get_numpy_function(self, image_1d):
+        assert _get_numpy_function(image_1d, "mean") == np.mean
+        image_1d[0] = np.nan
+        assert _get_numpy_function(image_1d, "mean") == np.nanmean
+
+    def test_sigma_clip_1d(self, image_1d):
+        result = sigma_clip(image_1d, size=3, sigma=2, maxiters=1)
+        expected_values = np.array([2, 2, 2, 4, 5, 5, 1, 2, 2, 2], dtype=float)
+        expected_mask = np.array([0, 1, 0, 0, 1, 0, 0, 0, 1, 0], dtype=bool)
+        assert u.isclose(result.data, expected_values).all()
+        assert (result.mask == expected_mask).all()
+
+    def test_sigma_clip_2d(self, image_2d):
+        result = sigma_clip(image_2d, size=5, sigma=2, maxiters=1)
+        expected_values = np.ones_like(image_2d)
+        expected_mask = np.zeros_like(image_2d, dtype=bool)
+        expected_mask[3, 4] = True
+        assert u.isclose(result.data, expected_values).all()
+        assert (result.mask == expected_mask).all()
+        result = sigma_clip(image_2d, size=5, sigma=2, maxiters=1, masked=False)
+        assert u.isclose(result, expected_values).all()