diff --git a/changelog/259.feature.rst b/changelog/259.feature.rst new file mode 100644 index 0000000..458c76c --- /dev/null +++ b/changelog/259.feature.rst @@ -0,0 +1 @@ +Migrated `e-CALLISTO` and `WAVES` spectrogram creation to the new `from_raw` factory pattern, which extracts parsing logic out of the main `SpectrogramFactory`. diff --git a/radiospectra/spectrogram/sources/callisto.py b/radiospectra/spectrogram/sources/callisto.py index e7a41bf..3f49948 100644 --- a/radiospectra/spectrogram/sources/callisto.py +++ b/radiospectra/spectrogram/sources/callisto.py @@ -1,6 +1,9 @@ import astropy.units as u from astropy.coordinates import EarthLocation, SkyCoord +from sunpy.net import attrs as a +from sunpy.time import parse_time + from radiospectra.spectrogram.meta import SpectrogramMeta from radiospectra.spectrogram.spectrogrambase import GenericSpectrogram @@ -60,5 +63,41 @@ def observatory_location(self): return self.meta.observer_coordinate @classmethod - def is_datasource_for(cls, data, meta, **kwargs): - return meta["instrument"] == "e-CALLISTO" or meta["detector"] == "e-CALLISTO" + def is_datasource_for(cls, data_or_header, meta_or_raw, **kwargs): + meta = data_or_header if hasattr(data_or_header, "get") else meta_or_raw + if not hasattr(meta, "get"): + return False + return "e-CALLISTO" in meta.get("CONTENT", "") or meta.get("instrument") == "e-CALLISTO" + + @classmethod + def from_raw(cls, header, raw_object): + hd_pairs = raw_object + data = hd_pairs[0].data + times = hd_pairs[1].data["TIME"].flatten() * u.s + freqs = hd_pairs[1].data["FREQUENCY"].flatten() * u.MHz + start_time = parse_time(hd_pairs[0].header["DATE-OBS"] + " " + hd_pairs[0].header["TIME-OBS"]) + try: + end_time = parse_time(hd_pairs[0].header["DATE-END"] + " " + hd_pairs[0].header["TIME-END"]) + except ValueError: + # See https://github.com/sunpy/radiospectra/issues/74 + time_comps = hd_pairs[0].header["TIME-END"].split(":") + time_comps[0] = "00" + fixed_time = ":".join(time_comps) + date_offset = parse_time(hd_pairs[0].header["DATE-END"] + " " + fixed_time) + end_time = date_offset + 1 * u.day + + times = start_time + times + meta = CALISTOMeta( + { + "fits_meta": hd_pairs[0].header, + "detector": "e-CALLISTO", + "instrument": "e-CALLISTO", + "observatory": hd_pairs[0].header["INSTRUME"], + "start_time": start_time, + "end_time": end_time, + "wavelength": a.Wavelength(freqs.min(), freqs.max()), + "times": times, + "freqs": freqs, + } + ) + return cls(data, meta) diff --git a/radiospectra/spectrogram/sources/eovsa.py b/radiospectra/spectrogram/sources/eovsa.py index ef9dd6a..64e6a5a 100644 --- a/radiospectra/spectrogram/sources/eovsa.py +++ b/radiospectra/spectrogram/sources/eovsa.py @@ -1,3 +1,9 @@ +import astropy.units as u +from astropy.time import Time + +from sunpy.net import attrs as a +from sunpy.time import parse_time + from radiospectra.spectrogram.meta import SpectrogramMeta from radiospectra.spectrogram.spectrogrambase import GenericSpectrogram @@ -44,8 +50,35 @@ def polarisation(self): return self.meta.polarisation @classmethod - def is_datasource_for(cls, data, meta, **kwargs): - return meta["instrument"] == "EOVSA" or meta["detector"] == "EOVSA" + def is_datasource_for(cls, data_or_header, meta_or_raw, **kwargs): + meta = data_or_header if hasattr(data_or_header, "get") else meta_or_raw + if not hasattr(meta, "get"): + return False + return ( + meta.get("TELESCOP", "") == "EOVSA" + or meta.get("instrument", "") == "EOVSA" + or meta.get("detector", "") == "EOVSA" + ) - # TODO fix time gaps for plots need to render them as gaps - # can prob do when generateing proper pcolormesh grid but then prob doesn't belong here + @classmethod + def from_raw(cls, header, raw_object): + hd_pairs = raw_object + times = Time(hd_pairs[2].data["mjd"] + hd_pairs[2].data["time"] / 1000.0 / 86400.0, format="mjd") + freqs = hd_pairs[1].data["sfreq"] * u.GHz + data = hd_pairs[0].data + start_time = parse_time(hd_pairs[0].header["DATE_OBS"]) + end_time = parse_time(hd_pairs[0].header["DATE_END"]) + meta = EOVSAMeta( + { + "fits_meta": hd_pairs[0].header, + "detector": "EOVSA", + "instrument": "EOVSA", + "observatory": "Owens Valley", + "start_time": start_time, + "end_time": end_time, + "wavelength": a.Wavelength(freqs.min(), freqs.max()), + "times": times, + "freqs": freqs, + } + ) + return cls(data, meta) diff --git a/radiospectra/spectrogram/sources/ilofar357.py b/radiospectra/spectrogram/sources/ilofar357.py index b277efb..7cd5236 100644 --- a/radiospectra/spectrogram/sources/ilofar357.py +++ b/radiospectra/spectrogram/sources/ilofar357.py @@ -1,5 +1,14 @@ +import numpy as np + +import astropy.units as u +from astropy.time import Time + +from sunpy import log +from sunpy.net import attrs as a + from radiospectra.spectrogram.meta import SpectrogramMeta from radiospectra.spectrogram.spectrogrambase import GenericSpectrogram +from radiospectra.utils import subband_to_freq __all__ = [ "ILOFARMode357Spectrogram", @@ -40,5 +49,52 @@ def polarisation(self): return self.meta.polarisation @classmethod - def is_datasource_for(cls, data, meta, **kwargs): - return meta["instrument"] == "ILOFAR" + def is_datasource_for(cls, data_or_header, meta_or_raw, **kwargs): + meta = data_or_header if hasattr(data_or_header, "get") else meta_or_raw + if not hasattr(meta, "get"): + return False + return meta.get("instrument") == "ILOFAR" + + @classmethod + def from_raw(cls, header, raw_object): + file = header["file_path"] + subbands = (np.arange(54, 454, 2), np.arange(54, 454, 2), np.arange(54, 230, 2)) + num_subbands = 488 + + data = np.fromfile(file) + polarisation = file.stem[-1] + + num_times = data.shape[0] / num_subbands + if not num_times.is_integer(): + log.warning("BST file seems incomplete dropping incomplete frequencies") + num_times = np.floor(num_times).astype(int) + truncate = num_times * num_subbands + data = data[:truncate] + data = data.reshape(-1, num_subbands).T # (Freq x Time).T = (Time x Freq) + dt = np.arange(num_times) * 1 * u.s + start_time = Time.strptime(file.name.split("_bst")[0], "%Y%m%d_%H%M%S") + times = start_time + dt + + obs_mode = (3, 5, 7) + freqs = [subband_to_freq(sb, mode) for sb, mode in zip(subbands, obs_mode)] + + # 1st 200 sbs mode 3, next 200 sbs mode 5, last 88 sbs mode 7 + spec = {0: data[:200, :], 1: data[200:400, :], 2: data[400:, :]} + results = [] + for i in range(3): + meta = ILOFARMeta( + { + "instrument": "ILOFAR", + "observatory": "Birr (IE613)", + "start_time": times[0], + "mode": obs_mode[i], + "wavelength": a.Wavelength(freqs[i][0], freqs[i][-1]), + "freqs": freqs[i], + "times": times, + "end_time": times[-1], + "detector": "ILOFAR", + "polarisation": polarisation, + } + ) + results.append(cls(spec[i], meta)) + return results diff --git a/radiospectra/spectrogram/sources/nda.py b/radiospectra/spectrogram/sources/nda.py index 0cdfee5..a4b928f 100644 --- a/radiospectra/spectrogram/sources/nda.py +++ b/radiospectra/spectrogram/sources/nda.py @@ -45,7 +45,38 @@ def observatory_location(self): return None @classmethod - def is_datasource_for(cls, data, meta, **kwargs): - telescope = meta.get("fits_meta", {}).get("TELESCOP", "") + def is_datasource_for(cls, data_or_header, meta_or_raw, **kwargs): + meta = data_or_header if hasattr(data_or_header, "get") else meta_or_raw + if not hasattr(meta, "get"): + return False + telescope = meta.get("fits_meta", {}).get("TELESCOP", "") if "fits_meta" in meta else meta.get("TELESCOP", "") instrument = meta.get("instrument", "") return telescope == "NDA" or instrument == "NDA" + + @classmethod + def from_raw(cls, header, raw_object): + from astropy.time import Time + + from sunpy.net import attrs as a + + hd_pairs = raw_object + times = Time(hd_pairs[2].data["jd"], format="jd") + freqs = hd_pairs[1].data["frq"].flatten() * u.MHz + data = hd_pairs[2].data["data"] + + res = [] + for i, channel in enumerate(["LL", "RR"]): + meta = { + "fits_meta": hd_pairs[0].header, + "detector": hd_pairs[0].header.get("INSTRUME", "newroutine"), + "instrument": "NDA", + "observatory": "ORN", + "start_time": times[0], + "end_time": times[-1], + "wavelength": a.Wavelength(freqs.min(), freqs.max()), + "times": times, + "freqs": freqs, + "polarisation": channel, + } + res.append(cls(data[:, :, i].T, meta)) + return res diff --git a/radiospectra/spectrogram/sources/psp_rfs.py b/radiospectra/spectrogram/sources/psp_rfs.py index a2d41d3..a8e7bd5 100644 --- a/radiospectra/spectrogram/sources/psp_rfs.py +++ b/radiospectra/spectrogram/sources/psp_rfs.py @@ -1,3 +1,8 @@ +import astropy.units as u +from astropy.time import Time + +from sunpy.net import attrs as a + from radiospectra.spectrogram.meta import SpectrogramMeta from radiospectra.spectrogram.spectrogrambase import GenericSpectrogram @@ -66,7 +71,55 @@ def version(self): return self.meta.version @classmethod - def is_datasource_for(cls, data, meta, **kwargs): + def is_datasource_for(cls, data_or_header, meta_or_raw, **kwargs): + meta = data_or_header if hasattr(data_or_header, "get") else meta_or_raw + if not hasattr(meta, "get"): + return False + + if ( + meta.get("observatory") == "PSP" + and meta.get("instrument") == "FIELDS/RFS" + and meta.get("detector") in ("lfr", "hfr") + ): + return True + + cdf_globals = meta.get("cdf_globals") + if not cdf_globals: + return False return ( - meta["observatory"] == "PSP" and meta["instrument"] == "FIELDS/RFS" and meta["detector"] in ("lfr", "hfr") + cdf_globals.get("Project", "")[0] == "PSP" + and cdf_globals.get("Source_name", [""])[0] == "PSP_FLD>Parker Solar Probe FIELDS" + and "Radio Frequency Spectrometer" in cdf_globals.get("Descriptor", [""])[0] + ) + + @classmethod + def from_raw(cls, header, raw_object): + cdf = raw_object + cdf_globals = header["cdf_globals"] + short, _long = cdf_globals["Descriptor"][0].split(">") + detector = short[4:].lower() + times, data, freqs = ( + cdf.varget(name) + for name in [ + f"epoch_{detector}_auto_averages_ch0_V1V2", + f"psp_fld_l2_rfs_{detector}_auto_averages_ch0_V1V2", + f"frequency_{detector}_auto_averages_ch0_V1V2", + ] + ) + times = Time("J2000.0", scale="tt") + (times << u.ns) + freqs = freqs[0, :] << u.Hz + data = data.T << u.Unit("Volt**2/Hz") + meta = RFSMeta( + { + "cdf_globals": cdf_globals, + "detector": detector, + "instrument": "FIELDS/RFS", + "observatory": "PSP", + "start_time": times[0], + "end_time": times[-1], + "wavelength": a.Wavelength(freqs.min(), freqs.max()), + "times": times, + "freqs": freqs, + } ) + return cls(data, meta) diff --git a/radiospectra/spectrogram/sources/rpw.py b/radiospectra/spectrogram/sources/rpw.py index 4e04721..52985f9 100644 --- a/radiospectra/spectrogram/sources/rpw.py +++ b/radiospectra/spectrogram/sources/rpw.py @@ -1,3 +1,11 @@ +import numpy as np + +import astropy.units as u +from astropy.time import Time + +from sunpy.net import attrs as a +from sunpy.sun.constants import sfu + from radiospectra.spectrogram.meta import SpectrogramMeta from radiospectra.spectrogram.spectrogrambase import GenericSpectrogram @@ -67,5 +75,168 @@ def __init__(self, data, meta, **kwargs): super().__init__(meta=meta, data=data, **kwargs) @classmethod - def is_datasource_for(cls, data, meta, **kwargs): - return meta["instrument"] == "RPW" + def is_datasource_for(cls, data_or_header, meta_or_raw, **kwargs): + meta = data_or_header if hasattr(data_or_header, "get") else meta_or_raw + if not hasattr(meta, "get"): + return False + if meta.get("instrument") == "RPW": + return True + cdf_globals = meta.get("cdf_globals") + if not cdf_globals: + return False + return "SOLO" in cdf_globals.get("Project", "")[0] + + @classmethod + def from_raw(cls, header, raw_object): + cdf = raw_object + cdf_globals = header["cdf_globals"] + data_type = cdf_globals.get("Data_type", [""])[0] + data_descriptor = cdf_globals.get("Descriptor", "")[0] + if "RPW-HFR-SURV" not in data_descriptor and "RPW-TNR-SURV-FLUX" not in data_descriptor: + raise ValueError( + f"Currently radiospectra supports Level 2 HFR survey data " + "and Level 3 HFR, TNR survey data the file is " + f"{cdf_globals.get('Logical_source_description', [''])[0]}" + ) + if "L3" in data_type: + epoch = cdf.varget("Epoch") + times = Time("J2000.0") + epoch * u.ns + freqs = cdf.varget("FREQUENCY") << u.Unit(cdf.varattsget("FREQUENCY")["UNITS"]) + data = cdf.varget("PSD_SFU") + data = np.squeeze(data).T << sfu + detector = cdf_globals.get("Instrument", [""])[0].split(">")[0] + meta = RPWMeta( + { + "cdf_globals": cdf_globals, + "detector": detector, + "instrument": "RPW", + "observatory": "SOLO", + "start_time": times[0], + "end_time": times[-1], + "wavelength": a.Wavelength(freqs.min(), freqs.max()), + "times": times, + "freqs": freqs, + } + ) + return cls(data, meta) + elif "L2" in data_type: + # FREQUENCY_BAND_LABELS = ["HF1", "HF2"] + # SURVEY_MODE_LABELS = ["SURVEY_NORMAL", "SURVEY_BURST"] + # CHANNEL_LABELS = ["1", "2"] + SENSOR_MAPPING = { + 1: "V1", + 2: "V2", + 3: "V3", + 4: "V1-V2", + 5: "V2-V3", + 6: "V3-V1", + 7: "B_MF", + 9: "HF_V1-V2", + 10: "HF_V2-V3", + 11: "HF_V3-V1", + } + + # Extract variables + all_times = Time("J2000.0") + cdf.varget("EPOCH") * u.Unit(cdf.varattsget("EPOCH")["UNITS"]) + all_freqs = cdf.varget("FREQUENCY") << u.Unit(cdf.varattsget("FREQUENCY")["UNITS"]) + + sweep_start_indices = np.asarray(np.diff(cdf.varget("SWEEP_NUM")) != 0).nonzero() + sweep_start_indices = np.insert((sweep_start_indices[0] + 1), 0, 0) + times = all_times[sweep_start_indices] + + sensor = cdf.varget("SENSOR_CONFIG") + np.unique(cdf.varget("FREQUENCY")) + band = cdf.varget("HFR_BAND") + + u.Unit(cdf.varattsget("AGC1").get("UNIT", "V^2/Hz")) + agc1 = cdf.varget("AGC1") + agc2 = cdf.varget("AGC2") + + # Define number of records + n_rec = band.shape[0] + # Get Epoch times of first sample of each sweep in the file + sweep_times = times + nt = len(sweep_times) + # Get complete list of HFR frequency values + hfr_frequency = 375 + 50 * np.arange(321) # This is a guess something between 320 and 324 + nf = len(hfr_frequency) + + # Initialize output 2D array containing voltage spectral power values in V^2/Hz + # Dims = (channels[2], time of the first sweep sample[len(time)], frequency[192]) + specs = np.empty((2, nt, nf)) + # Fill 2D array with NaN for HRF frequencies not actually measured in the file + specs[:] = np.nan + + # Get list of first index of sweeps + isweep = sweep_start_indices[:] + # Get number of sweeps + n_sweeps = len(isweep) + # Insert an element in the end of the isweep list + # containing the end of the latest sweep + # (required for the loop below, in order to have + # a start/end index range for each sweep) + isweep = np.insert(isweep, n_sweeps, n_rec) + + # Initialize sensor_config + sensor_config = np.zeros((2, nt), dtype=object) + tm = [] + # Perform a loop on each sweep + for i in range(n_sweeps): + # Get first and last index of the sweep + i0 = isweep[i] + i1 = isweep[i + 1] + + ts = all_times[i0] + te = all_times[i1 - 1] + tt = (te - ts) * 0.5 + ts + tm.append(tt) + + # Get indices of the actual frequency channels in the frequency vector + freq_indices = ((all_freqs[i0:i1].value - 375) / 50).astype(int) + + # fill output 2D array + specs[0, i, freq_indices] = agc1[i0:i1] + specs[1, i, freq_indices] = agc2[i0:i1] + + # Fill sensor config + sensor_config[0, i] = SENSOR_MAPPING[sensor[i0, 0]] + sensor_config[1, i] = SENSOR_MAPPING[sensor[i0, 1]] + + # Define hfr bands + hfc = np.array(["HF1", "HF2"]) + hfc[band[:100] - 1] + + hfr_frequency = hfr_frequency << u.kHz + + res = [] + if np.any(agc1): + meta1 = RPWMeta( + { + "cdf_globals": cdf_globals, + "detector": "RPW-AGC1", + "instrument": "RPW", + "observatory": "SOLO", + "start_time": times[0], + "end_time": times[-1], + "wavelength": a.Wavelength(hfr_frequency.min(), hfr_frequency.max()), + "times": times, + "freqs": hfr_frequency, + } + ) + res.append(cls(specs[0].T, meta1)) + if np.any(agc2): + meta2 = RPWMeta( + { + "cdf_globals": cdf_globals, + "detector": "RPW-AGC2", + "instrument": "RPW", + "observatory": "SOLO", + "start_time": times[0], + "end_time": times[-1], + "wavelength": a.Wavelength(hfr_frequency.min(), hfr_frequency.max()), + "times": times, + "freqs": hfr_frequency, + } + ) + res.append(cls(specs[1].T, meta2)) + return res diff --git a/radiospectra/spectrogram/sources/rstn.py b/radiospectra/spectrogram/sources/rstn.py index 91729b3..ddec6c7 100644 --- a/radiospectra/spectrogram/sources/rstn.py +++ b/radiospectra/spectrogram/sources/rstn.py @@ -1,3 +1,14 @@ +import gzip +import struct + +import numpy as np +import pandas as pd + +import astropy.units as u +from astropy.time import Time + +from sunpy.net import attrs as a + from radiospectra.spectrogram.meta import SpectrogramMeta from radiospectra.spectrogram.spectrogrambase import GenericSpectrogram @@ -34,5 +45,93 @@ def __init__(self, data, meta, **kwargs): super().__init__(meta=meta, data=data, **kwargs) @classmethod - def is_datasource_for(cls, data, meta, **kwargs): - return meta["instrument"] == "RSTN" + def is_datasource_for(cls, data_or_header, meta_or_raw, **kwargs): + meta = data_or_header if hasattr(data_or_header, "get") else meta_or_raw + if not hasattr(meta, "get"): + return False + return meta.get("instrument") == "RSTN" + + @classmethod + def from_raw(cls, header, raw_object): + file = header["file_path"] + with file.open("rb") as buff: + data = buff.read() + if file.suffixes[-1] == ".gz": + data = gzip.decompress(data) + # Data is store as a series of records made of different numbers of bytes + # General header information + # 1 Year (last 2 digits) Byte integer (unsigned) + # 2 Month number (1 to 12) " + # 3 Day (1 to 31) " + # 4 Hour (0 to 23 UT) " + # 5 Minute (0 to 59) " + # 6 Second at start of scan (0 to 59) " + # 7 Site Number (0 to 255) " + # 8 Number of bands in the record (2) " + # + # Band 1 (A-band) header information + # 9,10 Start Frequency (MHz) Word integer (16 bits) + # 11,12 End Frequency (MHz) " + # 13,14 Number of bytes in data record (401)" + # 15 Analyser reference level Byte integer + # 16 Analyser attenuation (dB) " + # + # Band 2 (B-band) header information + # 17-24 As for band 1 + # + # Spectrum Analyser data + # 25-425 401 data bytes for band 1 (A-band) + # 426-826 401 data bytes for band 2 (B-band) + record_struc = struct.Struct("B" * 8 + "H" * 3 + "B" * 2 + "H" * 3 + "B" * 2 + "B" * 401 + "B" * 401) + records = record_struc.iter_unpack(data) + # Map of numeric records to locations + site_map = {1: "Palehua", 2: "Holloman", 3: "Learmonth", 4: "San Vito"} + df = pd.DataFrame([(*r[:18], np.array(r[18:419]), np.array(r[419:820])) for r in records]) + df.columns = [ + "year", + "month", + "day", + "hour", + "minute", + "second", + "site", + " num_bands", + "start_freq1", + "end_freq1", + "num_bytes1", + "analyser_ref1", + "analyser_atten1", + "start_freq2", + "end_freq2", + "num_bytes2", + "analyser_ref2", + "analyser_atten2", + "spec1", + "spec2", + ] + # Hack to make to_datetime work - earliest dates seem to be 2000 and won't be + # around in 3000! + df["year"] = df["year"] + 2000 + df["time"] = pd.to_datetime(df[["year", "month", "day", "hour", "minute", "second"]]) + # Equations taken from document + n = np.arange(1, 402) + freq_a = (25 + 50 * (n - 1) / 400) * u.MHz + freq_b = (75 + 105 * (n - 1) / 400) * u.MHz + freqs = np.hstack([freq_a, freq_b]) + data = np.hstack([np.vstack(df[name].to_numpy()) for name in ["spec1", "spec2"]]).T + times = Time( + Time(df["time"]), format="iso" + ) # TODO update once datetime format is supported by current plotters + meta = RSTNMeta( + { + "instrument": "RSTN", + "observatory": site_map[df["site"][0]], + "start_time": times[0], + "end_time": times[-1], + "detector": "RSTN", + "wavelength": a.Wavelength(freqs[0], freqs[-1]), + "freqs": freqs, + "times": times, + } + ) + return cls(data, meta) diff --git a/radiospectra/spectrogram/sources/swaves.py b/radiospectra/spectrogram/sources/swaves.py index faed71f..0bf1a9f 100644 --- a/radiospectra/spectrogram/sources/swaves.py +++ b/radiospectra/spectrogram/sources/swaves.py @@ -1,3 +1,10 @@ +import numpy as np + +import astropy.units as u +from astropy.time import Time + +from sunpy.net import attrs as a + from radiospectra.spectrogram.meta import SpectrogramMeta from radiospectra.spectrogram.spectrogrambase import GenericSpectrogram @@ -36,5 +43,39 @@ def __init__(self, data, meta, **kwargs): super().__init__(meta=meta, data=data, **kwargs) @classmethod - def is_datasource_for(cls, data, meta, **kwargs): - return meta["instrument"] == "swaves" + def is_datasource_for(cls, data_or_header, meta_or_raw, **kwargs): + meta = data_or_header if hasattr(data_or_header, "get") else meta_or_raw + if not hasattr(meta, "get"): + return False + return meta.get("instrument", "").lower() == "swaves" + + @classmethod + def from_raw(cls, header, raw_object): + file = header["file_path"] + name, prod, date, spacecraft, receiver = file.stem.split("_") + # frequency range + freqs = np.genfromtxt(file, max_rows=1) * u.kHz + # bg which is already subtracted from data + bg = np.genfromtxt(file, skip_header=1, max_rows=1) + # data + data = np.genfromtxt(file, skip_header=2) + times = data[:, 0] * u.min + data = data[:, 1:].T + start_time = Time.strptime(date, "%Y%m%d") + end_time = start_time + times[-1] + times = start_time + times + meta = SWAVESMeta( + { + "instrument": name, + "observatory": f"STEREO {spacecraft.upper()}", + "product": prod, + "start_time": start_time, + "end_time": end_time, + "wavelength": a.Wavelength(freqs[0], freqs[-1]), + "detector": receiver, + "freqs": freqs, + "times": times, + "background": bg, + } + ) + return cls(data, meta) diff --git a/radiospectra/spectrogram/sources/tests/test_callisto.py b/radiospectra/spectrogram/sources/tests/test_callisto.py index 2843bbd..29844c4 100644 --- a/radiospectra/spectrogram/sources/tests/test_callisto.py +++ b/radiospectra/spectrogram/sources/tests/test_callisto.py @@ -23,6 +23,12 @@ def test_callisto(parse_path_moc): start_time = Time("2011-06-07 06:24:00.213") meta = { "fits_meta": { + "CONTENT": "e-CALLISTO", + "DATE-OBS": "2011/06/07", + "TIME-OBS": "06:24:00.213", + "DATE-END": "2011/06/07", + "TIME-END": "06:39:00", + "INSTRUME": "BIR", "OBS_LAC": "N", "OBS_LAT": 53.0941390991211, "OBS_LOC": "E", @@ -241,7 +247,15 @@ def test_callisto(parse_path_moc): * u.MHz, } array = np.zeros((200, 3600)) - parse_path_moc.return_value = [(array, meta)] + hdu0 = MagicMock() + hdu1 = MagicMock() + hdu0.header = meta["fits_meta"] + hdu0.data = array + hdu1.data = { + "TIME": np.arange(3600) * 0.25, + "FREQUENCY": np.array(meta["freqs"].value), + } + parse_path_moc.return_value = [(meta["fits_meta"], [hdu0, hdu1])] file = Path("fake.fit.gz") spec = Spectrogram(file) assert isinstance(spec, CALISTOSpectrogram) diff --git a/radiospectra/spectrogram/sources/tests/test_waves.py b/radiospectra/spectrogram/sources/tests/test_waves.py index 045f1b6..4cba969 100644 --- a/radiospectra/spectrogram/sources/tests/test_waves.py +++ b/radiospectra/spectrogram/sources/tests/test_waves.py @@ -7,9 +7,6 @@ import pytest import astropy.units as u -from astropy.time import Time - -from sunpy.net import attrs as a from radiospectra.spectrogram import Spectrogram from radiospectra.spectrogram.sources import WAVESSpectrogram @@ -20,18 +17,10 @@ @mock.patch("radiospectra.spectrogram.spectrogram_factory.parse_path") def test_waves_rad1(parse_path_moc): - meta = { - "instrument": "WAVES", - "observatory": "wind", - "start_time": Time("2020-11-28 00:00:00"), - "end_time": Time("2020-11-28 23:59:00"), - "wavelength": a.Wavelength(20 * u.kHz, 1040 * u.kHz), - "detector": "rad1", - "freqs": np.linspace(20, 1040, 256) * u.kHz, - "times": np.arange(1440) * u.min, - } - array = np.zeros((256, 1440)) - parse_path_moc.return_value = [(array, meta)] + header = {"instrument": "waves", "file_type": "idl_sav", "file_path": Path("wind_waves_rad1_20201128.R1")} + array = np.zeros((256, 1441)) + raw_object = {"arrayb": array} + parse_path_moc.return_value = [(header, raw_object)] file = Path("fake.r1") spec = Spectrogram(file) assert isinstance(spec, WAVESSpectrogram) @@ -39,25 +28,17 @@ def test_waves_rad1(parse_path_moc): assert spec.instrument == "WAVES" assert spec.detector == "RAD1" assert spec.start_time.datetime == datetime(2020, 11, 28, 0, 0) - assert spec.end_time.datetime == datetime(2020, 11, 28, 23, 59) + assert spec.end_time.datetime == datetime(2020, 11, 28, 23, 59, 59) assert spec.wavelength.min == 20.0 * u.kHz assert spec.wavelength.max == 1040.0 * u.kHz @mock.patch("radiospectra.spectrogram.spectrogram_factory.parse_path") def test_waves_rad2(parse_path_moc): - meta = { - "instrument": "WAVES", - "observatory": "WIND", - "start_time": Time("2020-11-28 00:00:00"), - "end_time": Time("2020-11-28 23:59:00"), - "wavelength": a.Wavelength(1.075 * u.MHz, 13.825 * u.MHz), - "detector": "RAD2", - "freqs": np.linspace(1.075, 13.825, 256) * u.MHz, - "times": np.arange(1440) * u.min, - } - array = np.zeros((319, 1440)) - parse_path_moc.return_value = [(array, meta)] + header = {"instrument": "waves", "file_type": "idl_sav", "file_path": Path("wind_waves_rad2_20201128.R2")} + array = np.zeros((256, 1441)) + raw_object = {"arrayb": array} + parse_path_moc.return_value = [(header, raw_object)] file = Path("fake.dat") spec = Spectrogram(file) assert isinstance(spec, WAVESSpectrogram) @@ -65,7 +46,7 @@ def test_waves_rad2(parse_path_moc): assert spec.instrument == "WAVES" assert spec.detector == "RAD2" assert spec.start_time.datetime == datetime(2020, 11, 28, 0, 0) - assert spec.end_time.datetime == datetime(2020, 11, 28, 23, 59) + assert spec.end_time.datetime == datetime(2020, 11, 28, 23, 59, 59) assert spec.wavelength.min == 1.075 * u.MHz assert spec.wavelength.max == 13.825 * u.MHz @@ -75,9 +56,10 @@ def test_waves_prefixed_filename_parses_date(readsav_mock): data_array = np.zeros((256, 1441)) readsav_mock.return_value = {"arrayb": data_array} - _, meta = SpectrogramFactory._read_idl_sav(Path("wind_waves_rad1_20200711.R1"), instrument="waves") + header, raw_object = SpectrogramFactory._read_idl_sav(Path("wind_waves_rad1_20200711.R1"), instrument="waves") + spec = WAVESSpectrogram.from_raw(header, raw_object) - assert meta["start_time"].isot == "2020-07-11T00:00:00.000" + assert spec.start_time.isot == "2020-07-11T00:00:00.000" @pytest.mark.remote_data diff --git a/radiospectra/spectrogram/sources/waves.py b/radiospectra/spectrogram/sources/waves.py index a09f526..50b6c34 100644 --- a/radiospectra/spectrogram/sources/waves.py +++ b/radiospectra/spectrogram/sources/waves.py @@ -1,3 +1,10 @@ +import numpy as np + +import astropy.units as u +from astropy.time import Time + +from sunpy.net import attrs as a + from radiospectra.spectrogram.meta import SpectrogramMeta from radiospectra.spectrogram.spectrogrambase import GenericSpectrogram @@ -36,5 +43,44 @@ def __init__(self, data, meta, **kwargs): super().__init__(meta=meta, data=data, **kwargs) @classmethod - def is_datasource_for(cls, data, meta, **kwargs): - return meta.get("instrument", None) == "WAVES" + def is_datasource_for(cls, data_or_header, meta_or_raw, **kwargs): + meta = data_or_header if hasattr(data_or_header, "get") else meta_or_raw + if not hasattr(meta, "get"): + return False + if meta.get("instrument") == "WAVES": + return True + return meta.get("file_type") == "idl_sav" and meta.get("instrument") == "waves" + + @classmethod + def from_raw(cls, header, raw_object): + file = header.get("file_path") + data = raw_object + data_array = data["arrayb"] + if file.suffix == ".R1": + freqs = np.linspace(20, 1040, 256) * u.kHz + receiver = "RAD1" + elif file.suffix == ".R2": + freqs = np.linspace(1.075, 13.825, 256) * u.MHz + receiver = "RAD2" + else: + raise ValueError(f"Unknown WIND/WAVES file type: {file.suffix}") + + bg = data_array[:, -1] + data_vals = data_array[:, :-1] + start_time = Time.strptime(file.stem.split("_")[-1], "%Y%m%d") + end_time = start_time + 86399 * u.s + times = start_time + (np.arange(1440) * 60 + 30) * u.s + meta = WAVESMeta( + { + "instrument": "WAVES", + "observatory": "WIND", + "start_time": start_time, + "end_time": end_time, + "wavelength": a.Wavelength(freqs[0], freqs[-1]), + "detector": receiver, + "freqs": freqs, + "times": times, + "background": bg, + } + ) + return cls(data_vals, meta) diff --git a/radiospectra/spectrogram/spectrogram_factory.py b/radiospectra/spectrogram/spectrogram_factory.py index 4e88f0a..6a795f6 100644 --- a/radiospectra/spectrogram/spectrogram_factory.py +++ b/radiospectra/spectrogram/spectrogram_factory.py @@ -1,5 +1,3 @@ -import gzip -import struct import pathlib import warnings import functools @@ -9,18 +7,14 @@ import cdflib import numpy as np -import pandas as pd from scipy.io import readsav import astropy.units as u from astropy.io import fits from astropy.io.fits import Header -from astropy.time import Time -from sunpy import log from sunpy.data import cache from sunpy.net import attrs as a -from sunpy.sun.constants import sfu from sunpy.time import parse_time from sunpy.util.datatype_factory_base import ( BasicRegistrationFactory, @@ -30,12 +24,10 @@ ) from sunpy.util.exceptions import SunpyUserWarning, warn_user from sunpy.util.io import is_url, parse_path, possibly_a_path -from sunpy.util.metadata import MetaDict from sunpy.util.util import expand_list from radiospectra.exceptions import NoSpectrogramInFileError, SpectraMetaValidationError from radiospectra.spectrogram.spectrogrambase import GenericSpectrogram -from radiospectra.utils import subband_to_freq SUPPORTED_ARRAY_TYPES = (np.ndarray,) try: @@ -189,16 +181,16 @@ def __call__(self, *args, silence_errors=False, **kwargs): """ data_header_pairs = self._parse_args(*args, silence_errors=silence_errors, **kwargs) new_maps = list() - # Loop over each registered type and check to see if WidgetType - # matches the arguments. If it does, use that type. for pair in data_header_pairs: if isinstance(pair, GenericSpectrogram): new_maps.append(pair) continue - data, header = pair - meta = MetaDict(header) + # Detect whether the pair is (header, raw_object) or (data, meta). + # If the first element is a dict or FITS Header, it's a raw pair. + first = pair[0] + is_raw = isinstance(first, (dict, Header)) try: - new_map = self._check_registered_widgets(data, meta, **kwargs) + new_map = self._check_registered_widgets(pair[0], pair[1], from_raw=is_raw, **kwargs) new_maps.append(new_map) except (NoMatchError, MultipleMatchError, ValidationFunctionError, SpectraMetaValidationError) as e: if not silence_errors: @@ -210,12 +202,15 @@ def __call__(self, *args, silence_errors=False, **kwargs): return new_maps[0] return new_maps - def _check_registered_widgets(self, data, meta, **kwargs): + def _check_registered_widgets(self, data_or_header, meta_or_raw, from_raw=False, **kwargs): candidate_widget_types = list() for key in self.registry: - # Call the registered validation function for each registered class - if self.registry[key](data, meta, **kwargs): - candidate_widget_types.append(key) + try: + if self.registry[key](data_or_header, meta_or_raw, **kwargs): + candidate_widget_types.append(key) + except (KeyError, TypeError, IndexError, AttributeError): + # Validation function crashed — this widget doesn't match + continue n_matches = len(candidate_widget_types) if n_matches == 0: @@ -231,9 +226,10 @@ def _check_registered_widgets(self, data, meta, **kwargs): "identification." ) - # Only one is found WidgetType = candidate_widget_types[0] - return WidgetType(data, meta, **kwargs) + if from_raw: + return WidgetType.from_raw(data_or_header, meta_or_raw) + return WidgetType(data_or_header, meta_or_raw, **kwargs) def _read_file(self, file, **kwargs): file = Path(file) @@ -260,406 +256,35 @@ def _read_file(self, file, **kwargs): @staticmethod def _read_dat(file): if "swaves" in file.name: - name, prod, date, spacecraft, receiver = file.stem.split("_") - # frequency range - freqs = np.genfromtxt(file, max_rows=1) * u.kHz - # bg which is already subtracted from data - bg = np.genfromtxt(file, skip_header=1, max_rows=1) - # data - data = np.genfromtxt(file, skip_header=2) - times = data[:, 0] * u.min - data = data[:, 1:].T - meta = { - "instrument": name, - "observatory": f"STEREO {spacecraft.upper()}", - "product": prod, - "start_time": Time.strptime(date, "%Y%m%d"), - "wavelength": a.Wavelength(freqs[0], freqs[-1]), - "detector": receiver, - "freqs": freqs, - "background": bg, - } - meta["times"] = meta["start_time"] + times - meta["end_time"] = meta["start_time"] + times[-1] - return data, meta + header = {"file_type": "dat", "instrument": "swaves", "file_path": file} + return [(header, None)] elif "bst" in file.name: - subbands = (np.arange(54, 454, 2), np.arange(54, 454, 2), np.arange(54, 230, 2)) - num_subbands = 488 - - data = np.fromfile(file) - polarisation = file.stem[-1] - - num_times = data.shape[0] / num_subbands - if not num_times.is_integer(): - log.warning("BST file seems incomplete dropping incomplete frequencies") - num_times = np.floor(num_times).astype(int) - truncate = num_times * num_subbands - data = data[:truncate] - data = data.reshape(-1, num_subbands).T # (Freq x Time).T = (Time x Freq) - dt = np.arange(num_times) * 1 * u.s - start_time = Time.strptime(file.name.split("_bst")[0], "%Y%m%d_%H%M%S") - times = start_time + dt - - obs_mode = (3, 5, 7) - - freqs = [subband_to_freq(sb, mode) for sb, mode in zip(subbands, obs_mode)] - - # 1st 200 sbs mode 3, next 200 sbs mode 5, last 88 sbs mode 7 - spec = {0: data[:200, :], 1: data[200:400, :], 2: data[400:, :]} - data_header_pairs = [] - for i in range(3): - meta = { - "instrument": "ILOFAR", - "observatory": "Birr (IE613)", - "start_time": times[0], - "mode": obs_mode[i], - "wavelength": a.Wavelength(freqs[i][0], freqs[i][-1]), - "freqs": freqs[i], - "times": times, - "end_time": times[-1], - "detector": "ILOFAR", - "polarisation": polarisation, - } - - data_header_pairs.append((spec[i], meta)) - return data_header_pairs + header = {"file_type": "dat", "instrument": "ILOFAR", "file_path": file} + return [(header, None)] else: raise ValueError(f"File {file} not supported.") @staticmethod def _read_srs(file): - with file.open("rb") as buff: - data = buff.read() - if file.suffixes[-1] == ".gz": - data = gzip.decompress(data) - # Data is store as a series of records made of different numbers of bytes - # General header information - # 1 Year (last 2 digits) Byte integer (unsigned) - # 2 Month number (1 to 12) " - # 3 Day (1 to 31) " - # 4 Hour (0 to 23 UT) " - # 5 Minute (0 to 59) " - # 6 Second at start of scan (0 to 59) " - # 7 Site Number (0 to 255) " - # 8 Number of bands in the record (2) " - # - # Band 1 (A-band) header information - # 9,10 Start Frequency (MHz) Word integer (16 bits) - # 11,12 End Frequency (MHz) " - # 13,14 Number of bytes in data record (401)" - # 15 Analyser reference level Byte integer - # 16 Analyser attenuation (dB) " - # - # Band 2 (B-band) header information - # 17-24 As for band 1 - # - # Spectrum Analyser data - # 25-425 401 data bytes for band 1 (A-band) - # 426-826 401 data bytes for band 2 (B-band) - record_struc = struct.Struct("B" * 8 + "H" * 3 + "B" * 2 + "H" * 3 + "B" * 2 + "B" * 401 + "B" * 401) - records = record_struc.iter_unpack(data) - # Map of numeric records to locations - site_map = {1: "Palehua", 2: "Holloman", 3: "Learmonth", 4: "San Vito"} - df = pd.DataFrame([(*r[:18], np.array(r[18:419]), np.array(r[419:820])) for r in records]) - df.columns = [ - "year", - "month", - "day", - "hour", - "minute", - "second", - "site", - " num_bands", - "start_freq1", - "end_freq1", - "num_bytes1", - "analyser_ref1", - "analyser_atten1", - "start_freq2", - "end_freq2", - "num_bytes2", - "analyser_ref2", - "analyser_atten2", - "spec1", - "spec2", - ] - # Hack to make to_datetime work - earliest dates seem to be 2000 and won't be - # around in 3000! - df["year"] = df["year"] + 2000 - df["time"] = pd.to_datetime(df[["year", "month", "day", "hour", "minute", "second"]]) - # Equations taken from document - n = np.arange(1, 402) - freq_a = (25 + 50 * (n - 1) / 400) * u.MHz - freq_b = (75 + 105 * (n - 1) / 400) * u.MHz - freqs = np.hstack([freq_a, freq_b]) - data = np.hstack([np.vstack(df[name].to_numpy()) for name in ["spec1", "spec2"]]).T - times = Time( - Time(df["time"]), format="iso" - ) # TODO update once datetime format is supported by current plotters - meta = { - "instrument": "RSTN", - "observatory": site_map[df["site"][0]], - "start_time": times[0], - "end_time": times[-1], - "detector": "RSTN", - "wavelength": a.Wavelength(freqs[0], freqs[-1]), - "freqs": freqs, - "times": times, - } - return data, meta + header = {"file_type": "srs", "instrument": "RSTN", "file_path": file} + return header, None @staticmethod def _read_cdf(file): cdf = cdflib.CDF(file) - cdf_globals = cdf.globalattsget() - - if ( - cdf_globals.get("Project", "")[0] == "PSP" - and cdf_globals.get("Source_name")[0] == "PSP_FLD>Parker Solar Probe FIELDS" - and "Radio Frequency Spectrometer" in cdf_globals.get("Descriptor")[0] - ): - short, _long = cdf_globals["Descriptor"][0].split(">") - - detector = short[4:].lower() - times, data, freqs = ( - cdf.varget(name) - for name in [ - f"epoch_{detector}_auto_averages_ch0_V1V2", - f"psp_fld_l2_rfs_{detector}_auto_averages_ch0_V1V2", - f"frequency_{detector}_auto_averages_ch0_V1V2", - ] - ) - times = Time("J2000.0", scale="tt") + (times << u.ns) - freqs = freqs[0, :] << u.Hz - data = data.T << u.Unit("Volt**2/Hz") - meta = { - "cdf_globals": cdf_globals, - "detector": detector, - "instrument": "FIELDS/RFS", - "observatory": "PSP", - "start_time": times[0], - "end_time": times[-1], - "wavelength": a.Wavelength(freqs.min(), freqs.max()), - "times": times, - "freqs": freqs, - } - return data, meta - elif "SOLO" in cdf_globals.get("Project", "")[0]: - data_type = cdf_globals.get("Data_type", [""])[0] - data_descriptor = cdf_globals.get("Descriptor", "")[0] - if "RPW-HFR-SURV" not in data_descriptor and "RPW-TNR-SURV-FLUX" not in data_descriptor: - raise ValueError( - f"Currently radiospectra supports Level 2 HFR survey data " - "and Level 3 HFR, TNR survey data the file " - f"{file.name} is {cdf_globals.get('Logical_source_description', [''])[0]}" - ) - if "L3" in data_type: - epoch = cdf.varget("Epoch") - times = Time("J2000.0") + epoch * u.ns - freqs = cdf.varget("FREQUENCY") << u.Unit(cdf.varattsget("FREQUENCY")["UNITS"]) - data = cdf.varget("PSD_SFU") - data = np.squeeze(data).T << sfu - detector = cdf_globals.get("Instrument", [""])[0].split(">")[0] - meta = { - "cdf_globals": cdf_globals, - "detector": detector, - "instrument": "RPW", - "observatory": "SOLO", - "start_time": times[0], - "end_time": times[-1], - "wavelength": a.Wavelength(freqs.min(), freqs.max()), - "times": times, - "freqs": freqs, - } - return data, meta - elif "L2" in data_type: - # FREQUENCY_BAND_LABELS = ["HF1", "HF2"] - # SURVEY_MODE_LABELS = ["SURVEY_NORMAL", "SURVEY_BURST"] - # CHANNEL_LABELS = ["1", "2"] - SENSOR_MAPPING = { - 1: "V1", - 2: "V2", - 3: "V3", - 4: "V1-V2", - 5: "V2-V3", - 6: "V3-V1", - 7: "B_MF", - 9: "HF_V1-V2", - 10: "HF_V2-V3", - 11: "HF_V3-V1", - } - - # Extract variables - all_times = Time("J2000.0") + cdf.varget("EPOCH") * u.Unit(cdf.varattsget("EPOCH")["UNITS"]) - all_freqs = cdf.varget("FREQUENCY") << u.Unit(cdf.varattsget("FREQUENCY")["UNITS"]) - - sweep_start_indices = np.asarray(np.diff(cdf.varget("SWEEP_NUM")) != 0).nonzero() - sweep_start_indices = np.insert((sweep_start_indices[0] + 1), 0, 0) - times = all_times[sweep_start_indices] - - sensor = cdf.varget("SENSOR_CONFIG") - np.unique(cdf.varget("FREQUENCY")) - band = cdf.varget("HFR_BAND") - - u.Unit(cdf.varattsget("AGC1").get("UNIT", "V^2/Hz")) - agc1 = cdf.varget("AGC1") - agc2 = cdf.varget("AGC2") - - # Define number of records - n_rec = band.shape[0] - # Get Epoch times of first sample of each sweep in the file - sweep_times = times - nt = len(sweep_times) - # Get complete list of HFR frequency values - hfr_frequency = 375 + 50 * np.arange(321) # This is a guess something between 320 and 324 - nf = len(hfr_frequency) - - # Initialize output 2D array containing voltage spectral power values in V^2/Hz - # Dims = (channels[2], time of the first sweep sample[len(time)], frequency[192]) - specs = np.empty((2, nt, nf)) - # Fill 2D array with NaN for HRF frequencies not actually measured in the file - specs[:] = np.nan - - # Get list of first index of sweeps - isweep = sweep_start_indices[:] - # Get number of sweeps - n_sweeps = len(isweep) - # Insert an element in the end of the isweep list - # containing the end of the latest sweep - # (required for the loop below, in order to have - # a start/end index range for each sweep) - isweep = np.insert(isweep, n_sweeps, n_rec) - - # Initialize sensor_config - sensor_config = np.zeros((2, nt), dtype=object) - tm = [] - # Perform a loop on each sweep - for i in range(n_sweeps): - # Get first and last index of the sweep - i0 = isweep[i] - i1 = isweep[i + 1] - - ts = all_times[i0] - te = all_times[i1 - 1] - tt = (te - ts) * 0.5 + ts - tm.append(tt) - - # Get indices of the actual frequency channels in the frequency vector - freq_indices = ((all_freqs[i0:i1].value - 375) / 50).astype(int) - - # fill output 2D array - specs[0, i, freq_indices] = agc1[i0:i1] - specs[1, i, freq_indices] = agc2[i0:i1] - - # Fill sensor config - sensor_config[0, i] = SENSOR_MAPPING[sensor[i0, 0]] - sensor_config[1, i] = SENSOR_MAPPING[sensor[i0, 1]] - - # Define hfr bands - hfc = np.array(["HF1", "HF2"]) - hfc[band[:100] - 1] - - hfr_frequency = hfr_frequency << u.kHz - - res = [] - if np.any(agc1): - meta1 = { - "cdf_globals": cdf_globals, - "detector": "RPW-AGC1", - "instrument": "RPW", - "observatory": "SOLO", - "start_time": times[0], - "end_time": times[-1], - "wavelength": a.Wavelength(hfr_frequency.min(), hfr_frequency.max()), - "times": times, - "freqs": hfr_frequency, - } - res.append((specs[0].T, meta1)) - if np.any(agc2): - meta2 = { - "cdf_globals": cdf_globals, - "detector": "RPW-AGC2", - "instrument": "RPW", - "observatory": "SOLO", - "start_time": times[0], - "end_time": times[-1], - "wavelength": a.Wavelength(hfr_frequency.min(), hfr_frequency.max()), - "times": times, - "freqs": hfr_frequency, - } - res.append((specs[1].T, meta2)) - return res + header = {"file_type": "cdf", "file_path": file, "cdf_globals": cdf_globals} + return header, cdf @staticmethod def _read_fits(file): hd_pairs = fits.open(file) if "e-CALLISTO" in hd_pairs[0].header.get("CONTENT", ""): - data = hd_pairs[0].data - times = hd_pairs[1].data["TIME"].flatten() * u.s - freqs = hd_pairs[1].data["FREQUENCY"].flatten() * u.MHz - start_time = parse_time(hd_pairs[0].header["DATE-OBS"] + " " + hd_pairs[0].header["TIME-OBS"]) - try: - end_time = parse_time(hd_pairs[0].header["DATE-END"] + " " + hd_pairs[0].header["TIME-END"]) - except ValueError: - # See https://github.com/sunpy/radiospectra/issues/74 - time_comps = hd_pairs[0].header["TIME-END"].split(":") - time_comps[0] = "00" - fixed_time = ":".join(time_comps) - date_offset = parse_time(hd_pairs[0].header["DATE-END"] + " " + fixed_time) - end_time = date_offset + 1 * u.day - - times = start_time + times - meta = { - "fits_meta": hd_pairs[0].header, - "detector": "e-CALLISTO", - "instrument": "e-CALLISTO", - "observatory": hd_pairs[0].header["INSTRUME"], - "start_time": start_time, - "end_time": end_time, - "wavelength": a.Wavelength(freqs.min(), freqs.max()), - "times": times, - "freqs": freqs, - } - return data, meta + return hd_pairs[0].header, hd_pairs elif hd_pairs[0].header.get("TELESCOP", "") == "EOVSA": - times = Time(hd_pairs[2].data["mjd"] + hd_pairs[2].data["time"] / 1000.0 / 86400.0, format="mjd") - freqs = hd_pairs[1].data["sfreq"] * u.GHz - data = hd_pairs[0].data - start_time = parse_time(hd_pairs[0].header["DATE_OBS"]) - end_time = parse_time(hd_pairs[0].header["DATE_END"]) - meta = { - "fits_meta": hd_pairs[0].header, - "detector": "EOVSA", - "instrument": "EOVSA", - "observatory": "Owens Valley", - "start_time": start_time, - "end_time": end_time, - "wavelength": a.Wavelength(freqs.min(), freqs.max()), - "times": times, - "freqs": freqs, - } - return data, meta + return hd_pairs[0].header, hd_pairs elif hd_pairs[0].header.get("TELESCOP", "") == "NDA": - times = Time(hd_pairs[2].data["jd"], format="jd") - freqs = hd_pairs[1].data["frq"].flatten() * u.MHz - data = hd_pairs[2].data["data"] - res = [] - for i, channel in enumerate(["LL", "RR"]): - meta = { - "fits_meta": hd_pairs[0].header, - "detector": hd_pairs[0].header.get("INSTRUME", "newroutine"), - "instrument": "NDA", - "observatory": "ORN", - "start_time": times[0], - "end_time": times[-1], - "wavelength": a.Wavelength(freqs.min(), freqs.max()), - "times": times, - "freqs": freqs, - "polarisation": channel, - } - res.append((data[:, :, i].T, meta)) - return res + return hd_pairs[0].header, hd_pairs # Semi standard - spec in primary and time and freq in 1st extension try: data = hd_pairs[0].data @@ -689,38 +314,9 @@ def _read_fits(file): @staticmethod def _read_idl_sav(file, instrument=None): data = readsav(file) - if instrument == "waves": - # See https://solar-radio.gsfc.nasa.gov/wind/one_minute_doc.html - data_array = data["arrayb"] - # frequency range - if file.suffix == ".R1": - freqs = np.linspace(20, 1040, 256) * u.kHz - receiver = "RAD1" - elif file.suffix == ".R2": - freqs = np.linspace(1.075, 13.825, 256) * u.MHz - receiver = "RAD2" - else: - raise ValueError(f"Unknown WIND/WAVES file type: {file.suffix}") - # bg which is already subtracted from data ? - bg = data_array[:, -1] - data = data_array[:, :-1] - start_time = Time.strptime(file.stem.split("_")[-1], "%Y%m%d") - end_time = start_time + 86399 * u.s - times = start_time + (np.arange(1440) * 60 + 30) * u.s - meta = { - "instrument": "WAVES", - "observatory": "WIND", - "start_time": start_time, - "end_time": end_time, - "wavelength": a.Wavelength(freqs[0], freqs[-1]), - "detector": receiver, - "freqs": freqs, - "times": times, - "background": bg, - } - return data, meta - else: - raise ValueError(f"Unrecognized IDL .save file: {file}") + # Return (header, raw_object) — parsing moved to WAVESSpectrogram.from_raw + header = {"file_type": "idl_sav", "instrument": instrument, "file_path": file} + return header, data Spectrogram = SpectrogramFactory(registry=GenericSpectrogram._registry, default_widget_type=GenericSpectrogram) diff --git a/radiospectra/spectrogram/spectrogrambase.py b/radiospectra/spectrogram/spectrogrambase.py index 200e01f..885c7bc 100644 --- a/radiospectra/spectrogram/spectrogrambase.py +++ b/radiospectra/spectrogram/spectrogrambase.py @@ -122,6 +122,13 @@ def _validate_meta(self, meta): if err_message: raise SpectraMetaValidationError("\n".join(err_message)) + @classmethod + def from_raw(cls, header, raw_object): + """ + Parse raw file objects into a spectrogram. Override in subclasses. + """ + raise NotImplementedError(f"{cls.__name__} does not implement from_raw") + @staticmethod def _time_axis_from_meta(meta): times = meta["times"]