Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Masks not handled in smoothing #415

Open
mpound opened this issue Oct 24, 2024 · 3 comments
Open

Masks not handled in smoothing #415

mpound opened this issue Oct 24, 2024 · 3 comments
Assignees
Labels
bug Something isn't working

Comments

@mpound
Copy link
Collaborator

mpound commented Oct 24, 2024

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

  • Dysh version 0.4.x
  • Python version 3.11
  • OS Ubuntu
@mpound mpound added the bug Something isn't working label Oct 24, 2024
@teuben
Copy link
Collaborator

teuben commented Oct 24, 2024

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
Copy link
Collaborator Author

mpound commented Oct 25, 2024

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?

@astrofle
Copy link
Collaborator

astrofle commented Dec 3, 2024

I'd preserve the numpy behavior. That has been useful for flagging from a starting channel all the way to the end channel without having to know where to stop (just set the end point of the range to a number of channels above the maximum allowed by the spectrometer).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

When branches are created from issues, their pull requests are automatically linked.

3 participants