diff --git a/changelog/148.feature.rst b/changelog/148.feature.rst new file mode 100644 index 0000000..80e68bd --- /dev/null +++ b/changelog/148.feature.rst @@ -0,0 +1 @@ +Added support for Solar Orbiter RPW Level 3 TNR and HFR survey data products. diff --git a/radiospectra/spectrogram/sources/tests/test_solo_rpw.py b/radiospectra/spectrogram/sources/tests/test_solo_rpw.py new file mode 100644 index 0000000..5858def --- /dev/null +++ b/radiospectra/spectrogram/sources/tests/test_solo_rpw.py @@ -0,0 +1,118 @@ +from pathlib import Path +from datetime import datetime +from unittest import mock + +import numpy as np + +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 RPWSpectrogram + + +@mock.patch("radiospectra.spectrogram.spectrogram_factory.parse_path") +def test_solo_rpw_tnr(parse_path_moc): + start_time = Time("2024-03-23 00:00:00.000") + end_time = Time("2024-03-24 00:00:00.000") + + meta = { + "cdf_meta": { + "Project": "SOLO>Solar Orbiter", + "Source_name": "SOLO>Solar Orbiter", + "Descriptor": "RPW-TNR-SURV-FLUX>Radio and Plasma Waves, Thermal Noise Receiver, Survey mode, spectral data in physical units", + "Data_type": "L3>Level 3 Data", + "Data_version": "01", + "Logical_file_id": "solo_L3_rpw-tnr-surv-flux_20240323_V01", + "Instrument": "RPW", + }, + "detector": "RPW-TNR", + "instrument": "RPW", + "observatory": "SOLO", + "start_time": start_time, + "end_time": end_time, + "wavelength": a.Wavelength(3.992000102996826 * u.kHz, 978.572021484375 * u.kHz), + "times": start_time + np.linspace(0, (end_time - start_time).to_value("s"), 7322) * u.s, + "freqs": np.array([ 3992, 4169, 4353, 4546, 4747, 4957, 5177, 5406, 5645, + 5895, 6156, 6429, 6713, 7011, 7321, 7645, 7984, 8337, + 8706, 9092, 9494, 9914, 10353, 10812, 11290, 11790, 12312, + 12857, 13427, 14021, 14642, 15290, 15967, 16674, 17412, 18183, + 18988, 19829, 20707, 21624, 22581, 23581, 24625, 25715, 26853, + 28042, 29284, 30580, 31934, 33348, 34825, 36366, 37976, 39658, + 41414, 43247, 45162, 47161, 49249, 51430, 53707, 56085, 58568, + 61161, 63869, 66696, 69649, 72733, 75953, 79316, 82827, 86494, + 90324, 94323, 98499, 102860, 107414, 112169, 117135, 122322, 127737, + 133393, 139298, 145466, 151906, 158631, 165655, 172989, 180648, 188646, + 196998, 205719, 214827, 224339, 234271, 244643, 255474, 266785, 278597, + 290931, 303812, 317263, 331309, 345977, 361295, 377291, 393995, 411439, + 429655, 448677, 468542, 489286, 510949, 533570, 557193, 581862, 607624, + 634525, 662618, 691955, 722590, 754582, 787990, 822878, 859310, 897355, + 937084, 978572]) * u.Hz, + } + + + array = np.zeros((128, 7322)) + parse_path_moc.return_value = [(array, meta)] + + file = Path("solo_L3_rpw-tnr-surv-flux_20240323_V01.cdf") + spec = Spectrogram(file) + + assert isinstance(spec, RPWSpectrogram) + assert spec.observatory == "SOLO" + assert spec.instrument == "RPW" + assert spec.detector == "RPW-TNR" + assert spec.start_time.datetime == datetime(2024, 3, 23, 0, 0, 0, 0) + assert spec.end_time.datetime == datetime(2024, 3, 24, 0, 0, 0, 0) + assert spec.wavelength.min == 3.992000102996826 * u.kHz + assert spec.wavelength.max == 978.572021484375 * u.kHz + + +@mock.patch("radiospectra.spectrogram.spectrogram_factory.parse_path") +def test_solo_rpw_hfr(parse_path_moc): + start_time = Time("2024-03-23 00:00:00.000") + end_time = Time("2024-03-24 00:00:00.000") + + meta = { + "cdf_meta": { + "Project": "SOLO>Solar Orbiter", + "Source_name": "SOLO>Solar Orbiter", + "Descriptor": "RPW-HFR-SURV-FLUX>Radio and Plasma Waves,High Frequency Receiver, Survey mode, spectral data in physical units", + "Data_type": "L3>Level 3 Data", + "Data_version": "01", + "Logical_file_id": "solo_L3_rpw-hfr-surv-flux_20240323_V01", + "Instrument": "RPW", + }, + "detector": "RPW-HFR", + "instrument": "RPW", + "observatory": "SOLO", + "start_time": start_time, + "end_time": end_time, + "wavelength": a.Wavelength(425 * u.kHz, 16325 * u.kHz), + "times": start_time + np.linspace(0, (end_time - start_time).to_value("s"), 16499) * u.s, + "freqs": np.array([425000, 525000, 625000, 675000, 775000, 875000, 975000, 1025000, + 1225000, 1375000, 1475000, 1725000, 1825000, 2075000, 2425000, + 2675000, 3225000, 3325000, 3525000, 3825000, 3925000, 4125000, + 4525000, 4875000, 5225000, 5475000, 5825000, 6175000, 6525000, + 6875000, 7375000, 7625000, 7975000, 8225000, 8575000, 9175000, + 10025000, 10125000, 11025000, 11375000, 12225000, 12475000, + 13375000, 13725000, 14375000, 14925000, 15275000, 15625000, + 16075000, 16325000]) * u.Hz, + } + + + array = np.zeros((50, 16499)) + parse_path_moc.return_value = [(array, meta)] + + file = Path("solo_L3_rpw-hfr-surv-flux_20240323_V01.cdf") + spec = Spectrogram(file) + + assert isinstance(spec, RPWSpectrogram) + assert spec.observatory == "SOLO" + assert spec.instrument == "RPW" + assert spec.detector == "RPW-HFR" + assert spec.start_time.datetime == datetime(2024, 3, 23, 0, 0, 0, 0) + assert spec.end_time.datetime == datetime(2024, 3, 24, 0, 0, 0, 0) + assert spec.wavelength.min == 425 * u.kHz + assert spec.wavelength.max == 16325 * u.kHz diff --git a/radiospectra/spectrogram/spectrogram_factory.py b/radiospectra/spectrogram/spectrogram_factory.py index 09a702f..f3d8b0b 100644 --- a/radiospectra/spectrogram/spectrogram_factory.py +++ b/radiospectra/spectrogram/spectrogram_factory.py @@ -20,6 +20,7 @@ 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, @@ -440,129 +441,154 @@ def _read_cdf(file): "freqs": freqs, } return data, meta - elif "SOLO" in cdf_globals.get("Project", "")[0]: - if "RPW-HFR-SURV" not in cdf_globals.get("Descriptor", "")[0]: - raise ValueError( - f"Currently radiospectra supports Level 2 HFR survey data the file" - f'{file.name} is {cdf_globals.get("Descriptor", "")}' + 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"] ) - - # 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 = { + 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": "RPW-AGC1", + "detector": detector, "instrument": "RPW", "observatory": "SOLO", "start_time": times[0], "end_time": times[-1], - "wavelength": a.Wavelength(hfr_frequency.min(), hfr_frequency.max()), + "wavelength": a.Wavelength(freqs.min(), freqs.max()), "times": times, - "freqs": hfr_frequency, + "freqs": freqs, } - 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, + 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", } - res.append((specs[1].T, meta2)) - return res + + # 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 @staticmethod def _read_fits(file):