diff --git a/imap_processing/mag/l0/decom_mag.py b/imap_processing/mag/l0/decom_mag.py index 49aade53e..342c3cded 100644 --- a/imap_processing/mag/l0/decom_mag.py +++ b/imap_processing/mag/l0/decom_mag.py @@ -21,7 +21,7 @@ logger = logging.getLogger(__name__) -def decom_packets(packet_file_path: str | Path) -> list[MagL0]: +def decom_packets(packet_file_path: str | Path) -> dict[str, list[MagL0]]: """Decom MAG data packets using MAG packet definition. Parameters @@ -31,9 +31,9 @@ def decom_packets(packet_file_path: str | Path) -> list[MagL0]: Returns ------- - data : list[MagL0] - A list of MAG L0 data classes, including both burst and normal packets. (the - packet type is defined in each instance of L0.) + data_dict : dict[str, list[MagL0]] + A dict with 2 keys pointing to lists of MAG L0 data classes. "norm" corresponds + to normal mode packets, "burst" corresponds to burst mode packets. """ # Define paths xtce_document = Path( @@ -43,7 +43,8 @@ def decom_packets(packet_file_path: str | Path) -> list[MagL0]: packet_definition = xtcedef.XtcePacketDefinition(xtce_document) mag_parser = parser.PacketParser(packet_definition) - data_list = [] + norm_data = [] + burst_data = [] with open(packet_file_path, "rb") as binary_data: mag_packets = mag_parser.generator(binary_data) @@ -57,59 +58,17 @@ def decom_packets(packet_file_path: str | Path) -> list[MagL0]: else item.raw_value for item in packet.data.values() ] - data_list.append(MagL0(CcsdsData(packet.header), *values)) - - return data_list - - -def export_to_xarray(l0_data: list[MagL0]): - """Generate xarray files for "raw" MAG CDF files from MagL0 data. - - Mag outputs "RAW" CDF files just after decomming. These have the immediately - post-decom data, with raw binary data for the vectors instead of vector values. - - Parameters - ---------- - l0_data: list[MagL0] - A list of MagL0 datapoints - - Returns - ------- - norm_data : xr.Dataset - xarray dataset for generating burst data CDFs - burst_data : xr.Dataset - xarray dataset for generating burst data CDFs - """ - # TODO split by mago and magi using primary sensor - norm_data = [] - burst_data = [] + if apid == Mode.NORMAL: + norm_data.append(MagL0(CcsdsData(packet.header), *values)) + else: + burst_data.append(MagL0(CcsdsData(packet.header), *values)) - for packet in l0_data: - if packet.ccsds_header.PKT_APID == Mode.NORMAL: - norm_data.append(packet) - if packet.ccsds_header.PKT_APID == Mode.BURST: - burst_data.append(packet) - - norm_dataset = None - burst_dataset = None - - if len(norm_data) > 0: - norm_dataset = generate_dataset( - norm_data, mag_cdf_attrs.mag_l1a_norm_raw_attrs.output() - ) - if len(burst_data) > 0: - burst_dataset = generate_dataset( - burst_data, mag_cdf_attrs.mag_l1a_burst_raw_attrs.output() - ) - - return norm_dataset, burst_dataset + return {"norm": norm_data, "burst": burst_data} def generate_dataset(l0_data: list[MagL0], dataset_attrs: dict) -> xr.Dataset: """ - Generate a CDF dataset from the sorted L0 MAG data. - - Used to create 2 similar datasets, for norm and burst data. + Generate a CDF dataset from the sorted raw L0 MAG data. Parameters ---------- @@ -124,6 +83,8 @@ def generate_dataset(l0_data: list[MagL0], dataset_attrs: dict) -> xr.Dataset: dataset : xr.Dataset xarray dataset with proper CDF attributes and shape. """ + # TODO: Correct CDF attributes from email + vector_data = np.zeros((len(l0_data), len(l0_data[0].VECTORS))) shcoarse_data = np.zeros(len(l0_data), dtype="datetime64[ns]") diff --git a/imap_processing/mag/l0/mag_l0_data.py b/imap_processing/mag/l0/mag_l0_data.py index edab08c77..f9c41b8b7 100644 --- a/imap_processing/mag/l0/mag_l0_data.py +++ b/imap_processing/mag/l0/mag_l0_data.py @@ -1,8 +1,12 @@ """Dataclasses for Level 0 MAG data.""" +from __future__ import annotations + from dataclasses import dataclass from enum import IntEnum +import numpy as np + from imap_processing.ccsds.ccsds_data import CcsdsData @@ -66,9 +70,10 @@ class MagL0: Secondary Coarse Time for first vector, seconds SEC_FNTM: int Secondary Fine Time for first vector, subseconds - VECTORS: bin + VECTORS: np.ndarray | str MAG Science Vectors - divide based on PRI_VECSEC and PUS_SSUBTYPE for vector - counts + counts. There is a post init call to convert a string into a numpy array - + the only place it is a string is in the class initialization. """ ccsds_header: CcsdsData @@ -90,7 +95,7 @@ class MagL0: PRI_FNTM: int SEC_COARSETM: int SEC_FNTM: int - VECTORS: bytearray + VECTORS: np.ndarray | str def __post_init__(self): """Convert Vectors attribute from string to bytearray if needed. @@ -98,11 +103,17 @@ def __post_init__(self): Also convert encoded "VECSEC" (vectors per second) into proper vectors per second values """ - if isinstance(self.VECTORS, str): - # Convert string output from space_packet_parser to bytearray - self.VECTORS = bytearray( - int(self.VECTORS, 2).to_bytes(len(self.VECTORS) // 8, "big") - ) + # Convert string output from space_packet_parser to numpy array of + # big-endian bytes + self.VECTORS = np.frombuffer( + int(self.VECTORS, 2).to_bytes(len(self.VECTORS) // 8, "big"), + dtype=np.dtype(">b"), + ) + + # Remove buffer from end of vectors. Vector data needs to be in 50 bit chunks, + # and may have an extra byte at the end from CCSDS padding. + if len(self.VECTORS) % 2: + self.VECTORS = self.VECTORS[:-1] self.PRI_VECSEC = 2**self.PRI_VECSEC self.SEC_VECSEC = 2**self.SEC_VECSEC diff --git a/imap_processing/mag/l1a/mag_l1a.py b/imap_processing/mag/l1a/mag_l1a.py index fe8f4b915..c45b4477c 100644 --- a/imap_processing/mag/l1a/mag_l1a.py +++ b/imap_processing/mag/l1a/mag_l1a.py @@ -3,7 +3,10 @@ import logging from imap_processing.cdf.utils import write_cdf +from imap_processing.mag import mag_cdf_attrs 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 logger = logging.getLogger(__name__) @@ -17,14 +20,115 @@ def mag_l1a(packet_filepath): packet_filepath : Packet files for processing """ - mag_l0 = decom_mag.decom_packets(packet_filepath) + packets = decom_mag.decom_packets(packet_filepath) - mag_norm, mag_burst = decom_mag.export_to_xarray(mag_l0) + norm_data = packets["norm"] + burst_data = packets["burst"] - if mag_norm is not None: - file = write_cdf(mag_norm) - logger.info(f"Created CDF file at {file}") + if norm_data is not None: + mag_norm_raw = decom_mag.generate_dataset( + norm_data, mag_cdf_attrs.mag_l1a_norm_raw_attrs.output() + ) + file = write_cdf(mag_norm_raw) + logger.info(f"Created RAW CDF file at {file}") - if mag_burst is not None: - file = write_cdf(mag_burst) - logger.info(f"Created CDF file at {file}") + if burst_data is not None: + mag_burst_raw = decom_mag.generate_dataset( + burst_data, mag_cdf_attrs.mag_l1a_burst_raw_attrs.output() + ) + file = write_cdf(mag_burst_raw) + logger.info(f"Created RAW CDF file at {file}") + + +def process_packets(mag_l0_list: list[MagL0]) -> dict[str, list[MagL1a]]: + """ + Given a list of MagL0 packets, process them into MagO and MagI L1A data classes. + + Parameters + ---------- + mag_l0_list : list[MagL0] + List of Mag L0 packets to process + + Returns + ------- + packet_dict: dict[str, list[MagL1a]] + Dictionary containing two keys: "mago" which points to a list of mago MagL1A + objects, and "magi" which points to a list of magi MagL1A objects. + + """ + magi = [] + mago = [] + + for mag_l0 in mag_l0_list: + if mag_l0.COMPRESSION: + raise NotImplementedError("Unable to process compressed data") + + primary_start_time = TimeTuple(mag_l0.PRI_COARSETM, mag_l0.PRI_FNTM) + secondary_start_time = TimeTuple(mag_l0.SEC_COARSETM, mag_l0.SEC_FNTM) + + # seconds of data in this packet is the SUBTYPE plus 1 + seconds_per_packet = mag_l0.PUS_SSUBTYPE + 1 + + # now we know the number of seconds 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 + + primary_vectors, secondary_vectors = MagL1a.process_vector_data( + mag_l0.VECTORS, total_primary_vectors, total_secondary_vectors + ) + + # Primary sensor is MAGO (most common expected case) + if mag_l0.PRI_SENS == 0: + mago_l1a = MagL1a( + True, + bool(mag_l0.MAGO_ACT), + primary_start_time, + mag_l0.PRI_VECSEC, + total_primary_vectors, + seconds_per_packet, + mag_l0.SHCOARSE, + primary_vectors, + ) + + magi_l1a = MagL1a( + False, + bool(mag_l0.MAGI_ACT), + secondary_start_time, + mag_l0.SEC_VECSEC, + total_secondary_vectors, + seconds_per_packet, + mag_l0.SHCOARSE, + secondary_vectors, + ) + # Primary sensor is MAGI + else: + magi_l1a = MagL1a( + False, + bool(mag_l0.MAGI_ACT), + primary_start_time, + mag_l0.PRI_VECSEC, + total_primary_vectors, + seconds_per_packet, + mag_l0.SHCOARSE, + primary_vectors, + ) + + mago_l1a = MagL1a( + True, + bool(mag_l0.MAGO_ACT), + secondary_start_time, + mag_l0.SEC_VECSEC, + total_secondary_vectors, + seconds_per_packet, + mag_l0.SHCOARSE, + secondary_vectors, + ) + + mago.append(mago_l1a) + magi.append(magi_l1a) + + return {"mago": mago, "magi": magi} diff --git a/imap_processing/mag/l1a/mag_l1a_data.py b/imap_processing/mag/l1a/mag_l1a_data.py new file mode 100644 index 000000000..4bdbb0ef8 --- /dev/null +++ b/imap_processing/mag/l1a/mag_l1a_data.py @@ -0,0 +1,254 @@ +"""Data classes for storing and processing MAG Level 1A data.""" + +from dataclasses import dataclass +from math import floor + +import numpy as np + +MAX_FINE_TIME = 65535 # maximum 16 bit unsigned int + + +@dataclass +class TimeTuple: + """ + Class for storing fine time/coarse time for MAG data. + + Course time is mission SCLK in seconds. Fine time is 16bit unsigned sub-second + counter. + """ + + coarse_time: int + fine_time: int + + def __add__(self, seconds: float): + """ + Add a number of seconds to the time tuple. + + Parameters + ---------- + seconds : float + Number of seconds to add + + Returns + ------- + time : TimeTuple + New time tuple with the current time tuple + seconds. + """ + # Add whole seconds to coarse time + coarse = self.coarse_time + floor(seconds) + # fine time is 1/65535th of a second + fine = self.fine_time + round((seconds % 1) * MAX_FINE_TIME) + + # If fine is larger than the max, move the excess into coarse time. + if fine > MAX_FINE_TIME: + coarse = coarse + floor(fine / MAX_FINE_TIME) + fine = fine % MAX_FINE_TIME + + return TimeTuple(coarse, fine) + + +@dataclass +class Vector: + """ + Data class for storing MAG vector data. + + Attributes + ---------- + timestamp : TimeTuple + Time of the vector sample + vectors : tuple[int, int, int, int] + Vector sample, containing x,y,z,range + """ + + # TODO: This timestamp should be in J2000/datetime64 + timestamp: TimeTuple + x: int + y: int + z: int + rng: int + + def __init__(self, vectors, time: TimeTuple): + self.timestamp = time + self.x = vectors[0] + self.y = vectors[1] + self.z = vectors[2] + self.rng = vectors[3] + + +@dataclass +class MagL1a: + """ + Data class for MAG Level 1A data. + + 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 + start_time : TimeTuple + The coarse and fine time for the sensor + vectors_per_second : int + Number of vectors per second + expected_vector_count : int + Expected number of vectors (vectors_per_second * seconds_of_data) + seconds_of_data : int + Number of seconds of data + SHCOARSE : int + Mission elapsed time + vectors : list[Vector] + List of magnetic vector samples, starting at start_time + + """ + + is_mago: bool + active: bool + start_time: TimeTuple + vectors_per_second: int + expected_vector_count: int + seconds_of_data: int + SHCOARSE: int + vectors: list + + def __post_init__(self): + """ + Convert the vector list to a vector list with timestamps associated. + + The first vector starts at start_time, then each subsequent vector time is + computed by adding 1/vectors_per_second to the previous vector's time. + + This replaces self.vectors with a list of Vector objects. + """ + sample_time_interval = 1 / self.vectors_per_second + current_time = self.start_time + for index, vector in enumerate(self.vectors): + self.vectors[index] = Vector(vector, current_time) + current_time = self.vectors[index].timestamp + sample_time_interval + + @staticmethod + def process_vector_data( + vector_data: np.ndarray, primary_count: int, secondary_count: int + ) -> (list[tuple], list[tuple]): + """ + Given raw packet data, process into Vectors. + + Vectors are grouped into primary sensor and secondary sensor, and returned as a + tuple (primary sensor vectors, secondary sensor vectors) + + Written by MAG instrument team + + Parameters + ---------- + vector_data : np.ndarray + Raw vector data, in bytes. Contains both primary and secondary vector data + (first primary, then secondary) + primary_count : int + Count of the number of primary vectors + secondary_count : int + Count of the number of secondary vectors + + Returns + ------- + (primary, secondary) + Two arrays, each containing tuples of (x, y, z, sample_range) for each + vector sample. + """ + + # TODO: error handling + def to_signed16(n): + n = n & 0xFFFF + return n | (-(n & 0x8000)) + + pos = 0 + primary_vectors = [] + secondary_vectors = [] + + # Since the vectors are stored as 50 bit chunks but accessed via hex (4 bit + # chunks) there is some shifting required for processing the bytes. + # However, from a bit processing perspective, the first 48 bits of each 50 bit + # chunk corresponds to 3 16 bit signed integers. The last 2 bits are the sensor + # range. + + for i in range(primary_count + secondary_count): # 0..63 say + x, y, z, rng = 0, 0, 0, 0 + if i % 4 == 0: # start at bit 0, take 8 bits + 8bits + # pos = 0, 25, 50... + x = ( + ((vector_data[pos + 0] & 0xFF) << 8) + | ((vector_data[pos + 1] & 0xFF) << 0) + ) & 0xFFFF + y = ( + ((vector_data[pos + 2] & 0xFF) << 8) + | ((vector_data[pos + 3] & 0xFF) << 0) + ) & 0xFFFF + z = ( + ((vector_data[pos + 4] & 0xFF) << 8) + | ((vector_data[pos + 5] & 0xFF) << 0) + ) & 0xFFFF + rng = (vector_data[pos + 6] >> 6) & 0x3 + pos += 6 + elif i % 4 == 1: # start at bit 2, take 6 bits, 8 bit, 2 bits per vector + # pos = 6, 31... + x = ( + ((vector_data[pos + 0] & 0x3F) << 10) + | ((vector_data[pos + 1] & 0xFF) << 2) + | ((vector_data[pos + 2] >> 6) & 0x03) + ) & 0xFFFF + y = ( + ((vector_data[pos + 2] & 0x3F) << 10) + | ((vector_data[pos + 3] & 0xFF) << 2) + | ((vector_data[pos + 4] >> 6) & 0x03) + ) & 0xFFFF + z = ( + ((vector_data[pos + 4] & 0x3F) << 10) + | ((vector_data[pos + 5] & 0xFF) << 2) + | ((vector_data[pos + 6] >> 6) & 0x03) + ) & 0xFFFF + rng = (vector_data[pos + 6] >> 4) & 0x3 + pos += 6 + elif i % 4 == 2: # start at bit 4, take 4 bits, 8 bits, 4 bits per vector + # pos = 12, 37... + x = ( + ((vector_data[pos + 0] & 0x0F) << 12) + | ((vector_data[pos + 1] & 0xFF) << 4) + | ((vector_data[pos + 2] >> 4) & 0x0F) + ) & 0xFFFF + y = ( + ((vector_data[pos + 2] & 0x0F) << 12) + | ((vector_data[pos + 3] & 0xFF) << 4) + | ((vector_data[pos + 4] >> 4) & 0x0F) + ) & 0xFFFF + z = ( + ((vector_data[pos + 4] & 0x0F) << 12) + | ((vector_data[pos + 5] & 0xFF) << 4) + | ((vector_data[pos + 6] >> 4) & 0x0F) + ) & 0xFFFF + rng = (vector_data[pos + 6] >> 2) & 0x3 + pos += 6 + elif i % 4 == 3: # start at bit 6, take 2 bits, 8 bits, 6 bits per vector + # pos = 18, 43... + x = ( + ((vector_data[pos + 0] & 0x03) << 14) + | ((vector_data[pos + 1] & 0xFF) << 6) + | ((vector_data[pos + 2] >> 2) & 0x3F) + ) & 0xFFFF + y = ( + ((vector_data[pos + 2] & 0x03) << 14) + | ((vector_data[pos + 3] & 0xFF) << 6) + | ((vector_data[pos + 4] >> 2) & 0x3F) + ) & 0xFFFF + z = ( + ((vector_data[pos + 4] & 0x03) << 14) + | ((vector_data[pos + 5] & 0xFF) << 6) + | ((vector_data[pos + 6] >> 2) & 0x3F) + ) & 0xFFFF + rng = (vector_data[pos + 6] >> 0) & 0x3 + pos += 7 + + vector = (to_signed16(x), to_signed16(y), to_signed16(z), rng) + if i < primary_count: + primary_vectors.append(vector) + else: + secondary_vectors.append(vector) + + return (primary_vectors, secondary_vectors) diff --git a/imap_processing/tests/mag/mag_l1a_test_output.csv b/imap_processing/tests/mag/mag_l1a_test_output.csv new file mode 100644 index 000000000..9cecd8273 --- /dev/null +++ b/imap_processing/tests/mag/mag_l1a_test_output.csv @@ -0,0 +1,97 @@ +sequence,x_pri,y_pri,z_pri,rng_pri,x_sec,y_sec,z_sec,rng_sec,pri_coarse,pri_fine,sec_coarse,sec_fine +5,13324,24467,20105,3,13324,24467,20104,3,439067318,64618,439067318,64612 +5,26623,-16643,-25350,3,26623,-16644,-25351,3,439067318,64618,439067318,64612 +5,26626,-16630,-25324,3,26626,-16631,-25325,3,439067318,64618,439067318,64612 +5,26629,-16617,-25298,3,26629,-16617,-25298,3,439067318,64618,439067318,64612 +5,26633,-16604,-25271,3,26633,-16604,-25272,3,439067318,64618,439067318,64612 +5,26636,-16591,-25245,3,26636,-16591,-25246,3,439067318,64618,439067318,64612 +5,26639,-16577,-25218,3,26639,-16578,-25219,3,439067318,64618,439067318,64612 +5,26643,-16564,-25192,3,26642,-16565,-25193,3,439067318,64618,439067318,64612 +5,26646,-16551,-25166,3,26646,-16552,-25167,3,439067318,64618,439067318,64612 +5,26649,-16538,-25139,3,26649,-16538,-25140,3,439067318,64618,439067318,64612 +5,26652,-16525,-25113,3,26652,-16525,-25114,3,439067318,64618,439067318,64612 +5,26656,-16512,-25087,3,26656,-16512,-25087,3,439067318,64618,439067318,64612 +5,26659,-16498,-25060,3,26659,-16499,-25061,3,439067318,64618,439067318,64612 +5,26662,-16485,-25034,3,26662,-16486,-25035,3,439067318,64618,439067318,64612 +5,26666,-16472,-25007,3,26666,-16472,-25008,3,439067318,64618,439067318,64612 +5,26669,-16459,-24981,3,26669,-16459,-24982,3,439067318,64618,439067318,64612 +6,26672,-16446,-24955,3,26672,-16446,-24956,3,439067307,64605,439067307,64600 +6,26676,-16432,-24928,3,26675,-16433,-24929,3,439067307,64605,439067307,64600 +6,26679,-16419,-24902,3,26679,-16420,-24903,3,439067307,64605,439067307,64600 +6,26682,-16406,-24876,3,26682,-16407,-24877,3,439067307,64605,439067307,64600 +6,26685,-16393,-24849,3,26685,-16393,-24850,3,439067307,64605,439067307,64600 +6,26689,-16380,-24823,3,26689,-16380,-24824,3,439067307,64605,439067307,64600 +6,26692,-16367,-24797,3,26692,-16367,-24797,3,439067307,64605,439067307,64600 +6,26695,-16353,-24770,3,26695,-16354,-24771,3,439067307,64605,439067307,64600 +6,26699,-16340,-24744,3,26698,-16341,-24745,3,439067307,64605,439067307,64600 +6,26702,-16327,-24717,3,26702,-16327,-24718,3,439067307,64605,439067307,64600 +6,26705,-16314,-24691,3,26705,-16314,-24692,3,439067307,64605,439067307,64600 +6,26708,-16301,-24665,3,26708,-16301,-24666,3,439067307,64605,439067307,64600 +6,26712,-16287,-24638,3,26712,-16288,-24639,3,439067307,64605,439067307,64600 +6,26715,-16274,-24612,3,26715,-16275,-24613,3,439067307,64605,439067307,64600 +6,26718,-16261,-24586,3,26718,-16261,-24586,3,439067307,64605,439067307,64600 +6,26722,-16248,-24559,3,26722,-16248,-24560,3,439067307,64605,439067307,64600 +7,26725,-16235,-24533,3,26725,-16235,-24534,3,439067315,64605,439067315,64600 +7,26728,-16222,-24507,3,26728,-16222,-24507,3,439067315,64605,439067315,64600 +7,26732,-16208,-24480,3,26731,-16209,-24481,3,439067315,64605,439067315,64600 +7,26735,-16195,-24454,3,26735,-16196,-24455,3,439067315,64605,439067315,64600 +7,26738,-16182,-24427,3,26738,-16182,-24428,3,439067315,64605,439067315,64600 +7,26741,-16169,-24401,3,26741,-16169,-24402,3,439067315,64605,439067315,64600 +7,26745,-16156,-24375,3,26745,-16156,-24376,3,439067315,64605,439067315,64600 +7,26748,-16142,-24348,3,26748,-16143,-24349,3,439067315,64605,439067315,64600 +7,26751,-16129,-24322,3,26751,-16130,-24323,3,439067315,64605,439067315,64600 +7,26755,-16116,-24296,3,26755,-16116,-24296,3,439067315,64605,439067315,64600 +7,26758,-16103,-24269,3,26758,-16103,-24270,3,439067315,64605,439067315,64600 +7,26761,-16090,-24243,3,26761,-16090,-24244,3,439067315,64605,439067315,64600 +7,26765,-16076,-24216,3,26764,-16077,-24217,3,439067315,64605,439067315,64600 +7,26768,-16063,-24190,3,26768,-16064,-24191,3,439067315,64605,439067315,64600 +7,26771,-16050,-24164,3,26771,-16051,-24165,3,439067315,64605,439067315,64600 +7,26774,-16037,-24137,3,26774,-16037,-24138,3,439067315,64605,439067315,64600 +8,26778,-16024,-24111,3,26778,-16024,-24112,3,439067323,64606,439067323,64601 +8,26781,-16011,-24085,3,26781,-16011,-24086,3,439067323,64606,439067323,64601 +8,26784,-15997,-24058,3,26784,-15998,-24059,3,439067323,64606,439067323,64601 +8,26788,-15984,-24032,3,26787,-15985,-24033,3,439067323,64606,439067323,64601 +8,26791,-15971,-24006,3,26791,-15971,-24006,3,439067323,64606,439067323,64601 +8,26794,-15958,-23979,3,26794,-15958,-23980,3,439067323,64606,439067323,64601 +8,26797,-15945,-23953,3,26797,-15945,-23954,3,439067323,64606,439067323,64601 +8,26801,-15931,-23926,3,26801,-15932,-23927,3,439067323,64606,439067323,64601 +8,26804,-15918,-23900,3,26804,-15919,-23901,3,439067323,64606,439067323,64601 +8,26807,-15905,-23874,3,26807,-15906,-23875,3,439067323,64606,439067323,64601 +8,26811,-15892,-23847,3,26811,-15892,-23848,3,439067323,64606,439067323,64601 +8,26814,-15879,-23821,3,26814,-15879,-23822,3,439067323,64606,439067323,64601 +8,26817,-15866,-23795,3,26817,-15866,-23795,3,439067323,64606,439067323,64601 +8,26821,-15852,-23768,3,26820,-15853,-23769,3,439067323,64606,439067323,64601 +8,26824,-15839,-23742,3,26824,-15840,-23743,3,439067323,64606,439067323,64601 +8,26827,-15826,-23715,3,26827,-15826,-23716,3,439067323,64606,439067323,64601 +9,26830,-15813,-23689,3,26830,-15813,-23690,3,439067331,64604,439067331,64601 +9,26834,-15800,-23663,3,26834,-15800,-23664,3,439067331,64604,439067331,64601 +9,26837,-15786,-23636,3,26837,-15787,-23637,3,439067331,64604,439067331,64601 +9,26840,-15773,-23610,3,26840,-15774,-23611,3,439067331,64604,439067331,64601 +9,26844,-15760,-23584,3,26843,-15761,-23585,3,439067331,64604,439067331,64601 +9,26847,-15747,-23557,3,26847,-15747,-23558,3,439067331,64604,439067331,64601 +9,26850,-15734,-23531,3,26850,-15734,-23532,3,439067331,64604,439067331,64601 +9,26853,-15721,-23505,3,26853,-15721,-23505,3,439067331,64604,439067331,64601 +9,26857,-15707,-23478,3,26857,-15708,-23479,3,439067331,64604,439067331,64601 +9,26860,-15694,-23452,3,26860,-15695,-23453,3,439067331,64604,439067331,64601 +9,26863,-15681,-23425,3,26863,-15681,-23426,3,439067331,64604,439067331,64601 +9,26867,-15668,-23399,3,26867,-15668,-23400,3,439067331,64604,439067331,64601 +9,26870,-15655,-23373,3,26870,-15655,-23374,3,439067331,64604,439067331,64601 +9,26873,-15641,-23346,3,26873,-15642,-23347,3,439067331,64604,439067331,64601 +9,26877,-15628,-23320,3,26876,-15629,-23321,3,439067331,64604,439067331,64601 +9,26880,-15615,-23294,3,26880,-15616,-23294,3,439067331,64604,439067331,64601 +10,26883,-15602,-23267,3,26883,-15602,-23268,3,439067339,64603,439067339,64601 +10,26886,-15589,-23241,3,26886,-15589,-23242,3,439067339,64603,439067339,64601 +10,26890,-15576,-23215,3,26890,-15576,-23215,3,439067339,64603,439067339,64601 +10,26893,-15562,-23188,3,26893,-15563,-23189,3,439067339,64603,439067339,64601 +10,26896,-15549,-23162,3,26896,-15550,-23163,3,439067339,64603,439067339,64601 +10,26900,-15536,-23135,3,26900,-15536,-23136,3,439067339,64603,439067339,64601 +10,26903,-15523,-23109,3,26903,-15523,-23110,3,439067339,64603,439067339,64601 +10,26906,-15510,-23083,3,26906,-15510,-23084,3,439067339,64603,439067339,64601 +10,26910,-15496,-23056,3,26909,-15497,-23057,3,439067339,64603,439067339,64601 +10,26913,-15483,-23030,3,26913,-15484,-23031,3,439067339,64603,439067339,64601 +10,26916,-15470,-23004,3,26916,-15470,-23004,3,439067339,64603,439067339,64601 +10,26919,-15457,-22977,3,26919,-15457,-22978,3,439067339,64603,439067339,64601 +10,26923,-15444,-22951,3,26923,-15444,-22952,3,439067339,64603,439067339,64601 +10,26926,-15430,-22924,3,26926,-15431,-22925,3,439067339,64603,439067339,64601 +10,26929,-15417,-22898,3,26929,-15418,-22899,3,439067339,64603,439067339,64601 +10,26933,-15404,-22872,3,26932,-15405,-22873,3,439067339,64603,439067339,64601 diff --git a/imap_processing/tests/mag/test_mag_decom.py b/imap_processing/tests/mag/test_mag_decom.py index d085b7451..f02c13e6a 100644 --- a/imap_processing/tests/mag/test_mag_decom.py +++ b/imap_processing/tests/mag/test_mag_decom.py @@ -5,13 +5,16 @@ from imap_processing.cdf import global_attrs from imap_processing.cdf.utils import write_cdf -from imap_processing.mag.l0.decom_mag import decom_packets, export_to_xarray +from imap_processing.mag import mag_cdf_attrs +from imap_processing.mag.l0.decom_mag import decom_packets, generate_dataset def test_mag_decom(): current_directory = Path(__file__).parent burst_test_file = current_directory / "mag_l0_test_data.pkts" - l0 = decom_packets(burst_test_file) + packets = decom_packets(str(burst_test_file)) + + l0 = packets["burst"] + packets["norm"] expected_output = pd.read_csv(current_directory / "mag_l0_test_output.csv") for index, test in enumerate(l0): @@ -33,15 +36,28 @@ def test_mag_decom(): assert test.SEC_COARSETM == expected_output["SEC_COARSETM"][index] assert test.SEC_FNTM == expected_output["SEC_FNTM"][index] + # Remove bytes for header and previous attributes from CCSDS_HEX, + # remaining bytes are vectors + assert ( + test.VECTORS.tobytes().hex() + == expected_output["CCSDS_HEX"][index][54:].lower() + ) + assert len(l0) == len(expected_output.index) def test_mag_raw_xarray(): current_directory = Path(__file__).parent burst_test_file = current_directory / "mag_l0_test_data.pkts" - l0 = decom_packets(str(burst_test_file)) + 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, burst_data = export_to_xarray(l0) required_attrs = list( global_attrs.GlobalInstrumentAttrs("", "", "").output().keys() ) @@ -62,9 +78,14 @@ def test_mag_raw_xarray(): def test_mag_raw_cdf_generation(): current_directory = Path(__file__).parent test_file = current_directory / "mag_l0_test_data.pkts" - l0 = decom_packets(str(test_file)) + packets = decom_packets(str(test_file)) + l0_norm = packets["norm"] + l0_burst = packets["burst"] - norm_data, burst_data = export_to_xarray(l0) + 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() + ) 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 new file mode 100644 index 000000000..42826c26a --- /dev/null +++ b/imap_processing/tests/mag/test_mag_l1a.py @@ -0,0 +1,115 @@ +from pathlib import Path + +import numpy as np +import pandas as pd + +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_data import ( + MAX_FINE_TIME, + MagL1a, + TimeTuple, + Vector, +) + + +def test_compare_validation_data(): + current_directory = Path(__file__).parent + test_file = current_directory / "mag_l1_test_data.pkts" + # Test file contains only normal packets + l0 = decom_packets(test_file) + + l1 = process_packets(l0["norm"]) + l1_mago = l1["mago"] + l1_magi = l1["magi"] + + assert len(l1_mago) == 6 + assert len(l1_magi) == 6 + + validation_data = pd.read_csv(current_directory / "mag_l1a_test_output.csv") + + vector_index = 0 + + # Validation data does not have differing timestamps + for index in validation_data.index: + # Sequence in validation data starts at 5 + # Mago is primary, Magi is secondary in test data + l1_pri = l1_mago[validation_data["sequence"][index] - 5] + l1_sec = l1_magi[validation_data["sequence"][index] - 5] + + assert l1_pri.vectors[vector_index].x == validation_data["x_pri"][index] + assert l1_pri.vectors[vector_index].y == validation_data["y_pri"][index] + assert l1_pri.vectors[vector_index].z == validation_data["z_pri"][index] + assert l1_pri.vectors[vector_index].rng == validation_data["rng_pri"][index] + + assert l1_sec.vectors[vector_index].x == validation_data["x_sec"][index] + assert l1_sec.vectors[vector_index].y == validation_data["y_sec"][index] + assert l1_sec.vectors[vector_index].z == validation_data["z_sec"][index] + assert l1_sec.vectors[vector_index].rng == validation_data["rng_sec"][index] + + vector_index = ( + 0 if vector_index == l1_pri.expected_vector_count - 1 else vector_index + 1 + ) + + +def test_process_vector_data(): + expected_vector_data = [[1001, 1002, -3001, 3], [2001, -2002, -3333, 1]] + + # 100 bits, created by hand by appending all bits from expected_vector_data into one + # hex string, with range being 2 bits (so the second half is offset from the hex + # values) + hex_string = "03E903EAF447C1F47E0BBCBED0" + input_data = np.frombuffer(bytes.fromhex(hex_string), dtype=np.dtype(">b")) + + total_primary_vectors = 1 + total_secondary_vectors = 1 + + # 36 bytes + (primary_vectors, secondary_vectors) = MagL1a.process_vector_data( + input_data, total_primary_vectors, total_secondary_vectors + ) + + assert primary_vectors[0][0] == expected_vector_data[0][0] + assert primary_vectors[0][1] == expected_vector_data[0][1] + assert primary_vectors[0][2] == expected_vector_data[0][2] + assert primary_vectors[0][3] == expected_vector_data[0][3] + + assert secondary_vectors[0][0] == expected_vector_data[1][0] + assert secondary_vectors[0][1] == expected_vector_data[1][1] + assert secondary_vectors[0][2] == expected_vector_data[1][2] + assert secondary_vectors[0][3] == expected_vector_data[1][3] + + +def test_time_tuple(): + test_time_tuple = TimeTuple(439067318, 64618) + + test_add = test_time_tuple + 2 + + assert test_add == TimeTuple(439067320, 64618) + + # 1 / MAX_FINE_TIME + test_add = test_time_tuple + 1 / MAX_FINE_TIME + + assert test_add == TimeTuple(439067318, 64619) + + test_add = test_time_tuple + (1000 / MAX_FINE_TIME) + + assert test_add == TimeTuple(439067319, 83) + + +def test_vector_time(): + test_data = MagL1a( + True, + True, + TimeTuple(10000, 0), + 2, # 2 vectors per second + 4, + 2, + 0, + [(1, 2, 3, 4), (1, 2, 3, 4), (2, 2, 2, 3), (3, 3, 3, 4)], + ) + + assert test_data.vectors[0] == Vector((1, 2, 3, 4), TimeTuple(10000, 0)) + assert test_data.vectors[1] == Vector((1, 2, 3, 4), TimeTuple(10000, 32768)) + assert test_data.vectors[2] == Vector((2, 2, 2, 3), TimeTuple(10001, 1)) + assert test_data.vectors[3] == Vector((3, 3, 3, 4), TimeTuple(10001, 32769))