Skip to content

Commit fea70ab

Browse files
tech3371tech3371
authored
SWE L1B implementation and updates to L1A (IMAP-Science-Operations-Center#253)
SWE L1B implementation and updates to L1A: - Added project wide utils function and swe specific utils - added logics to populate data based on ESA lookup table number - refactored code to store l1a data as 180x7 instead of 15x12x7 array. - added new event test data and moved tests data into folders. - updated l1a code to be ISTP compliant and to be able to create cdf file - added new test and update test as needed - implemented SWE l1b algorithm for both science and other data types - creates cdf file for SWE l1b data product - added deadtime correction and place holder for in-flight calibration. converted counts to rates - added rest of l1b algorithm to convert raw data to engineering unit. - updated function that write data to cdf to be more generic - changed cdf attrs to use new CDF dataclasses - change energy unit to be unitless for now.
1 parent f234d23 commit fea70ab

56 files changed

Lines changed: 3438 additions & 515 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎imap_processing/cdf/global_attrs.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -250,6 +250,7 @@ class ScienceAttrs(AttrBase):
250250
depend_0: str = None
251251
depend_1: str = None
252252
depend_2: str = None
253+
depend_3: str = None
253254
variable_purpose: str = None
254255
var_notes: str = None
255256

@@ -276,6 +277,9 @@ def output(self):
276277
if self.depend_2 is not None:
277278
endval["DEPEND_2"] = self.depend_2
278279

280+
if self.depend_3 is not None:
281+
endval["DEPEND_3"] = self.depend_3
282+
279283
if self.variable_purpose is not None:
280284
endval["VARIABLE_PURPOSE"] = self.variable_purpose
281285

‎imap_processing/cdf/utils.py‎

Lines changed: 52 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,37 @@
55
from cdflib.xarray import xarray_to_cdf
66

77

8-
def write_cdf(data: xr.Dataset, description: str = "", directory: str = ""):
8+
def calc_start_time(shcoarse_time: int):
9+
"""Calculate the datetime64 from the CCSDS secondary header information.
10+
11+
Since all instrument has SHCOARSE or MET seconds, we need convert it to
12+
UTC. Took this from IDEX code.
13+
14+
Parameters
15+
----------
16+
shcoarse_time: int
17+
Number of seconds since epoch (nominally the launch time)
18+
19+
Returns
20+
-------
21+
np.datetime64
22+
The time of the event
23+
24+
TODO
25+
-----
26+
This conversion is temporary for now, and will need SPICE in the future.
27+
Nick Dutton mentioned that s/c clock start epoch is
28+
jan-1-2010-00:01:06.184 ET
29+
We will use this for now.
30+
"""
31+
# Get the datetime of Jan 1 2010 as the start date
32+
launch_time = np.datetime64("2010-01-01T00:01:06.184")
33+
return launch_time + np.timedelta64(shcoarse_time, "s")
34+
35+
36+
def write_cdf(
37+
data: xr.Dataset, description: str = "", mode: str = "", directory: str = ""
38+
):
939
"""Write the contents of "data" to a CDF file using cdflib.xarray_to_cdf.
1040
1141
This function determines the file name to use from the global attributes,
@@ -20,6 +50,7 @@ def write_cdf(data: xr.Dataset, description: str = "", directory: str = ""):
2050
data (xarray.Dataset): The dataset object to convert to a CDF
2151
description (str): The description to insert into the file name after the
2252
orbit, before the SPICE field. No underscores allowed.
53+
mode (str): Instrument mode
2354
directory (str): The directory to write the file to
2455
2556
Returns
@@ -29,7 +60,17 @@ def write_cdf(data: xr.Dataset, description: str = "", directory: str = ""):
2960
"""
3061
# Determine the start date of the data in the file,
3162
# based on the time of the first dust impact
32-
file_start_date = data["Epoch"][0].data
63+
file_start_date = None
64+
if "idex" in data.attrs["Logical_source"]:
65+
file_start_date = data["Epoch"][0].data
66+
elif "swe" in data.attrs["Logical_source"]:
67+
start_time = data["Epoch"].data[0]
68+
file_start_date = calc_start_time(start_time)
69+
if file_start_date is None:
70+
raise ValueError(
71+
"Unable to determine file start date. Check Logical_source value"
72+
)
73+
3374
date_string = np.datetime_as_string(file_start_date, unit="D").replace("-", "")
3475

3576
# Determine the optional "description" field
@@ -38,13 +79,21 @@ def write_cdf(data: xr.Dataset, description: str = "", directory: str = ""):
3879
if (description.startswith("_") or not description)
3980
else f"_{description}"
4081
)
82+
mode = mode if (mode.startswith("_") or not mode) else f"_{mode}"
4183

