From 60c54cab26050bee87800b86bedf8b26210c946c Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 25 Oct 2024 14:54:02 -0600 Subject: [PATCH 1/2] Added `named_arrays.despike()` function, a wrapper around `astroscrappy.detect_cosmics()`. --- docs/conf.py | 1 + docs/refs.bib | 17 ++ .../function_named_array_functions.py | 65 +++++++- named_arrays/_named_array_functions.py | 151 +++++++++++++++++- .../_scalars/scalar_named_array_functions.py | 86 +++++++++- .../uncertainties_named_array_functions.py | 78 ++++++++- named_arrays/tests/test_despike.py | 48 ++++++ pyproject.toml | 1 + 8 files changed, 442 insertions(+), 5 deletions(-) create mode 100644 named_arrays/tests/test_despike.py diff --git a/docs/conf.py b/docs/conf.py index 9deb0ca4..c4c2a97f 100644 --- a/docs/conf.py +++ b/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), diff --git a/docs/refs.bib b/docs/refs.bib index 1b20a6a1..ad1c5ff4 100644 --- a/docs/refs.bib +++ b/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} +} diff --git a/named_arrays/_functions/function_named_array_functions.py b/named_arrays/_functions/function_named_array_functions.py index 354ac1df..7cef8f91 100644 --- a/named_arrays/_functions/function_named_array_functions.py +++ b/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 @@ def colorsynth_colorbar( 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 diff --git a/named_arrays/_named_array_functions.py b/named_arrays/_named_array_functions.py index 5fd5c92f..0ba5a680 100644 --- a/named_arrays/_named_array_functions.py +++ b/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, + ) diff --git a/named_arrays/_scalars/scalar_named_array_functions.py b/named_arrays/_scalars/scalar_named_array_functions.py index 356004b4..c0321211 100644 --- a/named_arrays/_scalars/scalar_named_array_functions.py +++ b/named_arrays/_scalars/scalar_named_array_functions.py @@ -1,4 +1,4 @@ -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 @@ -6,6 +6,7 @@ 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: + 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 diff --git a/named_arrays/_scalars/uncertainties/uncertainties_named_array_functions.py b/named_arrays/_scalars/uncertainties/uncertainties_named_array_functions.py index 24fdb1bd..05b39c67 100644 --- a/named_arrays/_scalars/uncertainties/uncertainties_named_array_functions.py +++ b/named_arrays/_scalars/uncertainties/uncertainties_named_array_functions.py @@ -1,4 +1,4 @@ -from typing import Callable, Sequence +from typing import Callable, Sequence, Literal import numpy as np import numpy.typing as npt import matplotlib.axes @@ -830,4 +830,78 @@ def ndfilter( where=where.distribution, **kwargs, ) - ) \ No newline at end of file + ) + + +@_implements(na.despike) +def despike( + array: na.AbstractScalar, + axis: tuple[str, str], + where: None | bool | na.AbstractScalar, + inbkg: None | na.AbstractScalar, + invar: None | float | na.AbstractScalar, + 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: + + try: + array = uncertainties._normalize(array) + where = uncertainties._normalize(where) + inbkg = uncertainties._normalize(inbkg) + invar = uncertainties._normalize(invar) + psfk = uncertainties._normalize(psfk) + except uncertainties.UncertainScalarTypeError: + return NotImplemented + + kwargs = dict( + axis=axis, + 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, + psfbeta=psfbeta, + verbose=verbose, + ) + + result = array.copy_shallow() + result.nominal = na.despike( + array=array.nominal, + where=where.nominal, + inbkg=inbkg.nominal, + invar=invar.nominal, + psfk=psfk.nominal, + **kwargs, + ) + result.distribution = na.despike( + array=array.distribution, + where=where.distribution, + inbkg=inbkg.distribution, + invar=invar.distribution, + psfk=psfk.distribution, + **kwargs, + ) + + return result diff --git a/named_arrays/tests/test_despike.py b/named_arrays/tests/test_despike.py new file mode 100644 index 00000000..b6135443 --- /dev/null +++ b/named_arrays/tests/test_despike.py @@ -0,0 +1,48 @@ +import pytest +import named_arrays as na + +shape = dict(x=101, y=101) + +img = na.random.normal(0, 6.5, shape) +spikes = 1000 * na.random.poisson(0.001, shape) + + +@pytest.mark.parametrize( + argnames="array", + argvalues=[ + img + 100 * spikes, + img + na.NormalUncertainScalarArray(100, 5) * spikes, + na.FunctionArray( + inputs=None, + outputs=img + 100 * spikes, + ), + ], +) +@pytest.mark.parametrize( + argnames="axis", + argvalues=[ + tuple(shape), + ], +) +@pytest.mark.parametrize("where", [None]) +@pytest.mark.parametrize("inbkg", [None]) +@pytest.mark.parametrize("invar", [None]) +@pytest.mark.parametrize("psfk", [None]) +def test_despike( + array: na.AbstractArray, + axis: tuple[str, str], + where: None | bool | na.AbstractArray, + inbkg: None | na.AbstractArray, + invar: None | float | na.AbstractArray, + psfk: None | na.AbstractArray, +): + result = na.despike( + array=array, + axis=axis, + where=where, + inbkg=inbkg, + invar=invar, + psfk=psfk, + ) + + assert result.sum() != array.sum() diff --git a/pyproject.toml b/pyproject.toml index f39282ce..e0562bfd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,6 +20,7 @@ dependencies = [ "matplotlib", "scipy", 'astropy', + "astroscrappy", "ndfilters==0.3.0", "colorsynth==0.1.3", ] From 7f025e9ab7c5d3825470186c087d3f740e075209 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 1 Nov 2024 15:13:34 -0600 Subject: [PATCH 2/2] coverage --- named_arrays/_scalars/scalar_named_array_functions.py | 2 +- .../uncertainties/uncertainties_named_array_functions.py | 2 +- named_arrays/tests/test_despike.py | 5 ++++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/named_arrays/_scalars/scalar_named_array_functions.py b/named_arrays/_scalars/scalar_named_array_functions.py index c0321211..87964342 100644 --- a/named_arrays/_scalars/scalar_named_array_functions.py +++ b/named_arrays/_scalars/scalar_named_array_functions.py @@ -1334,7 +1334,7 @@ def despike( 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: + except scalars.ScalarTypeError: # pragma: nocover return NotImplemented shape = na.shape_broadcasted( diff --git a/named_arrays/_scalars/uncertainties/uncertainties_named_array_functions.py b/named_arrays/_scalars/uncertainties/uncertainties_named_array_functions.py index 05b39c67..f09e5d5b 100644 --- a/named_arrays/_scalars/uncertainties/uncertainties_named_array_functions.py +++ b/named_arrays/_scalars/uncertainties/uncertainties_named_array_functions.py @@ -864,7 +864,7 @@ def despike( inbkg = uncertainties._normalize(inbkg) invar = uncertainties._normalize(invar) psfk = uncertainties._normalize(psfk) - except uncertainties.UncertainScalarTypeError: + except uncertainties.UncertainScalarTypeError: # pragma: nocover return NotImplemented kwargs = dict( diff --git a/named_arrays/tests/test_despike.py b/named_arrays/tests/test_despike.py index b6135443..50fcd454 100644 --- a/named_arrays/tests/test_despike.py +++ b/named_arrays/tests/test_despike.py @@ -27,7 +27,10 @@ @pytest.mark.parametrize("where", [None]) @pytest.mark.parametrize("inbkg", [None]) @pytest.mark.parametrize("invar", [None]) -@pytest.mark.parametrize("psfk", [None]) +@pytest.mark.parametrize("psfk", [ + None, + na.ScalarArray.ones(dict(x=5, y=5)), +]) def test_despike( array: na.AbstractArray, axis: tuple[str, str],