diff --git a/imap_processing/cdf/config/imap_mag_global_cdf_attrs.yaml b/imap_processing/cdf/config/imap_mag_global_cdf_attrs.yaml index 141820195..c3c93eccf 100644 --- a/imap_processing/cdf/config/imap_mag_global_cdf_attrs.yaml +++ b/imap_processing/cdf/config/imap_mag_global_cdf_attrs.yaml @@ -13,14 +13,14 @@ default: &default imap_mag_l1a_norm-raw: <<: *default - Data_type: L1A_raw-norm>Level-1A raw normal rate + Data_type: L1A_norm-raw>Level-1A raw normal rate Logical_source: imap_mag_l1a_norm-raw Logical_source_description: IMAP Mission MAG Normal Rate Instrument Level-1A Data. Data_level: L1A imap_mag_l1a_burst-raw: <<: *default - Data_type: L1A_raw-burst>Level-1A raw burst rate + Data_type: L1A_burst-raw>Level-1A raw burst rate Logical_source: imap_mag_l1a_burst-raw Logical_source_description: IMAP Mission MAG Burst Rate Instrument Level-1A Data. Data_level: L1A diff --git a/imap_processing/cdf/config/imap_mag_l1_variable_attrs.yaml b/imap_processing/cdf/config/imap_mag_l1_variable_attrs.yaml index 0e8c2543c..dd7ed75a6 100644 --- a/imap_processing/cdf/config/imap_mag_l1_variable_attrs.yaml +++ b/imap_processing/cdf/config/imap_mag_l1_variable_attrs.yaml @@ -1,3 +1,19 @@ +epoch: + CATDESC: Time, number of nanoseconds since J2000 with leap seconds included + FIELDNAM: epoch + LABLAXIS: epoch + FILLVAL: -9223372036854775808 + FORMAT: " " # Supposedly not required, fails in xarray_to_cdf + VALIDMIN: -9223372036854775808 + VALIDMAX: 9223372036854775807 + UNITS: ns + VAR_TYPE: support_data + SCALETYP: linear + MONOTON: INCREASE + TIME_BASE: J2000 + TIME_SCALE: Terrestrial Time + REFERENCE_POSITION: Rotating Earth Geoid + default_attrs: &default # Assumed values for all variable attrs unless overwritten DEPEND_0: epoch @@ -7,25 +23,34 @@ default_attrs: &default VALIDMIN: -9223372036854775808 VALIDMAX: 9223372036854775807 VAR_TYPE: data + UNITS: ' ' + +default_coords: &default_coords + # Assumed values for all coordinate attrs unless overwritten + FORMAT: F2.4 # Float with 4 digits + VAR_TYPE: data + DISPLAY_TYPE: time_series + FILLVAL: -9223372036854775808 + VALIDMIN: -9223372036854775808 + VALIDMAX: 9223372036854775807 + UNITS: counts -mag_raw_vector_attrs: - <<: *default +raw_vector_attrs: + <<: *default_coords CATDESC: Raw unprocessed magnetic field vector data in bytes DEPEND_1: direction FIELDNAM: Magnetic Field Vector LABLAXIS: Raw binary magnetic field vector data FORMAT: I3 - UNITS: counts vector_attrs: - <<: *default + <<: *default_coords CATDESC: Magnetic field vector with x y z and sensor range varying by time DEPEND_1: direction FIELDNAM: Magnetic Field Vector LABLAXIS: Magnetic field vector data FILLVAL: 9223372036854775807 FORMAT: I3 - UNITS: counts mag_support_attrs: &support_default <<: *default @@ -44,13 +69,13 @@ mag_flag_attrs: FORMAT: I1 raw_direction_attrs: - <<: *default + <<: *default_coords CATDESC: Raw magnetic field vector binary length FIELDNAM: Raw magnetic field vector binary length LABLAXIS: Magnetic field vector directions direction_attrs: - <<: *default + <<: *default_coords CATDESC: Magnetic field vector FIELDNAM: \[xyz\] magnetic field vector FORMAT: I3 @@ -154,4 +179,39 @@ VECTORS: CATDESC: Raw binary value of MAG Science vectors before processing FIELDNAM: Raw vector binary - +# CCSDS attributes for MAG L1A RAW + +VERSION: + <<: *metadata_default + CATDESC: CCSDS Packet version number + FIELDNAM: CCSDS Packet version number + +TYPE: + <<: *metadata_default + CATDESC: CCSDS Packet type + FIELDNAM: CCSDS Packet type + +SEC_HDR_FLG: + <<: *metadata_default + CATDESC: CCSDS Packet Secondary header flag + FIELDNAM: CCSDS Packet Secondary header flag + +PKT_APID: + <<: *metadata_default + CATDESC: CCSDS Packet Application Process ID + FIELDNAM: CCSDS Packet Application Process ID + +SEQ_FLGS: + <<: *metadata_default + CATDESC: CCSDS Packet Grouping flags + FIELDNAM: CCSDS Packet Grouping flags + +SRC_SEQ_CTR: + <<: *metadata_default + CATDESC: CCSDS Packet Source Sequence Counter + FIELDNAM: CCSDS Packet Source Sequence Counter + +PKT_LEN: + <<: *metadata_default + CATDESC: CCSDS Packet Length + FIELDNAM: CCSDS Packet Length \ No newline at end of file diff --git a/imap_processing/cli.py b/imap_processing/cli.py index 80b1989fd..30024fba4 100644 --- a/imap_processing/cli.py +++ b/imap_processing/cli.py @@ -645,14 +645,13 @@ def do_processing(self, dependencies: list) -> list[xr.Dataset]: class Mag(ProcessInstrument): """Process MAG.""" - # TODO: fixed in later PR - def do_processing(self, dependencies: list) -> list[xr.Dataset]: + def do_processing(self, dependencies: list[Path]) -> list[xr.Dataset]: """ Perform MAG specific processing. Parameters ---------- - dependencies : list + dependencies : list[Path] List of dependencies to process. Returns diff --git a/imap_processing/mag/__init__.py b/imap_processing/mag/__init__.py index 27d3f0217..2828fe266 100644 --- a/imap_processing/mag/__init__.py +++ b/imap_processing/mag/__init__.py @@ -1 +1,2 @@ __version__ = "01" +cdf_format_version = "v0.1" diff --git a/imap_processing/mag/constants.py b/imap_processing/mag/constants.py new file mode 100644 index 000000000..34e00e68b --- /dev/null +++ b/imap_processing/mag/constants.py @@ -0,0 +1,57 @@ +"""Collection of constant types or values for MAG.""" + +from enum import Enum + + +class DataMode(Enum): + """ + Enum for MAG data modes: burst and normal (BURST + NORM). + + Attributes + ---------- + BURST: str + Burst data mode - higher frequency data + NORM: str + Normal data mode - lower frequency data (downsampled from burst) + """ + + BURST = "BURST" + NORM = "NORM" + + +class Sensor(Enum): + """ + Enum for MAG sensors: raw, MAGo, and MAGi (RAW, MAGO, MAGI). + + Attributes + ---------- + MAGO : str + MAGo sensor - for the outboard sensor. This is nominally expected to be the + primary sensor. + MAGI : str + MAGi sensor - for the inboard sensor. + RAW : str + RAW data - contains both sensors. Here, the vectors are unprocessed. + """ + + MAGO = "MAGO" + MAGI = "MAGI" + RAW = "RAW" + + +class PrimarySensor(Enum): + """ + Enum for primary sensor: MAGo and MAGi (MAGO, MAGI). + + This corresponds to the PRI_SENS field in the MAG Level 0 data. + + Attributes + ---------- + MAGO : int + Primary sensor is MAGo. + MAGI : int + Primary sensor is MAGi. + """ + + MAGO = 0 + MAGI = 1 diff --git a/imap_processing/mag/l0/decom_mag.py b/imap_processing/mag/l0/decom_mag.py index 4e0e726ce..4ea4d8e8f 100644 --- a/imap_processing/mag/l0/decom_mag.py +++ b/imap_processing/mag/l0/decom_mag.py @@ -13,9 +13,9 @@ from imap_processing import imap_module_directory from imap_processing.ccsds.ccsds_data import CcsdsData -from imap_processing.cdf.global_attrs import ConstantCoordinates +from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.cdf.utils import met_to_j2000ns -from imap_processing.mag import mag_cdf_attrs +from imap_processing.mag.constants import DataMode from imap_processing.mag.l0.mag_l0_data import MagL0, Mode logger = logging.getLogger(__name__) @@ -67,7 +67,9 @@ def decom_packets(packet_file_path: str | Path) -> dict[str, list[MagL0]]: return {"norm": norm_data, "burst": burst_data} -def generate_dataset(l0_data: list[MagL0], dataset_attrs: dict) -> xr.Dataset: +def generate_dataset( + l0_data: list[MagL0], mode: DataMode, attribute_manager: ImapCdfAttributes +) -> xr.Dataset: """ Generate a CDF dataset from the sorted raw L0 MAG data. @@ -76,8 +78,11 @@ def generate_dataset(l0_data: list[MagL0], dataset_attrs: dict) -> xr.Dataset: l0_data : list[MagL0] List of sorted L0 MAG data. - dataset_attrs : dict - Global attributes for the dataset. + mode : DataMode + The mode of the CDF file - burst or norm. + + attribute_manager : ImapCdfAttributes + Attribute manager for the dataset, including all MAG L1A attributes. Returns ------- @@ -115,33 +120,39 @@ def generate_dataset(l0_data: list[MagL0], dataset_attrs: dict) -> xr.Dataset: for key, value in dataclasses.asdict(datapoint).items(): if key not in ("ccsds_header", "VECTORS", "SHCOARSE"): support_data[key].append(value) + if key == "ccsds_header": + for ccsds_key, ccsds_value in value.items(): + support_data[ccsds_key].append(ccsds_value) # Used in L1A vectors direction = xr.DataArray( np.arange(vector_data.shape[1]), name="direction", dims=["direction"], - attrs=mag_cdf_attrs.raw_direction_attrs.output(), + attrs=attribute_manager.get_variable_attributes("raw_direction_attrs"), ) + # TODO: Epoch here refers to the start of the sample. Confirm that this is # what mag is expecting, and if it is, CATDESC needs to be updated. epoch_time = xr.DataArray( shcoarse_data, name="epoch", dims=["epoch"], - attrs=ConstantCoordinates.EPOCH, + attrs=attribute_manager.get_variable_attributes("epoch"), ) # TODO: raw vectors units raw_vectors = xr.DataArray( vector_data, name="raw_vectors", dims=["epoch", "direction"], - attrs=mag_cdf_attrs.mag_raw_vector_attrs.output(), + attrs=attribute_manager.get_variable_attributes("raw_vector_attrs"), ) + logical_id = f"imap_mag_l1a_{mode.value.lower()}-raw" + output = xr.Dataset( coords={"epoch": epoch_time, "direction": direction}, - attrs=dataset_attrs, + attrs=attribute_manager.get_global_attributes(logical_id), ) output["raw_vectors"] = raw_vectors @@ -153,14 +164,7 @@ def generate_dataset(l0_data: list[MagL0], dataset_attrs: dict) -> xr.Dataset: value, name=key.lower(), dims=["epoch"], - attrs=dataclasses.replace( - mag_cdf_attrs.mag_support_attrs, - catdesc=mag_cdf_attrs.catdesc_fieldname_l0[key][0], - fieldname=mag_cdf_attrs.catdesc_fieldname_l0[key][1], - # TODO: label_axis should be as close to 6 letters as possible - label_axis=key, - display_type="no_plot", - ).output(), + attrs=attribute_manager.get_variable_attributes(key), ) return output diff --git a/imap_processing/mag/l1a/mag_l1a.py b/imap_processing/mag/l1a/mag_l1a.py index c5e7d10bf..1838bf203 100644 --- a/imap_processing/mag/l1a/mag_l1a.py +++ b/imap_processing/mag/l1a/mag_l1a.py @@ -1,22 +1,27 @@ """Methods for decomming packets, processing to level 1A, and writing CDFs for MAG.""" +import dataclasses import logging from pathlib import Path import numpy as np import xarray as xr -from imap_processing.cdf.global_attrs import ConstantCoordinates +from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.cdf.utils import J2000_EPOCH, met_to_j2000ns -from imap_processing.mag import mag_cdf_attrs +from imap_processing.mag.constants import DataMode, PrimarySensor from imap_processing.mag.l0 import decom_mag from imap_processing.mag.l0.mag_l0_data import MagL0 -from imap_processing.mag.l1a.mag_l1a_data import MagL1a, TimeTuple +from imap_processing.mag.l1a.mag_l1a_data import ( + MagL1a, + MagL1aPacketProperties, + TimeTuple, +) logger = logging.getLogger(__name__) -def mag_l1a(packet_filepath: Path, data_version: str) -> list[Path]: +def mag_l1a(packet_filepath: Path, data_version: str) -> list[xr.Dataset]: """ Will process MAG L0 data into L1A CDF files at cdf_filepath. @@ -29,7 +34,7 @@ def mag_l1a(packet_filepath: Path, data_version: str) -> list[Path]: Returns ------- - generated_files : list[pathlib.Path] + generated_files : list[xarray.Dataset] A list of generated filenames. """ packets = decom_mag.decom_packets(packet_filepath) @@ -37,31 +42,31 @@ def mag_l1a(packet_filepath: Path, data_version: str) -> list[Path]: norm_data = packets["norm"] burst_data = packets["burst"] - generated_files = process_and_write_data( - norm_data, - mag_cdf_attrs.mag_l1a_norm_raw_attrs.output(), - mag_cdf_attrs.mag_l1a_norm_mago_attrs.output(), - mag_cdf_attrs.mag_l1a_norm_magi_attrs.output(), - data_version, - ) - generated_files += process_and_write_data( - burst_data, - mag_cdf_attrs.mag_l1a_burst_raw_attrs.output(), - mag_cdf_attrs.mag_l1a_burst_mago_attrs.output(), - mag_cdf_attrs.mag_l1a_burst_magi_attrs.output(), - data_version, + input_files = [packet_filepath.name] + + # Create attribute manager and add MAG L1A attributes and global variables + attribute_manager = ImapCdfAttributes() + attribute_manager.add_instrument_global_attrs("mag") + attribute_manager.add_instrument_variable_attrs("mag", "l1") + + attribute_manager.add_global_attribute("Data_version", data_version) + attribute_manager.add_global_attribute("Input_files", str(input_files)) + attribute_manager.add_global_attribute( + "Generation_date", + np.datetime64( + "now", + ).astype(str), ) - return generated_files + generated_datasets = create_l1a(norm_data, DataMode.NORM, attribute_manager) + generated_datasets += create_l1a(burst_data, DataMode.BURST, attribute_manager) + + return generated_datasets -def process_and_write_data( - packet_data: list[MagL0], - raw_attrs: dict, - mago_attrs: dict, - magi_attrs: dict, - data_version: str, -) -> list[Path]: +def create_l1a( + packet_data: list[MagL0], data_mode: DataMode, attribute_manager: ImapCdfAttributes +) -> list[xr.Dataset]: """ Will process MAG L0 data into L1A, then create and write out CDF files. @@ -71,40 +76,41 @@ def process_and_write_data( ---------- packet_data : list[MagL0] List of MagL0 packets to process, containing primary and secondary sensor data. - raw_attrs : dict - Attributes for MagL1A raw CDF files. - mago_attrs : dict - Attributes for MagL1A MAGo CDF files. - magi_attrs : dict - Attributes for MagL1A MAGi CDF files. - data_version : str - Data version. + + data_mode : DataMode + Enum for distinguishing between norm and burst mode data. + + attribute_manager : ImapCdfAttributes + Attribute manager for CDF files for MAG L1A. Returns ------- - generated_files : list[pathlib.Path] + generated_files : list[xarray.Dataset] A list of generated filenames. """ if not packet_data: return [] - # TODO: Rework attrs to be better - raw_attrs["Data_version"] = data_version - magi_attrs["Data_version"] = data_version - mago_attrs["Data_version"] = data_version - - mag_raw = decom_mag.generate_dataset(packet_data, raw_attrs) + mag_raw = decom_mag.generate_dataset(packet_data, data_mode, attribute_manager) generated_datasets = [mag_raw] l1a = process_packets(packet_data) + # TODO: Rearrange generate_dataset to combine these two for loops + # Split into MAGo and MAGi for _, mago in l1a["mago"].items(): - norm_mago_output = generate_dataset(mago, mago_attrs) + logical_file_id = f"imap_mag_l1a_{data_mode.value.lower()}-mago" + norm_mago_output = generate_dataset(mago, logical_file_id, attribute_manager) generated_datasets.append(norm_mago_output) for _, magi in l1a["magi"].items(): - norm_magi_output = generate_dataset(magi, magi_attrs) + logical_file_id = f"imap_mag_l1a_{data_mode.value.lower()}-magi" + norm_magi_output = generate_dataset( + magi, + logical_file_id, + attribute_manager, + ) generated_datasets.append(norm_magi_output) return generated_datasets @@ -116,6 +122,9 @@ def process_packets( """ Given a list of MagL0 packets, process them into MagO and MagI L1A data classes. + This splits the MagL0 packets into MagO and MagI data, returning a dictionary with + keys "mago" and "magi." + Parameters ---------- mag_l0_list : list[MagL0] @@ -139,6 +148,8 @@ def process_packets( primary_start_time = TimeTuple(mag_l0.PRI_COARSETM, mag_l0.PRI_FNTM) secondary_start_time = TimeTuple(mag_l0.SEC_COARSETM, mag_l0.SEC_FNTM) + mago_is_primary = mag_l0.PRI_SENS == PrimarySensor.MAGO.value + primary_day = ( J2000_EPOCH + met_to_j2000ns(primary_start_time.to_seconds()).astype("timedelta64[ns]") @@ -150,88 +161,98 @@ def process_packets( ) ).astype("datetime64[D]") - # seconds of data in this packet is the SUBTYPE plus 1 - seconds_per_packet = mag_l0.PUS_SSUBTYPE + 1 + primary_packet_data = MagL1aPacketProperties( + mag_l0.SHCOARSE, + primary_start_time, + mag_l0.PRI_VECSEC, + mag_l0.PUS_SSUBTYPE, + mag_l0.ccsds_header.SRC_SEQ_CTR, + mag_l0.COMPRESSION, + mago_is_primary, + ) + secondary_packet_data = dataclasses.replace( + primary_packet_data, + start_time=secondary_start_time, + vectors_per_second=mag_l0.SEC_VECSEC, + pus_ssubtype=mag_l0.PUS_SSUBTYPE, + ) # now we know the number of secs of data in the packet, and the data rates of # each sensor, we can calculate how much data is in this packet and where the # byte boundaries are. - # VECSEC is already decoded in mag_l0 - total_primary_vectors = seconds_per_packet * mag_l0.PRI_VECSEC - total_secondary_vectors = seconds_per_packet * mag_l0.SEC_VECSEC - - # The raw vectors are of type int8, but the output vectors should be at least - # int16. primary_vectors, secondary_vectors = MagL1a.process_vector_data( mag_l0.VECTORS.astype(dtype=np.int32), # type: ignore[union-attr] # TODO Maybe Change, Item "str" of "Union[Any, str]" # has no attribute "astype" # this is because mypy expects both to have the attributes - total_primary_vectors, - total_secondary_vectors, + primary_packet_data.total_vectors, + secondary_packet_data.total_vectors, ) primary_timestamped_vectors = MagL1a.calculate_vector_time( - primary_vectors, mag_l0.PRI_VECSEC, primary_start_time + primary_vectors, + primary_packet_data.vectors_per_second, + primary_packet_data.start_time, ) secondary_timestamped_vectors = MagL1a.calculate_vector_time( - secondary_vectors, mag_l0.SEC_VECSEC, secondary_start_time + secondary_vectors, + secondary_packet_data.vectors_per_second, + secondary_packet_data.start_time, ) # Sort primary and secondary into MAGo and MAGi by 24 hour chunks - # TODO: Individual vectors should be sorted by day, not the whole packet - mago_is_primary = mag_l0.PRI_SENS == 0 - mago_day = primary_day if mago_is_primary else secondary_day magi_day = primary_day if not mago_is_primary else secondary_day if mago_day not in mago: mago[mago_day] = MagL1a( True, - bool(mag_l0.MAGO_ACT), + mag_l0.MAGO_ACT, mag_l0.SHCOARSE, primary_timestamped_vectors if mago_is_primary else secondary_timestamped_vectors, + primary_packet_data if mago_is_primary else secondary_packet_data, ) else: - mago[mago_day].vectors = np.concatenate( - [ - mago[mago_day].vectors, - ( - primary_timestamped_vectors - if mago_is_primary - else secondary_timestamped_vectors - ), - ] + mago[mago_day].append_vectors( + ( + primary_timestamped_vectors + if mago_is_primary + else secondary_timestamped_vectors + ), + primary_packet_data if mago_is_primary else secondary_packet_data, ) if magi_day not in magi: magi[magi_day] = MagL1a( False, - bool(mag_l0.MAGI_ACT), + mag_l0.MAGI_ACT, mag_l0.SHCOARSE, primary_timestamped_vectors if not mago_is_primary else secondary_timestamped_vectors, + primary_packet_data if not mago_is_primary else secondary_packet_data, ) else: - magi[magi_day].vectors = np.concatenate( - [ - magi[magi_day].vectors, - ( - primary_timestamped_vectors - if not mago_is_primary - else secondary_timestamped_vectors - ), - ] + magi[magi_day].append_vectors( + ( + primary_timestamped_vectors + if not mago_is_primary + else secondary_timestamped_vectors + ), + primary_packet_data if not mago_is_primary else secondary_packet_data, ) return {"mago": mago, "magi": magi} -def generate_dataset(single_file_l1a: MagL1a, dataset_attrs: dict) -> xr.Dataset: +def generate_dataset( + single_file_l1a: MagL1a, + logical_file_id: str, + attribute_manager: ImapCdfAttributes, +) -> xr.Dataset: """ Generate a Xarray dataset for L1A data to output to CDF files. @@ -245,14 +266,21 @@ def generate_dataset(single_file_l1a: MagL1a, dataset_attrs: dict) -> xr.Dataset ---------- single_file_l1a : MagL1a L1A data covering one day to process into a xarray dataset. - dataset_attrs : dict - Global attributes for the dataset, as created by mag_attrs. + logical_file_id : str + Indicates which sensor (MagO or MAGi) and mode (burst or norm) the data is from. + This is used to retrieve the global attributes from attribute_manager. + attribute_manager : ImapCdfAttributes + Attributes for the dataset, as created by ImapCdfAttributes. Returns ------- dataset : xarray.Dataset One xarray dataset with proper CDF attributes and shape containing MAG L1A data. """ + # TODO: add: + # gaps_in_data global attr + # magl1avectordefinition data + # TODO: Just leave time in datetime64 type with vector as dtype object to avoid this # Get the timestamp from the end of the vector time_data = single_file_l1a.vectors[:, 4].astype( @@ -263,7 +291,7 @@ def generate_dataset(single_file_l1a: MagL1a, dataset_attrs: dict) -> xr.Dataset np.arange(4), name="direction", dims=["direction"], - attrs=mag_cdf_attrs.direction_attrs.output(), + attrs=attribute_manager.get_variable_attributes("direction_attrs"), ) # TODO: Epoch here refers to the start of the sample. Confirm that this is @@ -272,19 +300,19 @@ def generate_dataset(single_file_l1a: MagL1a, dataset_attrs: dict) -> xr.Dataset time_data, name="epoch", dims=["epoch"], - attrs=ConstantCoordinates.EPOCH, + attrs=attribute_manager.get_variable_attributes("epoch"), ) vectors = xr.DataArray( single_file_l1a.vectors[:, :4], name="vectors", dims=["epoch", "direction"], - attrs=mag_cdf_attrs.vector_attrs.output(), + attrs=attribute_manager.get_variable_attributes("vector_attrs"), ) output = xr.Dataset( coords={"epoch": epoch_time, "direction": direction}, - attrs=dataset_attrs, + attrs=attribute_manager.get_global_attributes(logical_file_id), ) output["vectors"] = vectors diff --git a/imap_processing/mag/l1a/mag_l1a_data.py b/imap_processing/mag/l1a/mag_l1a_data.py index 40f4f509d..10162293d 100644 --- a/imap_processing/mag/l1a/mag_l1a_data.py +++ b/imap_processing/mag/l1a/mag_l1a_data.py @@ -1,11 +1,13 @@ """Data classes for storing and processing MAG Level 1A data.""" -from dataclasses import dataclass +from __future__ import annotations + +from dataclasses import InitVar, dataclass, field from math import floor import numpy as np -from imap_processing.cdf.utils import met_to_j2000ns +from imap_processing.cdf.utils import J2000_EPOCH, met_to_j2000ns MAX_FINE_TIME = 65535 # maximum 16 bit unsigned int @@ -33,8 +35,7 @@ class TimeTuple: coarse_time: int fine_time: int - def __add__(self, seconds: float): # type: ignore[no-untyped-def] - # ruff is saying TimeTuple is undefined for this usage. + def __add__(self, seconds: float) -> TimeTuple: """ Add a number of seconds to the time tuple. @@ -72,13 +73,78 @@ def to_seconds(self) -> float: return self.coarse_time + self.fine_time / MAX_FINE_TIME +@dataclass +class MagL1aPacketProperties: + """ + Data class with Mag L1A per-packet data. + + This contains per-packet variations in L1a data that is passed into CDF + files. Since each L1a file contains multiple packets, the variables in this + class vary by time in the end CDF file. + + seconds_per_packet, and total_vectors are calculated from pus_ssubtype and + vecsec values, which are only passed in to the init method and cannot be + accessed from an instance as they are InitVars. + + To use the class, pass in pus_ssubtype and either PRI_VECSEC or SEC_VECSEC, + then you can access seconds_per_packet and total_vectors. + + Attributes + ---------- + shcoarse : int + Mission elapsed time for the packet + start_time : TimeTuple + Start time of the packet + vectors_per_second : int + Number of vectors per second + pus_ssubtype : int + PUS Service Subtype - used to calculate seconds_per_packet. This is an InitVar, + meaning it is only used when creating the class and cannot be accessed from an + instance of the class - instead seconds_per_packet should be used. + src_seq_ctr : int + Sequence counter from the ccsds header + compression : int + Science Data Compression Flag from level 0 + mago_is_primary : int + 1 if mago is designated the primary sensor, otherwise 0 + seconds_per_packet : int + Number of seconds of data in this packet - calculated as pus_ssubtype + 1 + total_vectors : int + Total number of vectors in this packet - calculated as + seconds_per_packet * vecsec + """ + + shcoarse: int + start_time: TimeTuple + vectors_per_second: int + pus_ssubtype: InitVar[int] + src_seq_ctr: int # From ccsds header + compression: int + mago_is_primary: int + seconds_per_packet: int = field(init=False) + total_vectors: int = field(init=False) + + def __post_init__(self, pus_ssubtype: int) -> None: + """ + Calculate seconds_per_packet and total_vectors. + + Parameters + ---------- + pus_ssubtype : int + PUS Service Subtype, used to determine the seconds of data in the packet. + """ + # seconds of data in this packet is the SUBTYPE plus 1 + self.seconds_per_packet = pus_ssubtype + 1 + + # VECSEC is already decoded in mag_l0 + self.total_vectors = self.seconds_per_packet * self.vectors_per_second + + @dataclass class MagL1a: """ Data class for MAG Level 1A data. - TODO: This object should end up having 1 days worth of data. - One MAG L1A object corresponds to part of one MAG L0 packet, which corresponds to one packet of data from the MAG instrument. Each L0 packet consists of data from two sensors, MAGO (outboard) and MAGI (inboard). One of these sensors is designated @@ -92,30 +158,93 @@ class MagL1a: Attributes ---------- is_mago : bool - True if the data is from MagO, False if data is from MagI. - active : bool - True if the sensor is active. - SHCOARSE : int - Mission elapsed time. + True if the data is from MagO, False if data is from MagI + is_active : int + 1 if the sensor is active, 0 if not + shcoarse : int + Mission elapsed time for the first packet, the start time for the whole day vectors : numpy.ndarray List of magnetic vector samples, starting at start_time. [x, y, z, range, time], - where time is np.datetime64[ns]. + where time is numpy.datetime64[ns] + starting_packet : InitVar[MagL1aPacketProperties] + The packet properties for the first packet in the day. As an InitVar, this + cannot be accessed from an instance of the class. Instead, packet_definitions + should be used. + packet_definitions : dict[numpy.datetime64, MagL1aPacketProperties] + Dictionary of packet properties for each packet in the day. The key is the start + time of the packet, and the value is a dataclass of packet properties. + most_recent_sequence : int + Sequence number of the most recent packet added to the object + missing_sequences : list[int] + List of missing sequence numbers in the day + start_time : numpy.datetime64 + Start time of the day, in ns since J2000 epoch Methods ------- - calculate_vector_time(vectors, vectors_per_second, start_time) - process_vector_data(vector_data, primary_count, secondary_count) + append_vectors() + calculate_vector_time() + process_vector_data() """ is_mago: bool - active: bool - SHCOARSE: int - vectors: np.ndarray + is_active: int + shcoarse: int + vectors: np.array + starting_packet: InitVar[MagL1aPacketProperties] + packet_definitions: dict[np.datetime64, MagL1aPacketProperties] = field(init=False) + most_recent_sequence: int = field(init=False) + missing_sequences: list[int] = field(default_factory=list) + start_time: np.datetime64 = field(init=False) + + def __post_init__(self, starting_packet: MagL1aPacketProperties) -> None: + """ + Initialize the packet_definition dictionary and most_recent_sequence. + + Parameters + ---------- + starting_packet : MagL1aPacketProperties + The packet properties for the first packet in the day, including start time. + """ + # TODO should this be from starting_packet + self.start_time = (J2000_EPOCH + met_to_j2000ns(self.shcoarse)).astype( + "datetime64[D]" + ) + self.packet_definitions = {self.start_time: starting_packet} + # most_recent_sequence is the sequence number of the packet used to initialize + # the object + self.most_recent_sequence = starting_packet.src_seq_ctr + + def append_vectors( + self, additional_vectors: np.array, packet_properties: MagL1aPacketProperties + ) -> None: + """ + Append additional vectors to the current vectors array. + + Parameters + ---------- + additional_vectors : numpy.array + New vectors to append. + packet_properties : MagL1aPacketProperties + Additional vector definition to add to the l0_packets dictionary. + """ + vector_sequence = packet_properties.src_seq_ctr + + self.vectors = np.concatenate([self.vectors, additional_vectors]) + self.packet_definitions[self.start_time] = packet_properties + + # Every additional packet should be the next one in the sequence, if not, add + # the missing sequence(s) to the gap data + if not self.most_recent_sequence + 1 == vector_sequence: + self.missing_sequences += list( + range(self.most_recent_sequence + 1, vector_sequence) + ) + self.most_recent_sequence = vector_sequence @staticmethod def calculate_vector_time( - vectors: np.ndarray, vectors_per_second: int, start_time: TimeTuple - ) -> np.ndarray: + vectors: np.ndarray, vectors_per_sec: int, start_time: TimeTuple + ) -> np.array: """ Add timestamps to the vector list, turning the shape from (n, 4) to (n, 5). @@ -124,12 +253,12 @@ def calculate_vector_time( Parameters ---------- - vectors : numpy.ndarray + vectors : numpy.array List of magnetic vector samples, starting at start_time. Shape of (n, 4). - vectors_per_second : int + vectors_per_sec : int Number of vectors per second. start_time : TimeTuple - The coarse and fine time for the sensor. + Start time of the vectors, the timestamp of the first vector. Returns ------- @@ -137,9 +266,8 @@ def calculate_vector_time( Vectors with timestamps added in seconds, calculated from cdf.utils.met_to_j2000ns. """ - # TODO: Move timestamps to J2000 - timedelta = np.timedelta64(int(1 / vectors_per_second * 1e9), "ns") - + timedelta = np.timedelta64(int(1 / vectors_per_sec * 1e9), "ns") + # TODO: validate that start_time from SHCOARSE is precise enough start_time_ns = met_to_j2000ns(start_time.to_seconds()) # Calculate time skips for each vector in ns diff --git a/imap_processing/mag/mag_cdf_attrs.py b/imap_processing/mag/mag_cdf_attrs.py deleted file mode 100644 index c29d04a0e..000000000 --- a/imap_processing/mag/mag_cdf_attrs.py +++ /dev/null @@ -1,249 +0,0 @@ -"""Shared attribute values for MAG CDF files.""" - -from imap_processing.cdf.defaults import GlobalConstants -from imap_processing.cdf.global_attrs import ( - AttrBase, - GlobalDataLevelAttrs, - GlobalInstrumentAttrs, - ScienceAttrs, -) -from imap_processing.mag import __version__ - -text = ( - "The IMAP magnetometer (MAG) consists of a pair of identical magnetometers " - "which each measure the magnetic field in three directions in the vicinity of " - "the spacecraft. " - "MAG will contribute to our understanding of the acceleration and transport " - "of charged particles in the heliosphere. " - "MAG design and assembly is led by Imperial College, London. See " - "https://imap.princeton.edu/instruments/mag for more details." -) - -mag_base = GlobalInstrumentAttrs( - __version__, - "MAG>Magnetometer", - text, - "Magnetic Fields (space)", -) - -mag_l1a_norm_raw_attrs = GlobalDataLevelAttrs( - "L1A_raw-norm>Level-1A raw normal rate", - logical_source="imap_mag_l1a_norm-raw", - logical_source_desc="IMAP Mission MAG Normal Rate Instrument Level-1A Data.", - instrument_base=mag_base, -) - -mag_l1a_burst_raw_attrs = GlobalDataLevelAttrs( - "L1A_raw-burst>Level-1A raw burst rate", - logical_source="imap_mag_l1a_burst-raw", - logical_source_desc="IMAP Mission MAG Burst Rate Instrument Level-1A Data.", - instrument_base=mag_base, -) - -mag_l1a_norm_mago_attrs = GlobalDataLevelAttrs( - "L1A_raw-norm-mago>Level 1A MAGo normal rate", - logical_source="imap_mag_l1a_norm-mago", - logical_source_desc="IMAP Mission MAGo Normal Rate Instrument Level-1A Data.", - instrument_base=mag_base, -) - -mag_l1a_norm_magi_attrs = GlobalDataLevelAttrs( - "L1A_raw-norm-magi>Level 1A MAGi normal rate", - logical_source="imap_mag_l1a_norm-magi", - logical_source_desc="IMAP Mission MAGi Normal Rate Instrument Level-1A Data.", - instrument_base=mag_base, -) - -mag_l1a_burst_mago_attrs = GlobalDataLevelAttrs( - "L1A_raw-burst-mago>Level 1A MAGo burst rate", - logical_source="imap_mag_l1a_burst-mago", - logical_source_desc="IMAP Mission MAGo Burst Rate Instrument Level-1A Data.", - instrument_base=mag_base, -) - -mag_l1a_burst_magi_attrs = GlobalDataLevelAttrs( - "L1A_raw-burst-magi>Level 1A MAGi burst rate", - logical_source="imap_mag_l1a_burst-magi", - logical_source_desc="IMAP Mission MAGi Burst Rate Instrument Level-1A Data.", - instrument_base=mag_base, -) - - -mag_l1b_attrs = GlobalDataLevelAttrs( - "L1B_SCI>Level-1B Science Data", - # TODO: replace "sci" with descriptor "norm" / "burst" - logical_source="imap_mag_l1b_sci", - logical_source_desc="IMAP Mission MAG Instrument Level-1B Data.", - instrument_base=mag_base, -) - -mag_l1c_attrs = GlobalDataLevelAttrs( - "L1C_SCI>Level-1C Science Data", - # TODO: replace "sci" with descriptor "norm" / "burst" - logical_source="imap_mag_l1c_sci", - logical_source_desc="IMAP Mission MAG Instrument Level-1C Data.", - instrument_base=mag_base, -) - -# TODO: Supporting data attributes? - -# TODO: display type, catdesc, units, format, label_axis - -# TODO: update descriptor to be more accurate for L1A raw -# TODO: does raw value need "counts" -mag_raw_vector_attrs = ScienceAttrs( - validmin=GlobalConstants.INT_FILLVAL, - validmax=GlobalConstants.INT_MAXVAL, - catdesc="Raw, unprocessed magnetic field vector data in bytes", - depend_0="epoch", - depend_1="direction", - display_type="time_series", - fieldname="Magnetic Field Vector", - label_axis="Raw binary magnetic field vector data", - fill_val=GlobalConstants.INT_MAXVAL, - format="I3", - units="counts", - var_type="data", -) - -vector_attrs = ScienceAttrs( - validmin=GlobalConstants.INT_FILLVAL, - validmax=GlobalConstants.INT_MAXVAL, - catdesc="Magnetic field vector with x, y, z, and sensor range, varying by time", - depend_0="epoch", - depend_1="direction", - display_type="time_series", - fieldname="Magnetic Field Vector", - label_axis="Magnetic field vector data", - fill_val=GlobalConstants.INT_MAXVAL, - format="I3", - units="counts", - var_type="data", -) - -mag_support_attrs = ScienceAttrs( - validmin=GlobalConstants.INT_FILLVAL, - validmax=GlobalConstants.INT_MAXVAL, - depend_0="epoch", - display_type="time_series", - fill_val=GlobalConstants.INT_FILLVAL, - format="I12", - var_type="support_data", -) - -mag_metadata_attrs = AttrBase( - validmin=GlobalConstants.INT_FILLVAL, - validmax=GlobalConstants.INT_MAXVAL, - display_type="time_series", - fill_val=GlobalConstants.INT_FILLVAL, - var_type="metadata", -) - - -mag_flag_attrs = ScienceAttrs( - validmin=0, - validmax=1, - depend_0="epoch", - display_type="time_series", - fill_val=255, - format="I1", -) - -raw_direction_attrs = AttrBase( - validmin=GlobalConstants.INT_FILLVAL, - validmax=GlobalConstants.INT_MAXVAL, - catdesc="Raw magnetic field vector binary length", - fieldname="Raw magnetic field vector binary length", - format="I3", - var_type="support_data", - display_type="time_series", - label_axis="Magnetic field vector directions", -) - -direction_attrs = AttrBase( - validmin=GlobalConstants.INT_FILLVAL, - validmax=GlobalConstants.INT_MAXVAL, - catdesc="Magnetic field vector", - fieldname="[x,y,z] magnetic field vector", - format="I3", - var_type="support_data", - display_type="time_series", - label_axis="Magnetic field vector", -) - -# Catdesc (<80 chars), Fieldnam (<30 chars) - -catdesc_fieldname_l0 = { - # TODO: Don't include PUS_SPARE1? - "PUS_SPARE1": ["Spare header from ESA Standard", "Spare header"], - "PUS_VERSION": ["PUS Version number", "PUS Version number"], - "PUS_SPARE2": ["Spare header from ESA Standard", "Spare header"], - "PUS_STYPE": ["PUS Service type", "PUS Service type"], - "PUS_SSUBTYPE": [ - "PUS Service subtype, the number of seconds of data minus 1", - "Number of seconds of data minus 1", - ], - "COMPRESSION": [ - "Indicates if the data is compressed, with 1 indicating the data is" - " compressed", - "Data is compressed", - ], - "MAGO_ACT": ["MAGO Active status", "MAGO Active status boolean"], - "MAGI_ACT": ["MAGI Active status", "MAGI Active status boolean"], - "PRI_SENS": [ - "Indicates which instrument is designated as primary. 0 is MAGo, 1 is MAGi", - "MAGi primary sensor boolean", - ], - "SPARE1": ["Spare", "Spare"], - "PRI_VECSEC": ["Primary vectors per second count", "Primary vectors per second"], - "SEC_VECSEC": [ - "Secondary vectors per second count", - "Secondary vectors per second", - ], - "SPARE2": ["Spare", "Spare"], - "PRI_COARSETM": [ - "Primary coarse time, mission SCLK in whole seconds", - "Primary coarse time (s)", - ], - "PRI_FNTM": [ - "Primary fine time, mission SCLK in 16bit subseconds", - "Primary fine time (16 bit subsecond)", - ], - "SEC_COARSETM": [ - "Secondary coarse time, mission SCLK in whole seconds", - "Secondary coarse time (s)", - ], - "SEC_FNTM": [ - "Secondary fine time, mission SCLK in 16bit subseconds", - "Secondary fine time (16 bit subsecond)", - ], - "VECTORS": [ - "Raw binary value of MAG Science vectors before processing", - "Raw vector binary", - ], -} - -catdesc_fieldname_l1a = { - "is_mago": [ - "Indicates if the data is from MAGo (True is MAGo, false is MAGi)", - "Data is from MAGo", - ], - "active": ["Indicates if the sensor is active", "Sensor is active"], - # TODO is this in CDF - "start_time": ["The coarse and fine time for the sensor", ""], - "vectors_per_second": [ - "Number of vectors measured per second, determined by the instrument mode", - "Vectors per second", - ], - "expected_vector_count": [ - "Expected number of vectors (vectors_per_second * seconds_of_data)", - "Expected vector count", - ], - "seconds_of_data": ["Number of seconds of data", "Seconds of data"], - "SHCOARSE": ["Mission elapsed time", "Mission elapsed time"], - "vectors": [ - "List of magnetic vector samples. Each sample is in the format [x,y,z,rng] " - "for the x, y, z coordinates of the field and the range of the instrument.", - "Magnetic field vectors, [x, y, z, range]", - ], -} diff --git a/imap_processing/tests/mag/test_mag_decom.py b/imap_processing/tests/mag/test_mag_decom.py index e2217275e..0e790b528 100644 --- a/imap_processing/tests/mag/test_mag_decom.py +++ b/imap_processing/tests/mag/test_mag_decom.py @@ -1,13 +1,23 @@ from pathlib import Path import pandas as pd +import pytest -from imap_processing.cdf import global_attrs +from imap_processing.cdf.imap_cdf_manager import ImapCdfAttributes from imap_processing.cdf.utils import load_cdf, write_cdf -from imap_processing.mag import mag_cdf_attrs +from imap_processing.mag.constants import DataMode from imap_processing.mag.l0.decom_mag import decom_packets, generate_dataset +@pytest.fixture() +def cdf_attrs(): + test_attrs = ImapCdfAttributes() + test_attrs.add_instrument_global_attrs("mag") + test_attrs.add_instrument_variable_attrs("mag", "l1") + test_attrs.add_global_attribute("Data_version", "v001") + return test_attrs + + def test_mag_decom(): current_directory = Path(__file__).parent burst_test_file = current_directory / "mag_l0_test_data.pkts" @@ -45,28 +55,34 @@ def test_mag_decom(): assert len(l0) == len(expected_output.index) -def test_mag_raw_xarray(): +def test_mag_raw_xarray(cdf_attrs): current_directory = Path(__file__).parent burst_test_file = current_directory / "mag_l0_test_data.pkts" packets = decom_packets(str(burst_test_file)) l0_norm = packets["norm"] l0_burst = packets["burst"] - norm_data = generate_dataset(l0_norm, mag_cdf_attrs.mag_l1a_norm_raw_attrs.output()) - burst_data = generate_dataset( - l0_burst, mag_cdf_attrs.mag_l1a_burst_raw_attrs.output() + norm_data = generate_dataset(l0_norm, DataMode.NORM, cdf_attrs) + burst_data = generate_dataset(l0_burst, DataMode.BURST, cdf_attrs) + + # Logical_file_id is filled in at file creation time. The rest of the required + # values should be included. + assert all( + [ + item is not None + for key, item in norm_data.attrs.items() + if key != "Logical_file_id" + ] ) - required_attrs = list( - global_attrs.GlobalInstrumentAttrs("", "", "").output().keys() + assert all( + [ + item is not None + for key, item in burst_data.attrs.items() + if key != "Logical_file_id" + ] ) - assert all([item in list(norm_data.attrs.keys()) for item in required_attrs]) - assert all([item is not None for _, item in norm_data.attrs.items()]) - - assert all([item in list(burst_data.attrs.keys()) for item in required_attrs]) - assert all([item is not None for _, item in burst_data.attrs.items()]) - expected_norm_len = 17 assert norm_data.sizes["epoch"] == expected_norm_len @@ -74,17 +90,15 @@ def test_mag_raw_xarray(): assert burst_data.sizes["epoch"] == expected_burst_len -def test_mag_raw_cdf_generation(): +def test_mag_raw_cdf_generation(cdf_attrs): current_directory = Path(__file__).parent test_file = current_directory / "mag_l0_test_data.pkts" packets = decom_packets(str(test_file)) l0_norm = packets["norm"] l0_burst = packets["burst"] - norm_data = generate_dataset(l0_norm, mag_cdf_attrs.mag_l1a_norm_raw_attrs.output()) - burst_data = generate_dataset( - l0_burst, mag_cdf_attrs.mag_l1a_burst_raw_attrs.output() - ) + norm_data = generate_dataset(l0_norm, DataMode.NORM, cdf_attrs) + burst_data = generate_dataset(l0_burst, DataMode.BURST, cdf_attrs) output = write_cdf(norm_data) assert output.exists() diff --git a/imap_processing/tests/mag/test_mag_l1a.py b/imap_processing/tests/mag/test_mag_l1a.py index 3584bfb02..ea7668f99 100644 --- a/imap_processing/tests/mag/test_mag_l1a.py +++ b/imap_processing/tests/mag/test_mag_l1a.py @@ -5,10 +5,11 @@ from imap_processing.cdf.utils import met_to_j2000ns from imap_processing.mag.l0.decom_mag import decom_packets -from imap_processing.mag.l1a.mag_l1a import process_packets +from imap_processing.mag.l1a.mag_l1a import mag_l1a, process_packets from imap_processing.mag.l1a.mag_l1a_data import ( MAX_FINE_TIME, MagL1a, + MagL1aPacketProperties, TimeTuple, ) @@ -22,7 +23,6 @@ def test_compare_validation_data(): l1 = process_packets(l0["norm"]) # Should have one day of data expected_day = np.datetime64("2023-11-30") - print(l1["mago"]) l1_mago = l1["mago"][expected_day] l1_magi = l1["magi"][expected_day] @@ -113,3 +113,70 @@ def test_calculate_vector_time(): ] ) assert (test_data == expected_data).all() + + +def test_mag_l1a_data(): + test_vectors = np.array( + [ + [1, 2, 3, 4, 10000], + [5, 6, 7, 8, 10050], + [9, 10, 11, 12, 10100], + [13, 13, 11, 12, 10150], + ], + dtype=np.uint, + ) + test_vecsec = 2 + start_time = TimeTuple(10000, 0) + + packet_properties = MagL1aPacketProperties( + 435954628, start_time, test_vecsec, 1, 0, 0, 1 + ) + mag_l1a = MagL1a(True, True, 10000, test_vectors, packet_properties) + + new_vectors = np.array( + [[13, 14, 15, 16, 10400], [16, 17, 18, 19, 10450]], dtype=np.uint + ) + + new_seq = 5 + new_properties = MagL1aPacketProperties( + 435954628, TimeTuple(10400, 0), test_vecsec, 0, new_seq, 0, 1 + ) + mag_l1a.append_vectors(new_vectors, new_properties) + + assert np.array_equal( + mag_l1a.vectors, + np.array( + [ + [1, 2, 3, 4, 10000], + [5, 6, 7, 8, 10050], + [9, 10, 11, 12, 10100], + [13, 13, 11, 12, 10150], + [13, 14, 15, 16, 10400], + [16, 17, 18, 19, 10450], + ], + dtype=np.uint, + ), + ) + assert mag_l1a.missing_sequences == [1, 2, 3, 4] + + +def test_mag_l1a(): + current_directory = Path(__file__).parent + test_file = current_directory / "mag_l1_test_data.pkts" + + output_data = mag_l1a(test_file, "v001") + + # Test data is one day's worth of NORM data, so it should return one raw, one MAGO + # and one MAGI dataset + assert len(output_data) == 3 + expected_logical_source = [ + "imap_mag_l1a_norm-raw", + "imap_mag_l1a_norm-mago", + "imap_mag_l1a_norm-magi", + ] + + for data_type in [data.attrs["Logical_source"] for data in output_data]: + assert data_type in expected_logical_source + + for data in output_data: + assert data.attrs["Data_version"] == "v001"