4284
# Determine the file name based on the attributes in the xarray
85+
# Set file name based on this convention:
86+
# imap_<instrument>_<datalevel>_<mode>_<descriptor>_<startdate>_
87+
# <version>.cdf
88+
# data.attrs["Logical_source"] has the mission, instrument, and level
89+
# like this:
90+
# imap_idex_l1
4391
filename = (
4492
data.attrs["Logical_source"]
93+
+ mode
94+
+ description
4595
+ "_"
4696
+ date_string
47-
+ description
4897
+ f"_v{data.attrs['Data_version']}.cdf"
4998
)
5099
filename_and_path = os.path.join(directory, filename)

‎imap_processing/idex/tests/test_l1_cdfs.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@ def test_descriptor_in_file_name(decom_test_data, temp_path):
7777
date_to_test = "20250724"
7878
assert file_name == os.path.join(
7979
temp_path,
80-
f"{decom_test_data.data.attrs['Logical_source']}_{date_to_test}_impact-lab-test001_v{idex.__version__}.cdf",
80+
f"{decom_test_data.data.attrs['Logical_source']}_impact-lab-test001_{date_to_test}_v{idex.__version__}.cdf",
8181
)
8282
assert Path(file_name).exists()
8383

‎imap_processing/swe/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
__version__ = "01"

‎imap_processing/swe/l1a/swe_l1a.py‎

Lines changed: 43 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,16 @@
11
import logging
22

3-
from imap_processing.swe.l0 import decom_swe
3+
from imap_processing.cdf.utils import write_cdf
44
from imap_processing.swe.l1a.swe_science import swe_science
5+
from imap_processing.swe.utils.swe_utils import (
6+
SWEAPID,
7+
create_dataset,
8+
filename_descriptors,
9+
)
10+
from imap_processing.utils import group_by_apid, sort_by_time
511

612

7-
def swe_l1a(packet_file: str):
13+
def swe_l1a(packets, cdf_filepath):
814
"""Process SWE l0 data into l1a data.
915
1016
Receive all L0 data file. Based on appId, it
@@ -13,15 +19,39 @@ def swe_l1a(packet_file: str):
1319
1420
Parameters
1521
----------
16-
packet_file : str
17-
The path and filename to the L0 file to read
22+
packets: list
23+
Decom data list that contains all appIds
24+
cdf_filepath: str
25+
Folder path of where to write CDF file
26+
27+
Returns
28+
-------
29+
str
30+
Path name of where CDF file was created.
31+
This is used to upload file from local to s3.
32+
TODO: test this later.
1833
"""
19-
decom_data = decom_swe.decom_packets(packet_file)
20-
logging.info(f"Unpacking data from {packet_file}")
21-
22-
# If appId is science, then the file should contain all data of science appId
23-
if decom_data[0].header["PKT_APID"].raw_value == 1344:
24-
logging.info("Processing science data")
25-
return swe_science(decom_data=decom_data)
26-
else:
27-
return decom_data
34+
# group data by appId
35+
grouped_data = group_by_apid(packets)
36+
37+
for apid in grouped_data.keys():
38+
# If appId is science, then the file should contain all data of science appId
39+
if apid == SWEAPID.SWE_SCIENCE:
40+
# sort data by acquisition time
41+
sorted_packets = sort_by_time(grouped_data[apid], "ACQ_START_COARSE")
42+
logging.debug(
43+
"Processing science data for [%s] packets", len(sorted_packets)
44+
)
45+
data = swe_science(decom_data=sorted_packets)
46+
else:
47+
# If it's not science, we unpack, organize and save it as a dataset.
48+
sorted_packets = sort_by_time(grouped_data[apid], "SHCOARSE")
49+
data = create_dataset(packets=sorted_packets)
50+
51+
# write data to CDF
52+
return write_cdf(
53+
data,
54+
mode=f"{data['APP_MODE'].data[0]}" if apid == SWEAPID.SWE_APP_HK else "",
55+
description=filename_descriptors.get(apid),
56+
directory=cdf_filepath,
57+
)

