Add support for SOLO L3

changelog/148.feature.rst

@@ -0,0 +1 @@
+Added support for Solar Orbiter RPW Level 3 TNR and HFR survey data products.

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

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):