From df2d3a142d704e2c351c46e263e7fc15ff15bfa8 Mon Sep 17 00:00:00 2001 From: astrofle Date: Thu, 5 Dec 2024 14:43:57 -0500 Subject: [PATCH 1/8] Fix: smooth now preserves observer property and _observer_location is no longer a property of Spectrum --- src/dysh/spectra/spectrum.py | 55 +++++++++++-------------- src/dysh/spectra/tests/test_scan.py | 6 +-- src/dysh/spectra/tests/test_spectrum.py | 8 +++- 3 files changed, 33 insertions(+), 36 deletions(-) diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index efad3927..b33ddb86 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -75,14 +75,8 @@ def __init__(self, *args, **kwargs): self._velocity_frame = self._target.frame.name else: self._velocity_frame = None - # @todo - have _observer_location attribute instead? - # and observer property returns getITRS(observer_location,obstime) - self._observer_location = kwargs.pop("observer_location", None) self._observer = kwargs.pop("observer", None) - if self._observer is not None and self._observer_location is not None: - raise Exception("You can only specify one of observer_location or observer") Spectrum1D.__init__(self, *args, **kwargs) - # super(Spectrum1D, self).__init__(*args, **kwargs) # Try making a target from meta. If it fails, don't worry about it. if False: if self._target is None: @@ -383,6 +377,7 @@ def stats(self, roll=0, qac=False): return out + @log_call_to_history def decimate(self, n): """ Decimate the `Spectrum` by n pixels. @@ -411,7 +406,7 @@ def decimate(self, n): new_meta["CRPIX1"] = 1.0 + (self.meta["CRPIX1"] - 1) / n + 0.5 * (n - 1) / n new_meta["CRVAL1"] += cell_shift - s = Spectrum.make_spectrum(new_data, meta=new_meta) + s = Spectrum.make_spectrum(new_data, meta=new_meta, observer_location="from_meta") if self._baseline_model is not None: s._baseline_model = None @@ -510,7 +505,7 @@ def smooth(self, method="hanning", width=1, decimate=0, kernel=None): new_data = s1 * self.flux.unit new_meta["FREQRES"] = np.sqrt((kwidth * self.meta["CDELT1"]) ** 2 + self.meta["FREQRES"] ** 2) - s = Spectrum.make_spectrum(new_data, meta=new_meta) + s = Spectrum.make_spectrum(new_data, meta=new_meta, observer_location="from_meta") s._baseline_model = self._baseline_model # Now decimate if needed. @@ -673,10 +668,7 @@ def observer(self): observer : `~astropy.coordinates.BaseCoordinateFrame` or derivative The coordinate frame of the observer if present. """ - if self._observer is None and self._observer_location is not None: - return SpectralCoord._validate_coordinate(self._observer_location.get_itrs(obstime=self._obstime)) - else: - return self._observer + return self._observer @property def velocity_frame(self): @@ -1058,7 +1050,7 @@ def fake_spectrum(cls, nchan=1024, seed=None, **kwargs): # @todo allow observer or observer_location. And/or sort this out in the constructor. @classmethod - def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): + def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None, observer=None): # , shift_topo=False): """Factory method to create a Spectrum object from a data and header. The the data are masked, the Spectrum mask will be set to the data mask. @@ -1078,13 +1070,15 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): This will be transformed to `~astropy.coordinates.ITRS` using the time of observation DATE-OBS or MJD-OBS in `meta`. If this parameter is given the special str value 'from_meta', then an observer_location will be created from SITELONG, SITELAT, and SITEELEV in the meta dictionary. + observer : `~astropy.coordinates.BaseCoordinateFrame` + Coordinate frame for the observer. + Will be ignored if DATE-OBS or MJD-OBS are present in `meta` and `observer_location` is not `None`. Returns ------- spectrum : `~dysh.spectra.Spectrum` The spectrum object """ - # @todo add resolution being the channel separation, unless we use the FREQRES column # @todo generic check_required method since I now have this code in two places (coordinates/core.py). # @todo requirement should be either DATE-OBS or MJD-OBS, but make_target() needs to be updated # in that case as well. @@ -1111,6 +1105,9 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): if not _required <= meta.keys(): raise ValueError(f"Header (meta) is missing one or more required keywords: {_required}") + if (observer or observer_location) is None: + raise Exception("One of `observer` or `observer_location` is required.") + # @todo WCS is expensive. # Possibly figure how to calculate spectral_axis instead. # @todo allow fix=False in WCS constructor? @@ -1144,15 +1141,12 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): # is_topo = is_topocentric(meta["CTYPE1"]) # GBT-specific to use CTYPE1 instead of VELDEF target = make_target(meta) vc = veldef_to_convention(meta["VELDEF"]) - # vf = astropy_frame_dict[decode_veldef(meta["VELDEF"])[1]] - # could be clever: - # obstime = Time(meta.get("DATE-OBS",meta.get("MJD-OBS",None))) - # if obstime is not None: blah blah - if "DATE-OBS" in meta: - obstime = Time(meta["DATE-OBS"]) - elif "MJD-OBS" in meta: - obstime = Time(meta["MJD-OBS"]) - if "DATE-OBS" in meta or "MJD-OBS" in meta: + + # Define an observer as needed. + if observer is not None: + obsitrs = observer + elif observer_location is not None and (meta.get("DATE-OBS") or meta.get("MJD-OBS")) is not None: + obstime = Time(meta.get("DATE-OBS") or meta.get("MJD-OBS")) if observer_location == "from_meta": try: observer_location = Observatory.get_earth_location( @@ -1160,12 +1154,9 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): ) except KeyError as ke: raise Exception(f"Not enough info to create observer_location: {ke}") - if observer_location is None: - obsitrs = None - else: - obsitrs = SpectralCoord._validate_coordinate( - attach_zero_velocities(observer_location.get_itrs(obstime=obstime)) - ) + obsitrs = SpectralCoord._validate_coordinate( + attach_zero_velocities(observer_location.get_itrs(obstime=obstime)) + ) else: warnings.warn( "'meta' does not contain DATE-OBS or MJD-OBS. Spectrum won't be convertible to certain coordinate" @@ -1199,7 +1190,7 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): # For some reason, Spectrum1D.spectral_axis created with WCS do not inherit # the radial velocity. In fact, they get no radial_velocity attribute at all! # This method creates a new spectral_axis with the given radial velocity. - if observer_location is None: + if observer_location is None and observer is None: s.set_radial_velocity_to(target.radial_velocity) # open return s @@ -1565,7 +1556,7 @@ def average_spectra(spectra, weights="tsys", align=False): tsyss = np.empty(nspec, dtype=float) xcoos = np.empty(nspec, dtype=float) ycoos = np.empty(nspec, dtype=float) - obs_location = spectra[0]._observer_location + observer = spectra[0].observer units = spectra[0].flux.unit for i, s in enumerate(spectra): @@ -1603,6 +1594,6 @@ def average_spectra(spectra, weights="tsys", align=False): new_meta["CRVAL2"] = xcoo new_meta["CRVAL3"] = ycoo - averaged = Spectrum.make_spectrum(Masked(data * units, data.mask), meta=new_meta, observer_location=obs_location) + averaged = Spectrum.make_spectrum(Masked(data * units, data.mask), meta=new_meta, observer=observer) return averaged diff --git a/src/dysh/spectra/tests/test_scan.py b/src/dysh/spectra/tests/test_scan.py index df2444ec..60abaa53 100644 --- a/src/dysh/spectra/tests/test_scan.py +++ b/src/dysh/spectra/tests/test_scan.py @@ -398,7 +398,7 @@ def test_getfs_with_args(self, data_dir): print("MWP: NO FOLD") fsscan = sdf.getfs(scan=20, ifnum=0, plnum=1, fdnum=0, fold=False, debug=True) ta = fsscan.timeaverage(weights="tsys") - # we're using astropy access here, and ujse sdfitsload.SDFITSLoad() in the other test + # we're using astropy access here, and use gbtfitsload.GBTFITSLoad() in the other test hdu = fits.open(gbtidl_file_nofold) table = hdu[1].data data = table["DATA"] @@ -410,9 +410,9 @@ def test_getfs_with_args(self, data_dir): print("MWP: FOLD") fsscan = sdf.getfs(scan=20, ifnum=0, plnum=1, fdnum=0, fold=True) ta = fsscan.timeaverage(weights="tsys") - # we will be using SDFITSLoad() here instead of astropy + # We will be using GBTFITSLoad() here instead of astropy. if True: - sdf2 = sdfitsload.SDFITSLoad(gbtidl_file) + sdf2 = gbtfitsload.GBTFITSLoad(gbtidl_file) sp = sdf2.getspec(1).flux.value.astype(np.float32) else: hdu = fits.open(gbtidl_file) diff --git a/src/dysh/spectra/tests/test_spectrum.py b/src/dysh/spectra/tests/test_spectrum.py index 09f24558..08aae5a7 100644 --- a/src/dysh/spectra/tests/test_spectrum.py +++ b/src/dysh/spectra/tests/test_spectrum.py @@ -259,7 +259,7 @@ def test_slice(self, mock_show, tmp_path): For units we only consider frequencies for now. """ meta_ignore = ["CRPIX1", "CRVAL1"] - spec_pars = ["_target", "_velocity_frame", "_observer", "_obstime", "_observer_location"] + spec_pars = ["_target", "_velocity_frame", "_observer", "_obstime"] s = slice(1000, 1100, 1) trimmed = self.ps0[s] assert trimmed.flux[0] == self.ps0.flux[s.start] @@ -328,6 +328,9 @@ def test_smooth(self): assert ss.meta["FREQRES"] == pytest.approx(abs(self.ps0.meta["CDELT1"]) * width) assert np.diff(ss.spectral_axis).mean().value == ss.meta["CDELT1"] assert ss._resolution == pytest.approx(1) + assert ss.velocity_frame == self.ps0.velocity_frame + assert ss.doppler_convention == self.ps0.doppler_convention + assert ss.observer.frame_attributes == self.ps0.observer.frame_attributes def test_smooth_decimate(self): """Test for smooth with `decimate!=width`.""" @@ -352,6 +355,9 @@ def test_smooth_decimate(self): assert np.sqrt(fwhm**2 - sss.meta["FREQRES"] ** 2) == pytest.approx( abs(self.ss.meta["CDELT1"]) * self.ss.meta["FWHM"], abs=abs(self.ss.meta["CDELT1"]) / 9.0 ) + assert ss.velocity_frame == self.ps0.velocity_frame + assert ss.doppler_convention == self.ps0.doppler_convention + assert ss.observer.frame_attributes == self.ps0.observer.frame_attributes def test_smooth_nodecimate(self): """Test for smooth without decimation.""" From 7a6c429c28244e635a3b62ee9bb44368daeaa001 Mon Sep 17 00:00:00 2001 From: astrofle Date: Thu, 5 Dec 2024 15:08:44 -0500 Subject: [PATCH 2/8] Add: class method for spectrum averaging for ease of use --- src/dysh/spectra/spectrum.py | 34 +++++++++++++++++++++++++ src/dysh/spectra/tests/test_spectrum.py | 2 ++ 2 files changed, 36 insertions(+) diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index b33ddb86..4c12df5b 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -1366,6 +1366,40 @@ def wav2idx(wav, wcs, spectral_axis, coo, sto): observer_location=Observatory[meta["TELESCOP"]], ) + def average(self, spectra, weights="tsys", align=False): + """ + Average this `Spectrum` with `spectra`. + The resulting `average` will have an exposure equal to the sum of the exposures, + and coordinates and system temperature equal to the weighted average of the coordinates and system temperatures. + + Parameters + ---------- + spectra : list of `Spectrum` + Spectra to be averaged. They must have the same number of channels. + No checks are done to ensure they are aligned. + weights: str + 'tsys' or None. If 'tsys' the weight will be calculated as: + + :math:`w = t_{exp} \times \delta\nu/T_{sys}^2` + + Default: 'tsys' + align : bool + If `True` align the `spectra` to itself. + This uses `Spectrum.align_to`. + + Returns + ------- + average : `Spectrum` + Averaged spectra. + """ + + if type(spectra) is not list: + spectra = [spectra] + + spectra += [self] + + return average_spectra(spectra, weights=weights, align=align) + # @todo figure how how to document write() #################################################################### diff --git a/src/dysh/spectra/tests/test_spectrum.py b/src/dysh/spectra/tests/test_spectrum.py index 08aae5a7..dfc5bc6e 100644 --- a/src/dysh/spectra/tests/test_spectrum.py +++ b/src/dysh/spectra/tests/test_spectrum.py @@ -527,6 +527,8 @@ def test_average_spectra(self): ps1_org = self.ps1._copy() avg = average_spectra((self.ps0, self.ps1)) + avg2 = self.ps0.average((self.ps1)) + compare_spectrum(avg, avg2, ignore_history=True, ignore_comments=True) avg = average_spectra((self.ps0, self.ps1), align=True) compare_spectrum(ps0_org, self.ps0, ignore_history=True, ignore_comments=True) From 564eddf68e43901ba996176b49731abee9cf49e3 Mon Sep 17 00:00:00 2001 From: astrofle Date: Wed, 4 Dec 2024 15:40:03 -0500 Subject: [PATCH 3/8] Fix: removes changes to CTYPE1 when changing the rest frame --- src/dysh/spectra/spectrum.py | 5 ++++- src/dysh/spectra/tests/test_spectrum.py | 13 +++++++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index 4c12df5b..3588b692 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -742,7 +742,10 @@ def set_frame(self, toframe): else: actualframe = astropy_frame_dict.get(toframe, toframe) self._spectral_axis = self._spectral_axis.with_observer_stationary_relative_to(actualframe) - self._meta["CTYPE1"] = change_ctype(self._meta["CTYPE1"], toframe) + # This line is commented because: + # SDFITS defines CTYPE1 as always being the TOPO frequency. + # See Issue #373 on GitHub. + # self._meta["CTYPE1"] = change_ctype(self._meta["CTYPE1"], toframe) if isinstance(actualframe, str): self._velocity_frame = actualframe else: diff --git a/src/dysh/spectra/tests/test_spectrum.py b/src/dysh/spectra/tests/test_spectrum.py index dfc5bc6e..55bb9056 100644 --- a/src/dysh/spectra/tests/test_spectrum.py +++ b/src/dysh/spectra/tests/test_spectrum.py @@ -332,6 +332,19 @@ def test_smooth(self): assert ss.doppler_convention == self.ps0.doppler_convention assert ss.observer.frame_attributes == self.ps0.observer.frame_attributes + # Now, change the reference frame and see if it still works. + from dysh.coordinates import astropy_frame_dict + + s = Spectrum.fake_spectrum() + for frame in astropy_frame_dict.keys(): + try: + s.set_frame(frame) + except KeyError: + print(f"set_frame fails for: {frame}") + continue + print(f"frame set to: {frame}") + s.smooth("box", 3) + def test_smooth_decimate(self): """Test for smooth with `decimate!=width`.""" width = 10 From 525daf15f34169d348f28865d30d1ba063deb627 Mon Sep 17 00:00:00 2001 From: astrofle Date: Fri, 10 Jan 2025 14:44:05 -0500 Subject: [PATCH 4/8] Remove check for observer or observer_location --- src/dysh/spectra/spectrum.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index 3588b692..a40a9671 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -1108,9 +1108,6 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None, observe if not _required <= meta.keys(): raise ValueError(f"Header (meta) is missing one or more required keywords: {_required}") - if (observer or observer_location) is None: - raise Exception("One of `observer` or `observer_location` is required.") - # @todo WCS is expensive. # Possibly figure how to calculate spectral_axis instead. # @todo allow fix=False in WCS constructor? From cafc14ec0656a6a6a49f050bb7c872067fcd3d88 Mon Sep 17 00:00:00 2001 From: astrofle Date: Thu, 5 Dec 2024 14:43:57 -0500 Subject: [PATCH 5/8] Fix: smooth now preserves observer property and _observer_location is no longer a property of Spectrum --- src/dysh/spectra/spectrum.py | 55 +++++++++++-------------- src/dysh/spectra/tests/test_scan.py | 6 +-- src/dysh/spectra/tests/test_spectrum.py | 8 +++- 3 files changed, 33 insertions(+), 36 deletions(-) diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index 65da7cfe..6064e96d 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -76,14 +76,8 @@ def __init__(self, *args, **kwargs): self._velocity_frame = self._target.frame.name else: self._velocity_frame = None - # @todo - have _observer_location attribute instead? - # and observer property returns getITRS(observer_location,obstime) - self._observer_location = kwargs.pop("observer_location", None) self._observer = kwargs.pop("observer", None) - if self._observer is not None and self._observer_location is not None: - raise Exception("You can only specify one of observer_location or observer") Spectrum1D.__init__(self, *args, **kwargs) - # super(Spectrum1D, self).__init__(*args, **kwargs) # Try making a target from meta. If it fails, don't worry about it. if False: if self._target is None: @@ -384,6 +378,7 @@ def stats(self, roll=0, qac=False): return out + @log_call_to_history def decimate(self, n): """ Decimate the `Spectrum` by n pixels. @@ -412,7 +407,7 @@ def decimate(self, n): new_meta["CRPIX1"] = 1.0 + (self.meta["CRPIX1"] - 1) / n + 0.5 * (n - 1) / n new_meta["CRVAL1"] += cell_shift - s = Spectrum.make_spectrum(new_data, meta=new_meta) + s = Spectrum.make_spectrum(new_data, meta=new_meta, observer_location="from_meta") if self._baseline_model is not None: s._baseline_model = None @@ -511,7 +506,7 @@ def smooth(self, method="hanning", width=1, decimate=0, kernel=None): new_data = s1 * self.flux.unit new_meta["FREQRES"] = np.sqrt((kwidth * self.meta["CDELT1"]) ** 2 + self.meta["FREQRES"] ** 2) - s = Spectrum.make_spectrum(new_data, meta=new_meta) + s = Spectrum.make_spectrum(new_data, meta=new_meta, observer_location="from_meta") s._baseline_model = self._baseline_model # Now decimate if needed. @@ -674,10 +669,7 @@ def observer(self): observer : `~astropy.coordinates.BaseCoordinateFrame` or derivative The coordinate frame of the observer if present. """ - if self._observer is None and self._observer_location is not None: - return SpectralCoord._validate_coordinate(self._observer_location.get_itrs(obstime=self._obstime)) - else: - return self._observer + return self._observer @property def velocity_frame(self): @@ -1078,7 +1070,7 @@ def fake_spectrum(cls, nchan=1024, seed=None, **kwargs): # @todo allow observer or observer_location. And/or sort this out in the constructor. @classmethod - def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): + def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None, observer=None): # , shift_topo=False): """Factory method to create a Spectrum object from a data and header. The the data are masked, the Spectrum mask will be set to the data mask. @@ -1098,13 +1090,15 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): This will be transformed to `~astropy.coordinates.ITRS` using the time of observation DATE-OBS or MJD-OBS in `meta`. If this parameter is given the special str value 'from_meta', then an observer_location will be created from SITELONG, SITELAT, and SITEELEV in the meta dictionary. + observer : `~astropy.coordinates.BaseCoordinateFrame` + Coordinate frame for the observer. + Will be ignored if DATE-OBS or MJD-OBS are present in `meta` and `observer_location` is not `None`. Returns ------- spectrum : `~dysh.spectra.Spectrum` The spectrum object """ - # @todo add resolution being the channel separation, unless we use the FREQRES column # @todo generic check_required method since I now have this code in two places (coordinates/core.py). # @todo requirement should be either DATE-OBS or MJD-OBS, but make_target() needs to be updated # in that case as well. @@ -1131,6 +1125,9 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): if not _required <= meta.keys(): raise ValueError(f"Header (meta) is missing one or more required keywords: {_required}") + if (observer or observer_location) is None: + raise Exception("One of `observer` or `observer_location` is required.") + # @todo WCS is expensive. # Possibly figure how to calculate spectral_axis instead. # @todo allow fix=False in WCS constructor? @@ -1164,15 +1161,12 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): # is_topo = is_topocentric(meta["CTYPE1"]) # GBT-specific to use CTYPE1 instead of VELDEF target = make_target(meta) vc = veldef_to_convention(meta["VELDEF"]) - # vf = astropy_frame_dict[decode_veldef(meta["VELDEF"])[1]] - # could be clever: - # obstime = Time(meta.get("DATE-OBS",meta.get("MJD-OBS",None))) - # if obstime is not None: blah blah - if "DATE-OBS" in meta: - obstime = Time(meta["DATE-OBS"]) - elif "MJD-OBS" in meta: - obstime = Time(meta["MJD-OBS"]) - if "DATE-OBS" in meta or "MJD-OBS" in meta: + + # Define an observer as needed. + if observer is not None: + obsitrs = observer + elif observer_location is not None and (meta.get("DATE-OBS") or meta.get("MJD-OBS")) is not None: + obstime = Time(meta.get("DATE-OBS") or meta.get("MJD-OBS")) if observer_location == "from_meta": try: observer_location = Observatory.get_earth_location( @@ -1180,12 +1174,9 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): ) except KeyError as ke: raise Exception(f"Not enough info to create observer_location: {ke}") - if observer_location is None: - obsitrs = None - else: - obsitrs = SpectralCoord._validate_coordinate( - attach_zero_velocities(observer_location.get_itrs(obstime=obstime)) - ) + obsitrs = SpectralCoord._validate_coordinate( + attach_zero_velocities(observer_location.get_itrs(obstime=obstime)) + ) else: warnings.warn( "'meta' does not contain DATE-OBS or MJD-OBS. Spectrum won't be convertible to certain coordinate" @@ -1219,7 +1210,7 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None): # For some reason, Spectrum1D.spectral_axis created with WCS do not inherit # the radial velocity. In fact, they get no radial_velocity attribute at all! # This method creates a new spectral_axis with the given radial velocity. - if observer_location is None: + if observer_location is None and observer is None: s.set_radial_velocity_to(target.radial_velocity) # open return s @@ -1585,7 +1576,7 @@ def average_spectra(spectra, weights="tsys", align=False): tsyss = np.empty(nspec, dtype=float) xcoos = np.empty(nspec, dtype=float) ycoos = np.empty(nspec, dtype=float) - obs_location = spectra[0]._observer_location + observer = spectra[0].observer units = spectra[0].flux.unit for i, s in enumerate(spectra): @@ -1623,6 +1614,6 @@ def average_spectra(spectra, weights="tsys", align=False): new_meta["CRVAL2"] = xcoo new_meta["CRVAL3"] = ycoo - averaged = Spectrum.make_spectrum(Masked(data * units, data.mask), meta=new_meta, observer_location=obs_location) + averaged = Spectrum.make_spectrum(Masked(data * units, data.mask), meta=new_meta, observer=observer) return averaged diff --git a/src/dysh/spectra/tests/test_scan.py b/src/dysh/spectra/tests/test_scan.py index df2444ec..60abaa53 100644 --- a/src/dysh/spectra/tests/test_scan.py +++ b/src/dysh/spectra/tests/test_scan.py @@ -398,7 +398,7 @@ def test_getfs_with_args(self, data_dir): print("MWP: NO FOLD") fsscan = sdf.getfs(scan=20, ifnum=0, plnum=1, fdnum=0, fold=False, debug=True) ta = fsscan.timeaverage(weights="tsys") - # we're using astropy access here, and ujse sdfitsload.SDFITSLoad() in the other test + # we're using astropy access here, and use gbtfitsload.GBTFITSLoad() in the other test hdu = fits.open(gbtidl_file_nofold) table = hdu[1].data data = table["DATA"] @@ -410,9 +410,9 @@ def test_getfs_with_args(self, data_dir): print("MWP: FOLD") fsscan = sdf.getfs(scan=20, ifnum=0, plnum=1, fdnum=0, fold=True) ta = fsscan.timeaverage(weights="tsys") - # we will be using SDFITSLoad() here instead of astropy + # We will be using GBTFITSLoad() here instead of astropy. if True: - sdf2 = sdfitsload.SDFITSLoad(gbtidl_file) + sdf2 = gbtfitsload.GBTFITSLoad(gbtidl_file) sp = sdf2.getspec(1).flux.value.astype(np.float32) else: hdu = fits.open(gbtidl_file) diff --git a/src/dysh/spectra/tests/test_spectrum.py b/src/dysh/spectra/tests/test_spectrum.py index a606069a..dddd2b71 100644 --- a/src/dysh/spectra/tests/test_spectrum.py +++ b/src/dysh/spectra/tests/test_spectrum.py @@ -259,7 +259,7 @@ def test_slice(self, mock_show, tmp_path): For units we only consider frequencies for now. """ meta_ignore = ["CRPIX1", "CRVAL1"] - spec_pars = ["_target", "_velocity_frame", "_observer", "_obstime", "_observer_location"] + spec_pars = ["_target", "_velocity_frame", "_observer", "_obstime"] s = slice(1000, 1100, 1) trimmed = self.ps0[s] assert trimmed.flux[0] == self.ps0.flux[s.start] @@ -328,6 +328,9 @@ def test_smooth(self): assert ss.meta["FREQRES"] == pytest.approx(abs(self.ps0.meta["CDELT1"]) * width) assert np.diff(ss.spectral_axis).mean().value == ss.meta["CDELT1"] assert ss._resolution == pytest.approx(1) + assert ss.velocity_frame == self.ps0.velocity_frame + assert ss.doppler_convention == self.ps0.doppler_convention + assert ss.observer.frame_attributes == self.ps0.observer.frame_attributes # Now, change the reference frame and see if it still works. from dysh.coordinates import astropy_frame_dict @@ -365,6 +368,9 @@ def test_smooth_decimate(self): assert np.sqrt(fwhm**2 - sss.meta["FREQRES"] ** 2) == pytest.approx( abs(self.ss.meta["CDELT1"]) * self.ss.meta["FWHM"], abs=abs(self.ss.meta["CDELT1"]) / 9.0 ) + assert ss.velocity_frame == self.ps0.velocity_frame + assert ss.doppler_convention == self.ps0.doppler_convention + assert ss.observer.frame_attributes == self.ps0.observer.frame_attributes def test_smooth_nodecimate(self): """Test for smooth without decimation.""" From dc65e262156d79474aa70a081da416dacca8536f Mon Sep 17 00:00:00 2001 From: astrofle Date: Thu, 5 Dec 2024 15:08:44 -0500 Subject: [PATCH 6/8] Add: class method for spectrum averaging for ease of use --- src/dysh/spectra/spectrum.py | 34 +++++++++++++++++++++++++ src/dysh/spectra/tests/test_spectrum.py | 2 ++ 2 files changed, 36 insertions(+) diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index 6064e96d..e638364b 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -1386,6 +1386,40 @@ def wav2idx(wav, wcs, spectral_axis, coo, sto): observer_location=Observatory[meta["TELESCOP"]], ) + def average(self, spectra, weights="tsys", align=False): + """ + Average this `Spectrum` with `spectra`. + The resulting `average` will have an exposure equal to the sum of the exposures, + and coordinates and system temperature equal to the weighted average of the coordinates and system temperatures. + + Parameters + ---------- + spectra : list of `Spectrum` + Spectra to be averaged. They must have the same number of channels. + No checks are done to ensure they are aligned. + weights: str + 'tsys' or None. If 'tsys' the weight will be calculated as: + + :math:`w = t_{exp} \times \delta\nu/T_{sys}^2` + + Default: 'tsys' + align : bool + If `True` align the `spectra` to itself. + This uses `Spectrum.align_to`. + + Returns + ------- + average : `Spectrum` + Averaged spectra. + """ + + if type(spectra) is not list: + spectra = [spectra] + + spectra += [self] + + return average_spectra(spectra, weights=weights, align=align) + # @todo figure how how to document write() #################################################################### diff --git a/src/dysh/spectra/tests/test_spectrum.py b/src/dysh/spectra/tests/test_spectrum.py index dddd2b71..dc1f5b9b 100644 --- a/src/dysh/spectra/tests/test_spectrum.py +++ b/src/dysh/spectra/tests/test_spectrum.py @@ -540,6 +540,8 @@ def test_average_spectra(self): ps1_org = self.ps1._copy() avg = average_spectra((self.ps0, self.ps1)) + avg2 = self.ps0.average((self.ps1)) + compare_spectrum(avg, avg2, ignore_history=True, ignore_comments=True) avg = average_spectra((self.ps0, self.ps1), align=True) compare_spectrum(ps0_org, self.ps0, ignore_history=True, ignore_comments=True) From 93261c001c997f4278c893229324efabf864a125 Mon Sep 17 00:00:00 2001 From: astrofle Date: Wed, 4 Dec 2024 15:40:03 -0500 Subject: [PATCH 7/8] Fix: removes changes to CTYPE1 when changing the rest frame --- src/dysh/spectra/spectrum.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index e638364b..d015672c 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -764,7 +764,9 @@ def set_frame(self, toframe): self._velocity_frame = tfl else: self._velocity_frame = tfl.name - # While it is incorrect to change CTYPE1, it is reasonable to change VELDEF + # While it is incorrect to change CTYPE1, it is reasonable to change VELDEF. + # SDFITS defines CTYPE1 as always being the TOPO frequency. + # See Issue #373 on GitHub. self.meta["VELDEF"] = change_ctype(self.meta["VELDEF"], self._velocity_frame) def with_frame(self, toframe): From d74832767c13c6513f3d0e38a1abd3d3a9c0ffc8 Mon Sep 17 00:00:00 2001 From: astrofle Date: Fri, 10 Jan 2025 14:44:05 -0500 Subject: [PATCH 8/8] Remove check for observer or observer_location --- src/dysh/spectra/spectrum.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/dysh/spectra/spectrum.py b/src/dysh/spectra/spectrum.py index d015672c..6c573d99 100644 --- a/src/dysh/spectra/spectrum.py +++ b/src/dysh/spectra/spectrum.py @@ -1127,9 +1127,6 @@ def make_spectrum(cls, data, meta, use_wcs=True, observer_location=None, observe if not _required <= meta.keys(): raise ValueError(f"Header (meta) is missing one or more required keywords: {_required}") - if (observer or observer_location) is None: - raise Exception("One of `observer` or `observer_location` is required.") - # @todo WCS is expensive. # Possibly figure how to calculate spectral_axis instead. # @todo allow fix=False in WCS constructor?