diff --git a/changelog/98.breaking.rst b/changelog/98.breaking.rst new file mode 100644 index 0000000..8b6c8a1 --- /dev/null +++ b/changelog/98.breaking.rst @@ -0,0 +1 @@ +Update STIX spectrogram and srm loading and plotting functions for the new STIX file format. diff --git a/sunxspex/sunxspex_fitting/instruments.py b/sunxspex/sunxspex_fitting/instruments.py index 56a7a67..a5652e1 100644 --- a/sunxspex/sunxspex_fitting/instruments.py +++ b/sunxspex/sunxspex_fitting/instruments.py @@ -1110,7 +1110,7 @@ def _atimes2mdates(self, astrotimes): List of matplotlib dates. """ # convert astro time to datetime then use list comprehension to convert to matplotlib dates - return [mdates.date2num(dt.tt.datetime) for dt in astrotimes] + return [mdates.date2num(dt.utc.datetime) for dt in astrotimes] def _mdates_minute_locator(self, _obs_dt=None): """ Try to determine a nice tick separation for time axis on the lightcurve. @@ -1290,12 +1290,12 @@ def lightcurve(self, energy_ranges=None, axes=None, rebin_time=1): _y_pos = ax.get_ylim()[0] + (ax.get_ylim()[1]-ax.get_ylim()[0])*0.95 # stop region label overlapping axis spine if hasattr(self, "_start_background_time") and (type(self._start_background_time) != type(None)) and hasattr(self, "_end_background_time") and (type(self._end_background_time) != type(None)): ax.axvspan(*self._atimes2mdates([self._start_background_time, self._end_background_time]), alpha=0.1, color='orange') - ax.annotate("BG", (self._atimes2mdates([self._start_background_time])[0], _y_pos), color='orange', va="top", size=_def_fs-8) + ax.annotate("BG", (self._atimes2mdates([self._start_background_time])[0], _y_pos), color='orange', va="top", size=_def_fs-2) # plot event time range if hasattr(self, "_start_event_time") and hasattr(self, "_end_event_time"): ax.axvspan(*self._atimes2mdates([self._start_event_time, self._end_event_time]), alpha=0.1, color='purple') - ax.annotate("Evt", (self._atimes2mdates([self._start_event_time])[0], _y_pos), color='purple', va="top", size=_def_fs-8) + ax.annotate("Evt", (self._atimes2mdates([self._start_event_time])[0], _y_pos), color='purple', va="top", size=_def_fs-2) self._lightcurve_data = {"mdtimes": _ts, "lightcurves": _lcs, "lightcurve_error": _lcs_err, "energy_ranges": energy_ranges} @@ -1395,16 +1395,19 @@ def spectrogram(self, axes=None, rebin_time=1, rebin_energy=1, **kwargs): ax.set_title(self._instrument()+"Spectrogram [Counts s$^{-1}$]") + # change event and background start and end times from astropy dates to matplotlib dates + start_evt_time, end_evt_time, start_bg_time, end_bg_time = self._atimes2mdates([self._start_event_time, self._end_event_time, self._start_background_time, self._end_background_time]) + # plot background time range if there is one _y_pos = ax.get_ylim()[0] + (ax.get_ylim()[1]-ax.get_ylim()[0])*0.95 # stop region label overlapping axis spine if hasattr(self, "_start_background_time") and (type(self._start_background_time) != type(None)) and hasattr(self, "_end_background_time") and (type(self._end_background_time) != type(None)): - ax.plot(self._atimes2mdates([self._start_background_time, self._end_background_time]), [etop, etop], alpha=0.9, color='orange', lw=10) - ax.annotate("BG", (self._atimes2mdates([self._start_background_time])[0], _y_pos), color='orange', va="top", size=_def_fs-8) + ax.hlines(y=etop, xmin=start_bg_time, xmax=end_bg_time, alpha=0.9, color='orange', capstyle='butt', lw=10) + ax.annotate("BG", (start_bg_time, _y_pos), color='orange', va="top", size=_def_fs-2) # plot event time range if hasattr(self, "_start_event_time") and hasattr(self, "_end_event_time"): - ax.plot(self._atimes2mdates([self._start_event_time, self._end_event_time]), [etop, etop], alpha=0.9, color='#F37AFF', lw=10) - ax.annotate("Evt", (self._atimes2mdates([self._start_event_time])[0], _y_pos), color='#F37AFF', va="top", size=_def_fs-8) + ax.hlines(y=etop, xmin=start_evt_time, xmax=end_evt_time, alpha=0.9, color='#F37AFF', capstyle='butt', lw=10) + ax.annotate("Evt", (start_evt_time, _y_pos), color='#F37AFF', va="top", size=_def_fs-2) self._spectrogram_data = {"spectrogram": _spect, "extent": _ext} diff --git a/sunxspex/sunxspex_fitting/io.py b/sunxspex/sunxspex_fitting/io.py index 83b4801..029a4a3 100644 --- a/sunxspex/sunxspex_fitting/io.py +++ b/sunxspex/sunxspex_fitting/io.py @@ -1,9 +1,9 @@ """ The ``io`` module contains code to read instrument specific spectral data. """ -from astropy.io import fits +import numpy as np -from sunpy.io.special.genx import read_genx +from astropy.io import fits __all__ = ["_read_pha", "_read_arf", "_read_rmf", "_read_rhessi_spec_file", "_read_rhessi_srm_file", "_read_stix_spec_file", "_read_stix_srm_file"] @@ -150,6 +150,13 @@ def _read_stix_srm_file(srm_file): `dict` STIX SRM data (photon bins, count bins, and SRM in units of [counts/keV/photons]). """ - contents = read_genx(srm_file) - return {"photon_energy_bin_edges": contents["DRM"]['E_2D'], "count_energy_bin_edges": contents["DRM"]['EDGES_OUT'], - "drm": contents['DRM']['SMATRIX']} + with fits.open(srm_file) as hdul: + d0 = hdul[1].header + d1 = hdul[1].data + d3 = hdul[2].data + + pcb = np.concatenate((d1['ENERG_LO'][:, None], d1['ENERG_HI'][:, None]), axis=1) + + return {"photon_energy_bin_edges": pcb, + "count_energy_bin_edges": np.concatenate((d3['E_MIN'][:, None], d3['E_MAX'][:, None]), axis=1), + "drm": d1['MATRIX']*d0["GEOAREA"]} diff --git a/sunxspex/sunxspex_fitting/stix_spec_code.py b/sunxspex/sunxspex_fitting/stix_spec_code.py index dd51124..c76d02f 100644 --- a/sunxspex/sunxspex_fitting/stix_spec_code.py +++ b/sunxspex/sunxspex_fitting/stix_spec_code.py @@ -40,8 +40,8 @@ def _get_spec_file_info(spec_file): _plus_half_bin_width = np.ceil(time_deltas/2) t_hi = times_mids + _plus_half_bin_width - spec_stimes = [Time(sdict["0"][0]["DATE_BEG"], format='isot', scale='utc')+TimeDelta(time_diff_so2e * u.s)+TimeDelta(dt * u.ds) for dt in t_lo] - spec_etimes = [Time(sdict["0"][0]["DATE_BEG"], format='isot', scale='utc')+TimeDelta(time_diff_so2e * u.s)+TimeDelta(dt * u.ds) for dt in t_hi] + spec_stimes = [Time(sdict["0"][0]["DATE-BEG"], format='isot', scale='utc')+TimeDelta(time_diff_so2e * u.s)+TimeDelta(dt * u.cs) for dt in t_lo] + spec_etimes = [Time(sdict["0"][0]["DATE-BEG"], format='isot', scale='utc')+TimeDelta(time_diff_so2e * u.s)+TimeDelta(dt * u.cs) for dt in t_hi] time_bins = np.concatenate((np.array(spec_stimes)[:, None], np.array(spec_etimes)[:, None]), axis=1) channel_bins_inds, channel_bins = _return_masked_bins(sdict) @@ -54,21 +54,6 @@ def _get_spec_file_info(spec_file): return channel_bins, channel_bins_inds, time_bins, lvt, counts, counts_err, cts_rates, cts_rate_err -def _ds_times2s_times(ds_times): - """ STIX time might be in ds (deci-seconds). Convert to seconds. - - Parameters - ---------- - ds_times : time bins - A 2D array of time bins. - - Returns - ------- - A 2d array of the time bin edges. - """ - return ds_times - - def _return_masked_bins(sdict): """ Return the energy bins where there is data. @@ -85,9 +70,9 @@ def _return_masked_bins(sdict): e_bins = np.concatenate((sdict["4"][1]['e_low'][:, None], sdict["4"][1]['e_high'][:, None]), axis=1) # get all indices of the energy bins needed - mask_inds = sdict["1"][1]['energy_bin_mask'][0].astype(bool) + mask_inds = sdict["1"][1]['energy_bin_edge_mask'][0].astype(bool) - return mask_inds, e_bins[mask_inds] + return mask_inds, e_bins def _spec_file_units_check(stix_dict, time_dels): @@ -110,7 +95,7 @@ def _spec_file_units_check(stix_dict, time_dels): """ # stix can be saved out with counts, counts/sec, or counts/sec/cm^2/keV using counts, rate, or flux, respectively if stix_dict["0"][0]["BUNIT"] == "counts": - counts = stix_dict["2"][1]["counts"][:, 1:] + counts = stix_dict["2"][1]["counts"][:, :] counts_err = np.sqrt(counts) # should this be added in quadrature to the estimated count compression error cts_rates = counts / time_dels[:, None] cts_rate_err = counts_err / time_dels[:, None] diff --git a/sunxspex/sunxspex_fitting/tests/test_io.py b/sunxspex/sunxspex_fitting/tests/test_io.py index a7bf87d..04d65b5 100644 --- a/sunxspex/sunxspex_fitting/tests/test_io.py +++ b/sunxspex/sunxspex_fitting/tests/test_io.py @@ -104,15 +104,26 @@ def test_read_stix_spec_file(mock_open): assert np.array_equal(v[1][k], np.arange(k)) -@patch('sunxspex.sunxspex_fitting.io.read_genx') -def test_read_stix_srm_file(mock_read_genx): +@patch('astropy.io.fits.open') +def test_read_stix_srm_file(mock_open): + photon_bins = np.array([[1, 2], [2, 3]]) - count_bins = np.array([[2, 3], [4, 5]]) - drm = np.eye(2) - ret = {"DRM": {'E_2D': photon_bins, 'EDGES_OUT': count_bins, 'SMATRIX': drm}} + count_bins = np.array([[2, 3], [3, 4]]) + drm = 2*np.eye(2) + + hdul = [] + headers = [{}, {'GEOAREA': 2}, {}] + data = [{}, {'MATRIX': np.eye(2), 'ENERG_LO': np.array([1, 2]), 'ENERG_HI': np.array([2, 3])}, {'E_MIN': np.array([2, 3]), 'E_MAX':np.array([3, 4])}] + + for i in range(len(headers)): + m = MagicMock() + m.header = headers[i] + m.data = data[i] + hdul.append(m) + + mock_open.return_value.__enter__.return_value = hdul + res = _read_stix_srm_file('test.fits') - mock_read_genx.return_value = ret - res = _read_stix_srm_file('test.genx') assert np.array_equal(res['photon_energy_bin_edges'], photon_bins) assert np.array_equal(res['count_energy_bin_edges'], count_bins) assert np.array_equal(res['drm'], drm)