Added named_arrays.despike() function, a wrapper around astroscrappy.detect_cosmics().

docs/conf.py

@@ -115,6 +115,7 @@
     'scipy': ('https://docs.scipy.org/doc/scipy/', None),
     'matplotlib': ('https://matplotlib.org/stable', None),
     'astropy': ('https://docs.astropy.org/en/stable/', None),
+    'astroscrappy': ('https://astroscrappy.readthedocs.io/en/stable/', None),
     'xarray': ('https://docs.xarray.dev/en/stable/', None),
     'ndfilters': ('https://ndfilters.readthedocs.io/en/stable', None),
     'colorsynth': ('https://colorsynth.readthedocs.io/en/stable', None),

docs/refs.bib

@@ -19,3 +19,20 @@ @article{Goh2017
   url = {http://distill.pub/2017/momentum},
   doi = {10.23915/distill.00006}
 }
+@ARTICLE{vanDokkum2001,
+       author = {{van Dokkum}, Pieter G.},
+        title = "{Cosmic-Ray Rejection by Laplacian Edge Detection}",
+      journal = {\pasp},
+     keywords = {Instrumentation: Detectors, Methods: Data Analysis-techniques: image processing, Astrophysics},
+         year = 2001,
+        month = nov,
+       volume = {113},
+       number = {789},
+        pages = {1420-1427},
+          doi = {10.1086/323894},
+archivePrefix = {arXiv},
+       eprint = {astro-ph/0108003},
+ primaryClass = {astro-ph},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2001PASP..113.1420V},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}

named_arrays/_functions/function_named_array_functions.py

@@ -1,4 +1,4 @@
-from typing import Callable
+from typing import Callable, Literal
 import numpy as np
 import matplotlib
 import astropy.units as u
@@ -222,3 +222,66 @@
         wavelength_max=wavelength_max,
         wavelength_norm=wavelength_norm,
     )
+
+
+@_implements(na.despike)
+def despike(
+    array: na.AbstractScalar | na.AbstractFunctionArray,
+    axis: tuple[str, str],
+    where: None | bool | na.AbstractScalar | na.AbstractFunctionArray,
+    inbkg: None | na.AbstractScalar | na.AbstractFunctionArray,
+    invar: None | float | na.AbstractScalar | na.AbstractFunctionArray,
+    sigclip: float,
+    sigfrac: float,
+    objlim: float,
+    gain: float,
+    readnoise: float,
+    satlevel: float,
+    niter: int,
+    sepmed: bool,
+    cleantype: Literal["median", "medmask", "meanmask", "idw"],
+    fsmode: Literal["median", "convolve"],
+    psfmodel: Literal["gauss", "gaussx", "gaussy", "moffat"],
+    psffwhm: float,
+    psfsize: int,
+    psfk: None | na.AbstractScalar,
+    psfbeta: float,
+    verbose: bool,
+) -> na.ScalarArray:
+
+    result = array.copy_shallow()
+
+    if isinstance(array, na.AbstractFunctionArray):
+        array = array.outputs
+    if isinstance(where, na.AbstractFunctionArray):
+        where = where.outputs
+    if isinstance(inbkg, na.AbstractFunctionArray):
+        inbkg = inbkg.outputs
+    if isinstance(invar, na.AbstractFunctionArray):
+        invar = invar.outputs
+
+    result.outputs = na.despike(
+        array=array,
+        axis=axis,
+        where=where,
+        inbkg=inbkg,
+        invar=invar,
+        sigclip=sigclip,
+        sigfrac=sigfrac,
+        objlim=objlim,
+        gain=gain,
+        readnoise=readnoise,
+        satlevel=satlevel,
+        niter=niter,
+        sepmed=sepmed,
+        cleantype=cleantype,
+        fsmode=fsmode,
+        psfmodel=psfmodel,
+        psffwhm=psffwhm,
+        psfsize=psfsize,
+        psfk=psfk,
+        psfbeta=psfbeta,
+        verbose=verbose,
+    )
+
+    return result

named_arrays/_named_array_functions.py

@@ -1,5 +1,5 @@
 from __future__ import annotations
-from typing import Sequence, overload, Type, Any, Callable, TypeVar
+from typing import Sequence, overload, Type, Any, Callable, TypeVar, Literal
 import functools
 import numpy as np
 import astropy.units as u
@@ -25,6 +25,7 @@
     "interp",
     "histogram2d",
     'jacobian',
+    'despike',
 ]
 
 ArrayT = TypeVar("ArrayT")
@@ -854,3 +855,151 @@ def jacobian(
         dx=dx,
         like=f_x,
     )