‎imap_processing/swe/l1a/swe_science.py‎

Lines changed: 84 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,18 @@
11
import collections
2+
import dataclasses
23

34
import numpy as np
45
import xarray as xr
56

7+
from imap_processing.cdf.global_attrs import ConstantCoordinates
8+
from imap_processing.swe import swe_cdf_attrs
9+
from imap_processing.swe.utils.swe_utils import (
10+
add_metadata_to_array,
11+
)
612

7-
def uncompress_counts(cem_count):
8-
"""Uncompress counts from the CEMs.
13+
14+
def decompressed_counts(cem_count):
15+
"""Decompressed counts from the CEMs.
916
1017
Parameters
1118
----------
@@ -15,16 +22,16 @@ def uncompress_counts(cem_count):
1522
Returns
1623
-------
1724
int
18-
uncompressed count. Eg. 40959
25+
decompressed count. Eg. 40959
1926
"""
2027
# index is the first four bits of input data
2128
# multi is the last four bits of input data
2229
index = cem_count // 16
2330
multi = cem_count % 16
2431

2532
# This is look up table for the index to get
26-
# base and step_size to calculate the uncompressed count.
27-
uncompress_table = {
33+
# base and step_size to calculate the decompressed count.
34+
decompress_table = {
2835
0: {"base": 0, "step_size": 1},
2936
1: {"base": 16, "step_size": 1},
3037
2: {"base": 32, "step_size": 2},
@@ -43,36 +50,18 @@ def uncompress_counts(cem_count):
4350
15: {"base": 33792, "step_size": 2048},
4451
}
4552

46-
# uncompression formula from SWE algorithm document CN102D-D0001 and page 16.
53+
# decompression formula from SWE algorithm document CN102D-D0001 and page 16.
4754
# N = base[index] + multi * step_size[index] + (step_size[index] - 1) / 2
4855
# NOTE: for (step_size[index] - 1) / 2, we only keep the whole number part of
4956
# the quotient
5057

5158
return (
52-
uncompress_table[index]["base"]
53-
+ (multi * uncompress_table[index]["step_size"])
54-
+ ((uncompress_table[index]["step_size"] - 1) // 2)
59+
decompress_table[index]["base"]
60+
+ (multi * decompress_table[index]["step_size"])
61+
+ ((decompress_table[index]["step_size"] - 1) // 2)
5562
)
5663

5764

58-
def add_metadata_to_array(data_packet, metadata_arrays):
59-
"""Add metadata to the metadata_arrays.
60-
61-
Parameters
62-
----------
63-
data_packet : space_packet_parser.parser.Packet
64-
SWE data packet
65-
metadata_arrays : dict
66-
metadata arrays
67-
"""
68-
for key, value in data_packet.header.items():
69-
metadata_arrays.setdefault(key, []).append(value.raw_value)
70-
71-
for key, value in data_packet.data.items():
72-
if key != "SCIENCE_DATA":
73-
metadata_arrays.setdefault(key, []).append(value.raw_value)
74-
75-
7665
def swe_science(decom_data):
7766
"""SWE L1a science processing.
7867
@@ -116,7 +105,7 @@ def swe_science(decom_data):
116105

117106
# We know we can only have 8 bit numbers input, so iterate over all
118107
# possibilities once up front
119-
decompression_table = np.array([uncompress_counts(i) for i in range(256)])
108+
decompression_table = np.array([decompressed_counts(i) for i in range(256)])
120109

121110
for data_packet in decom_data:
122111
# read raw data
@@ -129,51 +118,98 @@ def swe_science(decom_data):
129118
# convert bytes to numpy array of uint8
130119
raw_counts = np.frombuffer(byte_data, dtype=np.uint8)
131120

132-
# Uncompress counts. Uncompressed data is a list of 1260
133-
# where 1260 = 15 seconds x 12 energy steps x 7 CEMs
121+
# Uncompress counts. Decompressed data is a list of 1260
122+
# where 1260 = 180 x 7 CEMs
134123
# Take the "raw_counts" indices/counts mapping from
135124
# decompression_table and then reshape the return
136-
uncompress_data = np.take(decompression_table, raw_counts).reshape(15, 12, 7)
137-
raw_counts = raw_counts.reshape(15, 12, 7)
125+
uncompress_data = np.take(decompression_table, raw_counts).reshape(180, 7)
126+
# Save raw counts data as well
127+
raw_counts = raw_counts.reshape(180, 7)
138128

