Skip to content

Commit 404d053

Browse files
committed
Initial commit mostly from original #212
But with ruff treatment.
1 parent 2eabdb5 commit 404d053

File tree

3 files changed

+244
-0
lines changed

3 files changed

+244
-0
lines changed

cheta/fetch.py

+19
Original file line numberDiff line numberDiff line change
@@ -1093,6 +1093,25 @@ def copy(self):
10931093

10941094
return deepcopy(self)
10951095

1096+
def set_hi_res_times(self, fmt_intervals=None):
1097+
"""Update MSID timestamps to have minor frame time resolution.
1098+
This makes the MSID ``times`` match MAUDE. (TODO: better explanation)
1099+
"""
1100+
from cheta.time_offsets import get_hi_res_times
1101+
1102+
if self.stat is not None:
1103+
return
1104+
1105+
data_sources = list(self.data_source.keys())
1106+
if data_sources == ["maude"]:
1107+
return
1108+
elif data_sources == ["cxc"]:
1109+
self.times, _ = get_hi_res_times(self, fmt_intervals)
1110+
else:
1111+
raise ValueError(
1112+
"cannot set hi_res times for query with mixed cxc and maude data"
1113+
)
1114+
10961115
def filter_bad(self, bads=None, copy=False):
10971116
"""Filter out any bad values.
10981117

cheta/tests/test_time_offsets.py

+53
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
import numpy as np
2+
3+
from cheta import fetch
4+
from cheta.time_offsets import MNF_TIME, adjust_time, get_hi_res_times
5+
6+
7+
def test_time_adjust_adjust_time():
8+
cpa2pwr = fetch.Msid("cpa2pwr", "2020:230:00:00:01", "2020:230:00:00:04")
9+
cpa2pwr_adj = adjust_time(cpa2pwr, "2020:230:00:00:01", "2020:230:00:00:04")
10+
for t, t_adj in zip(cpa2pwr.times, cpa2pwr_adj.times):
11+
assert np.isclose(t_adj - t, MNF_TIME) # MNF offset is 1 for this MSID
12+
13+
14+
def test_time_adjust_get_hi_res_times():
15+
# Cover an interval that has formats 1, 2, and 4.
16+
start = "2020:030:03:00:00"
17+
stop = "2020:030:05:00:00"
18+
dat = fetch.Msid("tephin", start, stop)
19+
20+
dat_adj = adjust_time(dat, start, stop)
21+
times_adj, fmt_intervals = get_hi_res_times(dat)
22+
23+
# Test re-using fmt_intervals
24+
times_adj2, fmt_intervals = get_hi_res_times(dat, fmt_intervals)
25+
26+
assert np.all(times_adj == dat_adj.times)
27+
assert np.all(times_adj2 == times_adj)
28+
29+
offsets = times_adj - dat.times
30+
assert np.allclose(offsets[:115], 59 * MNF_TIME) # Format 2 MNF offset is 59
31+
assert np.allclose(offsets[115:1011], 0 * MNF_TIME) # Format 4 MNF offset is 0
32+
assert np.allclose(offsets[1011:], 59 * MNF_TIME) # Format 1 MNF offset is 59
33+
34+
35+
def test_time_adjust_set_hi_res_times():
36+
start, stop = "2020:030:03:00:00", "2020:030:05:00:00"
37+
tephin = fetch.Msid("tephin", start, stop)
38+
times_adj, _ = get_hi_res_times(tephin)
39+
40+
tephin.set_hi_res_times()
41+
assert np.all(tephin.times == times_adj)
42+
43+
44+
def test_time_adjust_msidset_set_hi_res_times():
45+
start, stop = "2020:030:03:00:00", "2020:030:05:00:00"
46+
47+
msids = fetch.Msidset(["tephin", "tcylaft6"], start, stop)
48+
msids.set_hi_res_times()
49+
50+
for msid_name in ("tephin", "tcylaft6"):
51+
msid = fetch.Msid(msid_name, start, stop)
52+
times_adj, _ = get_hi_res_times(msid)
53+
assert np.all(msids[msid_name].times == times_adj)

cheta/time_offsets.py

+172
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
# Licensed under a 3-clause BSD style license - see LICENSE.rst
2+
"""
3+
Describe what this module does.
4+
"""
5+
6+
import copy
7+
8+
import numpy as np
9+
from Ska.engarchive import fetch_eng as fetch
10+
from Ska.engarchive.utils import logical_intervals, state_intervals
11+
from Ska.tdb import msids as tdb_msids
12+
from Ska.tdb import tables
13+
14+
MNF_TIME = 0.25625 # Minor frame time (s)
15+
16+
17+
def adjust_time(msid, start, stop):
18+
"""
19+
Apply an offset to cheta timestamps to achieve minor frame time resolution
20+
21+
Given an msid returned by fetch.MSID or fetch.Msid and a time interval,
22+
determine the time offset for each sample. Time offset is based on the TDB
23+
and the Telemetry format at the time. Only applies to msids retrieved with
24+
'cxc' as a data source
25+
26+
Notes/Improvements:
27+
28+
- TBD: Should accept MSIDSet in addition to MSID as input, but this could
29+
interfere with the current MSIDSet.interpolate() behavior
30+
- TBD: add an option for in-place adjustment
31+
- TBD: - maybe take msid as string, but manually set data source to avoid
32+
'maude' 'cxc' differences...
33+
- NOTE: Could use times from the msid.times span, but may create problems at
34+
the beginning and end since these are subject to adjustment. Should use
35+
time query that generated the MSID
36+
37+
:param msid: fetch ``MSID`` or ``Msid`` object
38+
MSID for which to calculate the time adjustments statistics for.
39+
:param start: str, ``CxoTime``, ``DateTime``
40+
Start time (CxoTime-compatible format)
41+
:param stop: str, ``CxoTime``, ``DateTime``
42+
Stop time (CxoTime-compatible format)
43+
44+
:returns: fetch ``MSID`` or ``Msid`` object
45+
"""
46+
samp_rate = tables["tsmpl"][msid.msid]["SAMPLE_RATE"]
47+
str_num = tables["tsmpl"][msid.msid]["STREAM_NUMBER"]
48+
start_minor_frame = tables["tloc"][msid.msid]["START_MINOR_FRAME"]
49+
# stream# to FMT name converter. Unsure if FMT6/SHuttle data needs to be handled
50+
fmt_dict = {
51+
1: "FMT1",
52+
2: "FMT2",
53+
3: "FMT3",
54+
4: "FMT4",
55+
5: "FMT5",
56+
72: "FMT6",
57+
6: "FMT6",
58+
}
59+
# for each format, generate time offset based on stream number, sample
60+
# rate and start minor frame...
61+
t_off = {}
62+
t_samp = {}
63+
str_num_2_idx = {}
64+
idx = 0
65+
66+
# Build stream number to index decoder ring. For each MSID there will be N
67+
# stream number entries in the table. There's probably a list comprehension
68+
# way to do this
69+
for stream in str_num:
70+
str_num_2_idx[stream] = idx
71+
idx += 1
72+
73+
for stream in str_num:
74+
# Now calculate the offset and sample rate for each stream
75+
# off_secs = 128/samp_rate -1 # magnitude of offset. For SR = 128 it's 0,
76+
# for SR = 64 it's 0,1, for SR = 32, it's 0,1,2,3 &c.
77+
fmt = fmt_dict[stream]
78+
# Note this will fail for very long data sets. Should use average mission
79+
# rate or clock look up table.
80+
off_phase = MNF_TIME * start_minor_frame[str_num_2_idx[stream]]
81+
t_off[fmt] = off_phase
82+
# 128 -> 0.25625, 64 -> 0.51250 , &c.
83+
t_samp[fmt] = (128 / samp_rate[str_num_2_idx[stream]]) * 0.25625
84+
85+
# Get Telemetry format for the time interval in question. CCSDSTMF is
86+
# updated 128x per MjF, so it's at max time resolution
87+
tmf = fetch.Msid("CCSDSTMF", start, stop)
88+
89+
# Generate list of intervals for each format using logical intervals
90+
fmts = ("FMT1", "FMT2", "FMT3", "FMT4", "FMT5", "FMT6")
91+
tmf_intervals = {}
92+
for fmt in fmts:
93+
tmf_intervals[fmt] = logical_intervals(
94+
tmf.times, tmf.vals == fmt, complete_intervals=False
95+
)
96+
97+
# Make scratchpad copy to avoid in-place effects (I think need to check by
98+
# ref or by val convention)
99+
times = msid.times.copy()
100+
for fmt in fmts:
101+
for interval in tmf_intervals[fmt]:
102+
# Now traverse each interval in the msid to be adjusted and add the appropriate offset.
103+
times[
104+
(msid.times >= interval["tstart"]) & (msid.times < interval["tstop"])
105+
] += t_off[fmt] # not positive on >= convention
106+
107+
# This is abusing the MSID class a bit. Not sure how to create 'empty' MSID object
108+
out_msid = copy.deepcopy(msid)
109+
out_msid.times = times
110+
return out_msid
111+
112+
113+
def get_hi_res_times(msid, fmt_intervals=None):
114+
"""
115+
Determine MSID timestamps to achieve minor frame time resolution.
116+
117+
Given an msid returned by fetch.MSID or fetch.Msid determine the time offset
118+
for each sample. Time offset is based on the TDB and the Telemetry format
119+
at the time. Only applies to msids retrieved with 'cxc' as a data source
120+
121+
:param msid: ``MSID`` or ``Msid`` object
122+
MSID for which to calculate the time adjustments statistics.
123+
:param fmt_intervals: TBD
124+
125+
:returns: (np.array, TDB)
126+
Return tuple of (hi-res times, telemetry format intervals)
127+
"""
128+
# TODO: check for data_source
129+
# TODO: check for stats MSID
130+
131+
try:
132+
tdb_msid = tdb_msids[msid.msid]
133+
except KeyError:
134+
raise ValueError(f"msid {msid.msid} is not in TDB")
135+
136+
# stream# to FMT name converter. Unsure if FMT6/SHuttle data needs to be handled
137+
fmt_dict = {
138+
1: "FMT1",
139+
2: "FMT2",
140+
3: "FMT3",
141+
4: "FMT4",
142+
5: "FMT5",
143+
72: "FMT6",
144+
6: "FMT6",
145+
}
146+
147+
# For each format, generate time offset based on stream number
148+
t_off = {}
149+
for stream, start_minor_frame in zip(
150+
tdb_msid.Tloc["STREAM_NUMBER"], tdb_msid.Tloc["START_MINOR_FRAME"]
151+
):
152+
fmt = fmt_dict[stream]
153+
# Note this will fail for very long data sets. Should use average mission
154+
# rate or clock look up table.
155+
t_off[fmt] = MNF_TIME * start_minor_frame
156+
157+
if fmt_intervals is None:
158+
msid_fmt = fetch.Msid("CCSDSTMF", msid.tstart, msid.tstop)
159+
fmt_intervals = state_intervals(msid_fmt.times, msid_fmt.vals)
160+
161+
# This gives start/stop indices that are equivalent to:
162+
# (msid.times >= fmt_interval['tstart']) & (msid.times < fmt_interval['tstop'])
163+
# Note see: Boolean masks and np.searchsorted() in
164+
# https://occweb.cfa.harvard.edu/twiki/bin/view/Aspect/SkaPython#Ska_idioms_and_style
165+
i_starts = np.searchsorted(msid.times, fmt_intervals["tstart"])
166+
i_stops = np.searchsorted(msid.times, fmt_intervals["tstop"])
167+
168+
times = msid.times.copy()
169+
for i_start, i_stop, fmt in zip(i_starts, i_stops, fmt_intervals["val"]):
170+
times[i_start:i_stop] += t_off[fmt]
171+
172+
return times, fmt_intervals

0 commit comments

Comments
 (0)