+
+
+def despike(
+    array: ArrayT,
+    axis: tuple[str, str],
+    where: None | bool | na.AbstractArray = None,
+    inbkg: None | na.AbstractArray = None,
+    invar: None | float | ArrayT = None,
+    sigclip: float = 4.5,
+    sigfrac: float = 0.3,
+    objlim: float = 5.0,
+    gain: float = 1.0,
+    readnoise: float = 6.5,
+    satlevel: float = 65536.0,
+    niter: int = 4,
+    sepmed: bool = True,
+    cleantype: Literal["median", "medmask", "meanmask", "idw"] = "meanmask",
+    fsmode: Literal["median", "convolve"] = "median",
+    psfmodel: Literal["gauss", "gaussx", "gaussy", "moffat"] = "gauss",
+    psffwhm: float = 2.5,
+    psfsize: int = 7,
+    psfk: None | na.AbstractArray = None,
+    psfbeta: float = 4.765,
+    verbose: bool = False,
+) -> ArrayT:
+    """
+    A thin wrapper around :func:`astroscrappy.detect_cosmics`
+    :cite:t:`vanDokkum2001`, which removes cosmic ray spikes from a series of
+    images.
+
+    Parameters
+    ----------
+    array
+        Input data array that will be used for cosmic ray detection. This
+        should include the sky background (or a mean background level, added
+        back in after sky subtraction), so that noise can be estimated
+        correctly from the data values. This should be in units of "counts".
+    axis
+        The two axes defining the logical axes of each image.
+    where
+        A boolean array of which pixels to consider during the cleaning
+        process. The inverse of `inmask` used in
+        :func:`astroscrappy.detect_cosmics`.
+    inbkg
+        A pre-determined background image, to be subtracted from `array`
+        before running the main detection algorithm.
+        This is used primarily with spectroscopic data, to remove
+        sky lines and the cross-section of an object continuum during
+        iteration, "protecting" them from spurious rejection (see the above
+        paper). This background is not removed from the final, cleaned output
+        (`cleanarr`). This should be in units of "counts", the same units of `array`.
+        This inbkg should be free from cosmic rays. When estimating the cosmic-ray
+        free noise of the image, we will treat ``inbkg`` as a constant Poisson
+        contribution to the variance.
+    invar
+        A pre-determined estimate of the data variance (ie. noise squared) in
+        each pixel, generated by previous processing of `array`. If provided,
+        this is used in place of an internal noise model based on `array`,
+        ``gain`` and ``readnoise``. This still gets median filtered and cleaned
+        internally, to estimate what the noise in each pixel *would* be in the
+        absence of cosmic rays. This should be in units of "counts" squared.
+    sigclip
+        Laplacian-to-noise limit for cosmic ray detection. Lower values will
+        flag more pixels as cosmic rays. Default: 4.5.
+    sigfrac
+        Fractional detection limit for neighboring pixels. For cosmic ray
+        neighbor pixels, a lapacian-to-noise detection limit of
+        sigfrac * sigclip will be used. Default: 0.3.
+    objlim
+        Minimum contrast between Laplacian image and the fine structure image.
+        Increase this value if cores of bright stars are flagged as cosmic
+        rays. Default: 5.0.
+    gain
+        Gain of the image (electrons / ADU). We always need to work in
+        electrons for cosmic ray detection. Default: 1.0
+    readnoise
+        Read noise of the image (electrons). Used to generate the noise model
+        of the image. Default: 6.5.
+    satlevel
+        Saturation of level of the image (electrons). This value is used to
+        detect saturated stars and pixels at or above this level are added to
+        the mask. Default: 65536.0.
+    niter
+        Number of iterations of the LA Cosmic algorithm to perform. Default: 4.
+    sepmed
+        Use the separable median filter instead of the full median filter.
+        The separable median is not identical to the full median filter, but
+        they are approximately the same and the separable median filter is
+        significantly faster and still detects cosmic rays well. Default: True
+    cleantype
+        Set which clean algorithm is used:\n
+        'median': An umasked 5x5 median filter\n
+        'medmask': A masked 5x5 median filter\n
+        'meanmask': A masked 5x5 mean filter\n
+        'idw': A masked 5x5 inverse distance weighted interpolation\n
+        Default: "meanmask".
+    fsmode
+        Method to build the fine structure image:\n
+        'median': Use the median filter in the standard LA Cosmic algorithm
+        'convolve': Convolve the image with the psf kernel to calculate the
+        fine structure image using a matched filter technique.
+        Default: 'median'.
+    psfmodel
+        Model to use to generate the psf kernel if fsmode == 'convolve' and
+        psfk is None. The current choices are Gaussian and Moffat profiles.
+        'gauss' and 'moffat' produce circular PSF kernels. The 'gaussx' and
+        'gaussy' produce Gaussian kernels in the x and y directions
+        respectively. Default: "gauss".
+    psffwhm
+        Full Width Half Maximum of the PSF to use to generate the kernel.
+        Default: 2.5.
+    psfsize
+        Size of the kernel to calculate. Returned kernel will have size
+        psfsize x psfsize. psfsize should be odd. Default: 7.
+    psfk
+        PSF kernel array to use for the fine structure image if
+        fsmode == 'convolve'. If None and fsmode == 'convolve', we calculate
+        the psf kernel using 'psfmodel'. Default: None.
+    psfbeta
+        Moffat beta parameter. Only used if fsmode=='convolve' and
+        psfmodel=='moffat'. Default: 4.765.
+    verbose
+        Print to the screen or not. Default: False.
+    """
+    return na._named_array_function(
+        despike,
+        array=array,
+        axis=axis,
+        where=where,
+        inbkg=inbkg,
+        invar=invar,
+        sigclip=sigclip,
+        sigfrac=sigfrac,
+        objlim=objlim,
+        gain=gain,
+        readnoise=readnoise,
+        satlevel=satlevel,
+        niter=niter,
+        sepmed=sepmed,
+        cleantype=cleantype,
+        fsmode=fsmode,
+        psfmodel=psfmodel,
+        psffwhm=psffwhm,
+        psfsize=psfsize,
+        psfk=psfk,
+        psfbeta=psfbeta,
+        verbose=verbose,
+    )

