diff --git a/src/dysh/spectra/core.py b/src/dysh/spectra/core.py index bd6ed2d5..a35b23cc 100644 --- a/src/dysh/spectra/core.py +++ b/src/dysh/spectra/core.py @@ -4,6 +4,7 @@ import warnings from copy import deepcopy +from functools import reduce import astropy.units as u import numpy as np @@ -109,12 +110,54 @@ def find_blanks(data): Returns ------- blanks : `~numpy.ndarray` - Array with indices of \blanked integrations. + Array with indices of blanked integrations. """ return np.where(integration_isnan(data)) +def find_nonblank_ints(cycle1, cycle2, cycle3=None, cycle4=None): + """ + Find the indices of integrations that are not blanked. + + Parameters + ---------- + cycle1 : `~numpy.ndarray` + Data for cycle 1. For example, signal with the noise diode off. + cycle2 : `~numpy.ndarray` + Data for cycle 2. For example, reference with the noise diode off. + cycle3 : `~numpy.ndarray` + Data for cycle 3. For example, signal with the noise diode on. + Default is `None`. + cycle4 : `~numpy.ndarray` + Data for cycle 4. For example, reference with the noise diode on. + Default is `None`. + + Returns + ------- + goodrows : `~numpy.array` + Indices of the non-blanked rows. + """ + + nb1 = find_non_blanks(cycle1) + nb2 = find_non_blanks(cycle2) + if cycle3 is not None: + nb3 = find_non_blanks(cycle3) + else: + nb3 = nb1 + if cycle4 is not None: + nb4 = find_non_blanks(cycle4) + else: + nb4 = nb2 + goodrows = reduce(np.intersect1d, (nb1, nb2, nb3, nb4)) + + if len(goodrows) != len(cycle1): + nblanks = len(cycle1) - len(goodrows) + logger.info(f"Ignoring {nblanks} blanked integration(s).") + + return goodrows + + def exclude_to_region(exclude, refspec, fix_exclude=False): """Convert an exclude list to a list of ~specutuls.SpectralRegion. diff --git a/src/dysh/spectra/scan.py b/src/dysh/spectra/scan.py index a50468f7..05685afd 100644 --- a/src/dysh/spectra/scan.py +++ b/src/dysh/spectra/scan.py @@ -21,6 +21,7 @@ from .core import ( # fft_shift, average, find_non_blanks, + find_nonblank_ints, mean_tsys, sq_weighted_avg, tsys_weight, @@ -667,20 +668,15 @@ def __init__( self._refoffrows = sorted(list(set(self._calrows["OFF"]).intersection(set(self._scanrows)))) self._refcalon = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._refonrows] self._refcaloff = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._refoffrows] - # now remove blanked integrations - # seems like this should be done for all Scan classes! - # PS: yes. - nb1 = find_non_blanks(self._refcalon) - nb2 = find_non_blanks(self._refcaloff) - goodrows = np.intersect1d(nb1, nb2) - # Tell the user about blank integration(s) that will be ignored. - if len(goodrows) != len(self._refcalon): - nblanks = len(self._refcalon) - len(goodrows) - print(f"Ignoring {nblanks} blanked integration(s).") + + # Catch blank integrations. + goodrows = find_nonblank_ints(self._refcalon, self._refcaloff) self._refcalon = self._refcalon[goodrows] self._refcaloff = self._refcaloff[goodrows] self._refonrows = [self._refonrows[i] for i in goodrows] self._refoffrows = [self._refoffrows[i] for i in goodrows] + self._nrows = len(self._refonrows) + len(self._refoffrows) + self._nchan = len(self._refcalon[0]) self._calibrate = calibrate self._data = None @@ -977,11 +973,26 @@ def __init__( self._sigoffrows = sorted(list(set(self._calrows["OFF"]).intersection(set(self._scanrows["ON"])))) self._refonrows = sorted(list(set(self._calrows["ON"]).intersection(set(self._scanrows["OFF"])))) self._refoffrows = sorted(list(set(self._calrows["OFF"]).intersection(set(self._scanrows["OFF"])))) - self._sigcalon = gbtfits.rawspectra(self._bintable_index)[self._sigonrows] - self._nchan = len(self._sigcalon[0]) + self._sigcalon = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._sigonrows] self._sigcaloff = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._sigoffrows] self._refcalon = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._refonrows] self._refcaloff = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._refoffrows] + + # Catch blank integrations. + goodrows = find_nonblank_ints(self._sigcaloff, self._refcaloff, self._sigcalon, self._refcalon) + self._refcalon = self._refcalon[goodrows] + self._refcaloff = self._refcaloff[goodrows] + self._refonrows = [self._refonrows[i] for i in goodrows] + self._refoffrows = [self._refoffrows[i] for i in goodrows] + self._sigcalon = self._sigcalon[goodrows] + self._sigcaloff = self._sigcaloff[goodrows] + self._sigonrows = [self._sigonrows[i] for i in goodrows] + self._sigoffrows = [self._sigoffrows[i] for i in goodrows] + # Update number of rows after removing blanks. + nsigrows = len(self._sigonrows) + len(self._sigoffrows) + self._nrows = nsigrows + + self._nchan = len(self._sigcalon[0]) self._tsys = None self._exposure = None self._calibrated = None @@ -1232,6 +1243,21 @@ def __init__( self._sigcaloff = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._refoffrows] self._refcalon = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._sigonrows] self._refcaloff = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._sigoffrows] + + # Catch blank integrations. + goodrows = find_nonblank_ints(self._sigcaloff, self._refcaloff, self._sigcalon, self._refcalon) + self._refcalon = self._refcalon[goodrows] + self._refcaloff = self._refcaloff[goodrows] + self._refonrows = [self._refonrows[i] for i in goodrows] + self._refoffrows = [self._refoffrows[i] for i in goodrows] + self._sigcalon = self._sigcalon[goodrows] + self._sigcaloff = self._sigcaloff[goodrows] + self._sigonrows = [self._sigonrows[i] for i in goodrows] + self._sigoffrows = [self._sigoffrows[i] for i in goodrows] + # Update number of rows after removing blanks. + nsigrows = len(self._sigonrows) + len(self._sigoffrows) + self._nrows = nsigrows + self._nchan = len(self._sigcalon[0]) self._tsys = None self._exposure = None @@ -1507,6 +1533,21 @@ def __init__( self._sigcaloff = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._sigoffrows] self._refcalon = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._refonrows] self._refcaloff = gbtfits.rawspectra(self._bintable_index, setmask=apply_flags)[self._refoffrows] + + # Catch blank integrations. + goodrows = find_nonblank_ints(self._sigcaloff, self._refcaloff, self._sigcalon, self._refcalon) + self._refcalon = self._refcalon[goodrows] + self._refcaloff = self._refcaloff[goodrows] + self._refonrows = [self._refonrows[i] for i in goodrows] + self._refoffrows = [self._refoffrows[i] for i in goodrows] + self._sigcalon = self._sigcalon[goodrows] + self._sigcaloff = self._sigcaloff[goodrows] + self._sigonrows = [self._sigonrows[i] for i in goodrows] + self._sigoffrows = [self._sigoffrows[i] for i in goodrows] + # Update number of rows after removing blanks. + nsigrows = len(self._sigonrows) + len(self._sigoffrows) + self._nrows = nsigrows + self._nchan = len(self._sigcalon[0]) self._tsys = None self._exposure = None