139129
# Save data with its metadata field to attrs and DataArray of xarray.
140-
science_array.append(uncompress_data)
141-
raw_science_array.append(raw_counts)
142-
add_metadata_to_array(data_packet, metadata_arrays)
130+
# Save data as np.int64 to be complaint with ISTP' FILLVAL
131+
science_array.append(uncompress_data.astype(np.int64))
132+
raw_science_array.append(raw_counts.astype(np.int64))
133+
metadata_arrays = add_metadata_to_array(data_packet, metadata_arrays)
143134

144-
met_time = xr.DataArray(
135+
epoch_time = xr.DataArray(
145136
metadata_arrays["SHCOARSE"],
146-
name="met_time",
147-
dims=["met_time"],
148-
attrs=dict(
149-
description="Mission elapsed time",
150-
units="seconds since start of the mission",
151-
),
137+
name="Epoch",
138+
dims=["Epoch"],
139+
attrs=ConstantCoordinates.EPOCH,
140+
)
141+
142+
# TODO: add more descriptive description
143+
energy = xr.DataArray(
144+
np.arange(180),
145+
name="Energy",
146+
dims=["Energy"],
147+
attrs=dataclasses.replace(
148+
swe_cdf_attrs.int_base,
149+
catdesc="Energy's index value in the lookup table",
150+
fieldname="Energy Bins",
151+
label_axis="Energy Bins",
152+
units="",
153+
).output(),
154+
)
155+
156+
counts = xr.DataArray(
157+
np.arange(7),
158+
name="Counts",
159+
dims=["Counts"],
160+
attrs=dataclasses.replace(
161+
swe_cdf_attrs.int_base,
162+
catdesc="Counts",
163+
fieldname="Counts",
164+
label_axis="Counts",
165+
units="int",
166+
).output(),
152167
)
153168

154169
science_xarray = xr.DataArray(
155170
science_array,
156-
dims=["met_time", "seconds", "energy_steps", "cem_counts"],
171+
dims=["Epoch", "Energy", "Counts"],
172+
attrs=swe_cdf_attrs.l1a_science_attrs.output(),
157173
)
174+
158175
raw_science_xarray = xr.DataArray(
159176
raw_science_array,
160-
dims=["met_time", "seconds", "energy_steps", "cem_counts"],
177+
dims=["Epoch", "Energy", "Counts"],
178+
attrs=swe_cdf_attrs.l1a_science_attrs.output(),
161179
)
162180

163181
dataset = xr.Dataset(
164-
{"SCIENCE_DATA": science_xarray},
165-
coords={"met_time": met_time},
182+
coords={
183+
"Epoch": epoch_time,
184+
"Energy": energy,
185+
"Counts": counts,
186+
},
187+
attrs=swe_cdf_attrs.swe_l1a_global_attrs.output(),
166188
)
167-
189+
dataset["SCIENCE_DATA"] = science_xarray
168190
dataset["RAW_SCIENCE_DATA"] = raw_science_xarray
169191

170192
# create xarray dataset for each metadata field
171193
for key, value in metadata_arrays.items():
172194
if key == "SHCOARSE":
173195
continue
196+
# TODO: figure out how to add more descriptive
197+
# description for each metadata field
198+
#
199+
# int_attrs["CATDESC"] = int_attrs["FIELDNAM"] = int_attrs["LABLAXIS"] = key
200+
# # get int32's max since most of metadata is under 32-bits
201+
# int_attrs["VALIDMAX"] = np.iinfo(np.int32).max
202+
# int_attrs["DEPEND_0"] = "Epoch"
174203
dataset[key] = xr.DataArray(
175204
value,
176-
dims=["met_time"],
205+
dims=["Epoch"],
206+
attrs=dataclasses.replace(
207+
swe_cdf_attrs.swe_metadata_attrs,
208+
catdesc=key,
209+
fieldname=key,
210+
label_axis=key,
211+
depend_0="Epoch",
212+
).output(),
177213
)
178214

179215
return dataset

0 commit comments

Comments
 (0)