named_arrays/_scalars/scalar_named_array_functions.py

@@ -1,11 +1,12 @@
-from typing import Callable, Sequence, Any
+from typing import Callable, Sequence, Any, Literal
 import numpy as np
 import numpy.typing as npt
 import matplotlib.axes
 import matplotlib.artist
 import matplotlib.pyplot as plt
 import matplotlib.animation
 import astropy.units as u
+import astroscrappy
 import ndfilters
 import colorsynth
 import named_arrays as na
@@ -1300,3 +1301,86 @@ def ndfilter(
         ),
         axes=axes,
     )
+
+
+@_implements(na.despike)
+def despike(
+    array: na.AbstractScalarArray,
+    axis: tuple[str, str],
+    where: None | bool | na.AbstractScalarArray,
+    inbkg: None | na.AbstractScalarArray,
+    invar: None | float | na.AbstractScalarArray,
+    sigclip: float,
+    sigfrac: float,
+    objlim: float,
+    gain: float,
+    readnoise: float,
+    satlevel: float,
+    niter: int,
+    sepmed: bool,
+    cleantype: Literal["median", "medmask", "meanmask", "idw"],
+    fsmode: Literal["median", "convolve"],
+    psfmodel: Literal["gauss", "gaussx", "gaussy", "moffat"],
+    psffwhm: float,
+    psfsize: int,
+    psfk: None | na.AbstractScalarArray,
+    psfbeta: float,
+    verbose: bool,
+) -> na.ScalarArray:
+
+    try:
+        array = scalars._normalize(array)
+        where = scalars._normalize(where) if where is not None else where
+        inbkg = scalars._normalize(inbkg) if inbkg is not None else inbkg
+        invar = scalars._normalize(invar) if invar is not None else invar
+        psfk = scalars._normalize(psfk) if psfk is not None else psfk
+    except scalars.ScalarTypeError:  # pragma: nocover
+        return NotImplemented
+
+    shape = na.shape_broadcasted(
+        array,
+        where,
+        inbkg,
+        invar,
+    )
+
+    array = array.broadcast_to(shape)
+    where = where.broadcast_to(shape) if where is not None else where
+    inbkg = inbkg.broadcast_to(shape) if inbkg is not None else inbkg
+    invar = invar.broadcast_to(shape) if invar is not None else invar
+
+    if psfk is not None:
+        shape_orthogonal = {ax: shape[ax] for ax in shape if ax not in axis}
+        shape_psfk = na.broadcast_shapes(shape_orthogonal, psfk.shape)
+        psfk = na.broadcast_to(psfk, shape_psfk)
+
+    result = array.copy()
+    inmask = ~where if where is not None else where
+
+    for index in na.ndindex(shape, axis_ignored=axis):
+        result_ndarray = astroscrappy.detect_cosmics(
+            indat=array[index].ndarray,
+            inmask=inmask[index].ndarray if inmask is not None else inmask,
+            inbkg=inbkg[index].ndarray if inbkg is not None else inbkg,
+            invar=invar[index].ndarray if invar is not None else invar,
+            sigclip=sigclip,
+            sigfrac=sigfrac,
+            objlim=objlim,
+            gain=gain,
+            readnoise=readnoise,
+            satlevel=satlevel,
+            niter=niter,
+            sepmed=sepmed,
+            cleantype=cleantype,
+            fsmode=fsmode,
+            psfmodel=psfmodel,
+            psffwhm=psffwhm,
+            psfsize=psfsize,
+            psfk=psfk[index].ndarray if psfk is not None else psfk,
+            psfbeta=psfbeta,
+            verbose=verbose,
+        )
+
+        result[index].value.ndarray[:] = result_ndarray[1]
+
+    return result