GreenBankObservatory / dysh

https://dysh.readthedocs.io
Other
9 stars 3 forks source link

Masks not handled in smoothing #415

Open mpound opened 1 month ago

mpound commented 1 month ago

Describe the bug With the new implementation of flagging PR#412, smoothing needs to be updated to handle masks.

How to Reproduce

from dysh.fits import GBTFITSLoad
from dysh.util import get_project_testdata
filename = get_project_testdata() / "AGBT05B_047_01/AGBT05B_047_01.raw.acs/AGBT05B_047_01.raw.acs.fits"
sdfits = GBTFITSLoad(filename)
sdfits.summary()
sdfits.flag_channel([[170,200],[2880,2980],[31000,32768]])
scan_block = sdfits.getps(ifnum=0, plnum=0)
ta = scan_block.timeaverage()
ta.plot(xaxis_unit="chan", yaxis_unit="mK", ymin=100, ymax=600, grid=True)
ta.baseline(model="chebyshev", degree=2, exclude=[(14000,18000)], remove=True)
ts = ta.smooth('gaussian',16)
ta.plot(xaxis_unit="chan", yaxis_unit="mK", ymin=100, ymax=600, grid=True)

The smooth spectrum does not mask the flagged channels.

Partial solution In spectrum.py smooth:

        # in core.smooth, we fill masked values with np.nan.
        # astropy.convolve does not return a new mask, so we recreate
        # a decimated mask where values are nan
        mask = np.full(s1.shape, False)
        mask[np.where(s1 == np.nan)] = True
        new_data = Masked(s1 * self.flux.unit, mask)

In spectrum.core.smooth:

new_data = convolve(data, kernel, boundary="extend",  nan_treatment="fill", fill_value=np.nan, mask=mask)

This fixes spectrum.smooth, but then getps(smoothref=...) fails pytest.

Environment

teuben commented 1 month ago

in your example, channel 32768 was quoted. Trying to be an ignorant user, is that a typo, or do channels run from 0..32767 but we use a python range notation where the right edge doesn't count? can be so confusing.

Plus I noted that if the right edge is picked at 33000, there was no warning. I'd probably like to see a warning, this also tells the user the answer to my question about 32768

mpound commented 1 month ago

The right edge is inclusive. This is in the docs for flag_channel. However, you are right that no warning is issued for channels out of range. It just silently truncates at the spectral length. This is a feature of numpy, apparently

import numpy as np
x = np.arange(10)
x[0:100] = 99
print(x)
[99 99 99 99 99 99 99 99 99 99]

Should we warn?