From 113394cfa7635436edb9aaf5ccd2b42e395243dd Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Fri, 11 Oct 2024 08:02:14 -0400 Subject: [PATCH 01/21] improvement to photon fits files --- padre_meddea/calibration/calibration.py | 16 +++++++++++++--- padre_meddea/data/fits_keywords_primaryhdu.csv | 2 ++ padre_meddea/io/file_tools.py | 6 +++--- 3 files changed, 18 insertions(+), 6 deletions(-) diff --git a/padre_meddea/calibration/calibration.py b/padre_meddea/calibration/calibration.py index a2cc766..32f90b9 100644 --- a/padre_meddea/calibration/calibration.py +++ b/padre_meddea/calibration/calibration.py @@ -49,20 +49,30 @@ def process_file(filename: Path, overwrite=False) -> list: if filename.suffix == ".bin": parsed_data = read_raw_file(filename) if parsed_data["photons"] is not None: # we have event list data - ph_list = parsed_data["photons"] + event_list = parsed_data["photons"] + hdu = fits.PrimaryHDU(data=None) + # fill in metadata hdu.header["DATE"] = (Time.now().fits, "FITS file creation date in UTC") + hdu.header["DATE-BEG"] = (event_list.time[0].fits, "Acquisition start time") + hdu.header["DATE-END"] = (event_list.time[-1].fits, "Acquisition end time") + hdu.header["DATE-AVG"] = (event_list.time[int(len(event_list)/2.)].fits, "Average time of acquisition") + hdu.header["DATEREF"] = (event_list.time[0].fits, "Reference date") + hdu.header["DSUN_AU"] = 1 + hdu.header["LEVEL"] = 1 + # add common fits keywords fits_meta = read_fits_keyword_file( padre_meddea._data_directory / "fits_keywords_primaryhdu.csv" ) for row in fits_meta: hdu.header[row["keyword"]] = (row["value"], row["comment"]) - bin_hdu = fits.BinTableHDU(data=Table(ph_list)) + + bin_hdu = fits.BinTableHDU(data=Table(event_list)) hdul = fits.HDUList([hdu, bin_hdu]) path = create_science_filename( "meddea", - ph_list["time"][0].fits, + event_list["time"][0].fits, "l1", descriptor="eventlist", test=True, diff --git a/padre_meddea/data/fits_keywords_primaryhdu.csv b/padre_meddea/data/fits_keywords_primaryhdu.csv index 1860e92..74223dc 100644 --- a/padre_meddea/data/fits_keywords_primaryhdu.csv +++ b/padre_meddea/data/fits_keywords_primaryhdu.csv @@ -1,6 +1,8 @@ keyword,value,comment AUTHOR,Steven D. Christe,Who designed the observation CREATOR,padre_meddea,Name of software pipeline that produced the FITS file +DETECTOR,meddea,Name of the detector +INFO_URL,https://padre-meddea.readthedocs.io/en/latest/user-guide/data.html,a human-readable web page describing the data OBSRVTRY,PADRE,Name of the observatory TIMESYS,UTC,Time scale of the time-related keywords TELESCOP,PADRE/MeDDEA,Telescope/Sensor name diff --git a/padre_meddea/io/file_tools.py b/padre_meddea/io/file_tools.py index 7c48ea8..f8d2332 100644 --- a/padre_meddea/io/file_tools.py +++ b/padre_meddea/io/file_tools.py @@ -162,10 +162,10 @@ def parse_ph_packets(filename: Path): time=ph_times, data={ "atod": hit_list[3, :], - "asic_num": hit_list[1, :], - "asic_channel": hit_list[2, :], + "asic": hit_list[1, :], + "channel": hit_list[2, :], "clock": hit_list[5, :], - "packet_num": hit_list[4, :], + "pktnum": hit_list[4, :], # "num": np.ones(len(hit_list[2, :])), }, meta={ From 0dd86b85d4a02e4d9bad7aed53fb5ec6bd5219ae Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Wed, 16 Oct 2024 13:35:02 -0400 Subject: [PATCH 02/21] add test for baseline --- padre_meddea/util/util.py | 36 +++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/padre_meddea/util/util.py b/padre_meddea/util/util.py index 2752c03..8bb407e 100644 --- a/padre_meddea/util/util.py +++ b/padre_meddea/util/util.py @@ -5,11 +5,16 @@ import os from datetime import datetime, timezone import time +from pathlib import Path +import numpy as np from astropy.time import Time +from ccsdspy.utils import split_packet_bytes, split_by_apid + __all__ = [ "create_science_filename", + "has_baseline" ] TIME_FORMAT = "%Y%m%dT%H%M%S" @@ -18,7 +23,6 @@ def create_science_filename( - instrument: str, time: str, level: str, version: str, @@ -99,3 +103,33 @@ def create_science_filename( filename = filename.replace("__", "_") # reformat if mode or descriptor not given return filename + FILENAME_EXTENSION + + +def has_baseline(filename: Path, packet_count=10): + """Given a stream of photon packets, check whether the baseline measurement is included. + Baseline packets have one extra word per photon for a total of 4 words (8 bytes). + + This function calculates the number of hits in the packet assuming 4 words per photon. + If the resultant is not an integer number then returns False. + + Parameters + ---------- + packet_bytes : byte string + Photon packet bytes, must be an integer number of whole packets and greaterh + + Returns + ------- + result : bool + + """ + HEADER_BYTES = 11 * 16 / 8 + BYTES_PER_PHOTON = 16 * 4 / 8 + + with open(filename, "rb") as mixed_file: + stream_by_apid = split_by_apid(mixed_file) + packet_stream = stream_by_apid[160] + packet_bytes = split_packet_bytes(packet_stream) + num_hits = np.zeros(packet_count) + for i in range(packet_count): + num_hits[i] = (len(packet_bytes[i]) - HEADER_BYTES) / BYTES_PER_PHOTON + return np.sum(num_hits - np.floor(num_hits)) == 1 \ No newline at end of file From 95f6a4132b89a4e6bf47dfafb27947959e3a479c Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Wed, 16 Oct 2024 13:35:27 -0400 Subject: [PATCH 03/21] baseline aware file processing --- padre_meddea/io/file_tools.py | 113 ++++++++++++++++++++++------------ 1 file changed, 72 insertions(+), 41 deletions(-) diff --git a/padre_meddea/io/file_tools.py b/padre_meddea/io/file_tools.py index f8d2332..1928fa7 100644 --- a/padre_meddea/io/file_tools.py +++ b/padre_meddea/io/file_tools.py @@ -9,6 +9,8 @@ import astropy.units as u from astropy.timeseries import TimeSeries from astropy.io import ascii +from astropy.table import Table +from astropy.time import Time import ccsdspy from ccsdspy import PacketField, PacketArray @@ -18,6 +20,7 @@ ) import padre_meddea +from padre_meddea.util.util import has_baseline __all__ = ["read_file", "read_raw_file"] @@ -122,20 +125,42 @@ def parse_ph_packets(filename: Path): with open(filename, "rb") as mixed_file: stream_by_apid = split_by_apid(mixed_file) packet_stream = stream_by_apid.get(APID["photon"], None) + print(packet_stream) if packet_stream is None: return None packet_definition = packet_definition_ph() pkt = ccsdspy.VariableLength(packet_definition) ph_data = pkt.load(packet_stream, include_primary_header=True) - integration_time = ph_data["INTEGRATION_TIME"] * 12.8 * 1e-6 - live_time = ph_data["LIVE_TIME"] * 12.8 * 1e-6 - dead_time = (1.0 - live_time / integration_time) * 100 packet_time = ph_data["TIME_S"] + ph_data["TIME_CLOCKS"] / 20.0e6 + if has_baseline(filename): + WORDS_PER_HIT = 4 + else: + WORDS_PER_HIT = 3 + + pkt_list = Table() + pkt_list['seqcount'] = ph_data["CCSDS_SEQUENCE_COUNT"] + pkt_list['time_s'] = ph_data["TIME_S"] + pkt_list['clocks'] = ph_data["TIME_CLOCKS"] + pkt_list['livetime'] = ph_data["LIVE_TIME"] + pkt_list['inttime'] = ph_data["INTEGRATION_TIME"] + total_hits = 0 for this_packet in ph_data["PIXEL_DATA"]: - total_hits += len(this_packet[0::3]) + total_hits += len(this_packet[0::WORDS_PER_HIT]) + + if (total_hits - np.floor(total_hits)) != 0: + raise ValueError(f"Got non-integer number of hits {total_hits - np.floor(total_hits)}.") - hit_list = np.zeros((6, total_hits), dtype="uint16") + # hit_list definition + # 0, raw id number + # 1, asic number, 0 to 7 + # 2, channel number, 0 to 32 + # 3, energy 12 bits + # 4, packet sequence number 12 bits + # 5, clock number + # 6, baseline (if exists) + + hit_list = np.zeros((7, total_hits), dtype="uint16") time_stamps = np.zeros(total_hits) # 0 time, 1 asic_num, 2 channel num, 3 hit channel i = 0 @@ -148,41 +173,47 @@ def parse_ph_packets(filename: Path): ids = this_ph_data[1::3] asic_num = (ids & 0b11100000) >> 5 channel_num = ids & 0b00011111 - hit_list[0, i : i + num_hits] = this_ph_data[0::3] - time_stamps[i : i + num_hits] = this_ph_data[0::3] * 12.8e-6 + this_time + hit_list[0, i : i + num_hits] = this_ph_data[0::WORDS_PER_HIT] + time_stamps[i : i + num_hits] = this_ph_data[0::WORDS_PER_HIT] * 12.8e-6 + this_time hit_list[1, i : i + num_hits] = asic_num hit_list[2, i : i + num_hits] = channel_num - hit_list[3, i : i + num_hits] = this_ph_data[2::3] hit_list[4, i : i + num_hits] = this_pkt_num - hit_list[5, i : i + num_hits] = this_ph_data[0::3] + hit_list[5, i : i + num_hits] = this_ph_data[0::WORDS_PER_HIT] + if WORDS_PER_HIT == 4: + hit_list[6, i : i + num_hits] = this_ph_data[2::WORDS_PER_HIT] # baseline + hit_list[3, i : i + num_hits] = this_ph_data[3::WORDS_PER_HIT] # hit energy + elif WORDS_PER_HIT == 3: + hit_list[3, i : i + num_hits] = this_ph_data[2::WORDS_PER_HIT] # hit energy i += num_hits - ph_times = [EPOCH + dt.timedelta(seconds=this_t) for this_t in time_stamps] - ph_list = TimeSeries( - time=ph_times, - data={ - "atod": hit_list[3, :], - "asic": hit_list[1, :], - "channel": hit_list[2, :], - "clock": hit_list[5, :], - "pktnum": hit_list[4, :], - # "num": np.ones(len(hit_list[2, :])), - }, - meta={ - "integration_times": integration_time, - "live_times": live_time, - "dead_times": dead_time, - "total_hits": total_hits, - "packet_time": packet_time, - "time_s": ph_data["TIME_S"], - "time_clocks": ph_data["TIME_CLOCKS"], - }, - ) - ph_list.sort() - int_time = (ph_list.time.max() - ph_list.time.min()).to("min") - ph_list.meta.update({"int_time": int_time}) - ph_list.meta.update({"avg rate": (total_hits * u.ct) / int_time}) - return ph_list + event_list = Table() + event_list["seqcount"] = hit_list[4, :] + event_list["clock"] = hit_list[5, :] + event_list["asic"] = hit_list[1, :].astype(np.uint8) + event_list["channel"] = hit_list[2, :].astype(np.uint8) + event_list["atod"] = hit_list[3, :] + event_list["baseline"] = hit_list[6, :] # if baseline not present then all zeros + + ph_times = Time([EPOCH + dt.timedelta(seconds=this_t) for this_t in time_stamps]) + event_list.meta.update({'DATE-BEG': ph_times.min().fits}) + event_list.meta.update({'DATE-END': ph_times.max().fits}) + event_list.meta.update({'DATE-AVG': ph_times[int(len(ph_times)/2.)].fits}) + #event_list = TimeSeries( + # time=ph_times, + # data={ + # "atod": hit_list[3, :], + # "asic": hit_list[1, :], + # "channel": hit_list[2, :], + # "clock": hit_list[5, :], + # "pktnum": hit_list[4, :], + # # "num": np.ones(len(hit_list[2, :])), + # } + #) + #event_list.sort() + # int_time = (event_list.time.max() - event_list.time.min()).to("min") + # event_list.meta.update({"int_time": int_time}) + # event_list.meta.update({"avg rate": (total_hits * u.ct) / int_time}) + return event_list, pkt_list def parse_hk_packets(filename: Path): @@ -200,12 +231,12 @@ def parse_hk_packets(filename: Path): """ with open(filename, "rb") as mixed_file: stream_by_apid = split_by_apid(mixed_file) - packet_stream = stream_by_apid.get(APID["housekeeping"], None) - if packet_stream is None: + packet_bytes = stream_by_apid.get(APID["housekeeping"], None) + if packet_bytes is None: return None packet_definition = packet_definition_hk() pkt = ccsdspy.FixedLength(packet_definition) - hk_data = pkt.load(packet_stream) + hk_data = pkt.load(packet_bytes) hk_timestamps = [ dt.timedelta(seconds=int(this_t)) + EPOCH for this_t in hk_data["TIMESTAMP"] ] @@ -228,12 +259,12 @@ def parse_spectrum_packets(filename: Path): """ with open(filename, "rb") as mixed_file: stream_by_apid = split_by_apid(mixed_file) - packet_stream = stream_by_apid.get(APID["spectrum"], None) - if packet_stream is None: + packet_bytes = stream_by_apid.get(APID["spectrum"], None) + if packet_bytes is None: return None packet_definition = packet_definition_hist2() pkt = ccsdspy.FixedLength(packet_definition) - data = pkt.load(packet_stream) + data = pkt.load(packet_bytes) timestamps = [ dt.timedelta(seconds=int(this_t)) + EPOCH for this_t in data["TIMESTAMPS"] ] From 44ba969026f78c7fedf57e69c9d24a7204f9b2dd Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Wed, 16 Oct 2024 13:35:49 -0400 Subject: [PATCH 04/21] adding more meta data to file files and new organization --- padre_meddea/calibration/calibration.py | 57 +++++++++++++++++-------- pyproject.toml | 1 + 2 files changed, 40 insertions(+), 18 deletions(-) diff --git a/padre_meddea/calibration/calibration.py b/padre_meddea/calibration/calibration.py index 32f90b9..426538b 100644 --- a/padre_meddea/calibration/calibration.py +++ b/padre_meddea/calibration/calibration.py @@ -11,6 +11,7 @@ from astropy.table import Table from swxsoc.util.util import record_timeseries +import git import padre_meddea from padre_meddea import log @@ -45,35 +46,55 @@ def process_file(filename: Path, overwrite=False) -> list: # Check if the LAMBDA_ENVIRONMENT environment variable is set lambda_environment = os.getenv("LAMBDA_ENVIRONMENT") output_files = [] + file_path = Path(filename) - if filename.suffix == ".bin": - parsed_data = read_raw_file(filename) + if file_path.suffix == ".bin": + parsed_data = read_raw_file(file_path) if parsed_data["photons"] is not None: # we have event list data - event_list = parsed_data["photons"] + event_list, pkt_list = parsed_data["photons"] + + primary_hdr = fits.Header() - hdu = fits.PrimaryHDU(data=None) # fill in metadata - hdu.header["DATE"] = (Time.now().fits, "FITS file creation date in UTC") - hdu.header["DATE-BEG"] = (event_list.time[0].fits, "Acquisition start time") - hdu.header["DATE-END"] = (event_list.time[-1].fits, "Acquisition end time") - hdu.header["DATE-AVG"] = (event_list.time[int(len(event_list)/2.)].fits, "Average time of acquisition") - hdu.header["DATEREF"] = (event_list.time[0].fits, "Reference date") - hdu.header["DSUN_AU"] = 1 - hdu.header["LEVEL"] = 1 + primary_hdr["DATE"] = (Time.now().fits, "FITS file creation date in UTC") + for this_keyword, this_str in zip(["DATE-BEG", "DATE-END", "DATE-AVG"], ["Acquisition start time", "Acquisition end time", "Average time of acquisition"]): + primary_hdr[this_keyword] = (event_list.meta.get(this_keyword, ""), this_str) + + primary_hdr["DATEREF"] = (primary_hdr["DATE-BEG"], "Reference date") + primary_hdr["LEVEL"] = 1 + + # add processing information + primary_hdr["PRSTEP1"] = ("PROCESS Raw to L1", "Processing step type") + primary_hdr["PRPROC1"] = ("padre_meddea.calibration.process", "Name of procedure performing PRSTEP1") + primary_hdr["PRPVER1"] = (padre_meddea.__version__, "Version of procedure PRPROC1") + primary_hdr["PRLIB1A"] = ("padre_meddea", "Software library containing PRPROC1") + primary_hdr["PRVER1A"] = (padre_meddea.__version__, "Version of PRLIB1A") + repo = git.Repo(search_parent_directories=True) + primary_hdr["PRHSH1A"] = (repo.head.object.hexsha, "GIT commit hash for PRLIB1A") + primary_hdr["PRBRA1A"] = (repo.active_branch.name, "GIT/SVN repository branch of PRLIB1A") + commits = list(repo.iter_commits("main", max_count=1)) + primary_hdr["PRVER1B"] = (Time(commits[0].committed_datetime).fits, "Date of last commit of PRLIB1B") + # primary_hdr["PRLOG1"] add log information, need to do this after the fact + # primary_hdr["PRENV1"] add information about processing env, need to do this after the fact + + # add common fits keywords fits_meta = read_fits_keyword_file( padre_meddea._data_directory / "fits_keywords_primaryhdu.csv" ) for row in fits_meta: - hdu.header[row["keyword"]] = (row["value"], row["comment"]) - - bin_hdu = fits.BinTableHDU(data=Table(event_list)) - hdul = fits.HDUList([hdu, bin_hdu]) + primary_hdr[row["keyword"]] = (row["value"], row["comment"]) + + empty_primary = fits.PrimaryHDU(header=primary_hdr) + pkt_hdu = fits.BinTableHDU(pkt_list, name="PKT") + pkt_hdu.add_checksum() + hit_hdu = fits.BinTableHDU(event_list, name="SCI") + hit_hdu.add_checksum() + hdul = fits.HDUList([empty_primary, hit_hdu, pkt_hdu]) path = create_science_filename( - "meddea", - event_list["time"][0].fits, - "l1", + time=primary_hdr["DATE-BEG"], + level="l1", descriptor="eventlist", test=True, version="0.1.0", diff --git a/pyproject.toml b/pyproject.toml index 4d7262d..7d934b9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,6 +33,7 @@ dependencies = [ 'ccsdspy==1.3.*', 'matplotlib==3.9.*', 'swxsoc @ git+https://github.com/swxsoc/swxsoc.git@main', + 'gitpython==3.1.* ] [project.optional-dependencies] From a8e2ec54a03d9bbd8e3c94a696db5213a1279f3a Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Thu, 17 Oct 2024 18:19:15 -0400 Subject: [PATCH 05/21] adding metadata and reading fits files --- padre_meddea/__init__.py | 10 +++ padre_meddea/calibration/calibration.py | 50 +++++++++-- padre_meddea/io/file_tools.py | 106 ++++++++++++++++-------- padre_meddea/util/util.py | 58 +++++++++++-- pyproject.toml | 2 +- 5 files changed, 173 insertions(+), 53 deletions(-) diff --git a/padre_meddea/__init__.py b/padre_meddea/__init__.py index 11d246c..90ef540 100644 --- a/padre_meddea/__init__.py +++ b/padre_meddea/__init__.py @@ -9,6 +9,7 @@ version_tuple = (0, 0, "unknown version") from swxsoc import config as swxsoc_config, log as swxsoc_log, print_config +from astropy.time import Time # Load user configuration config = swxsoc_config @@ -54,4 +55,13 @@ 10.73, ] +APID = { + "spectrum": 0xA2, # decimal 162 + "photon": 0xA0, # decimal 160 + "housekeeping": 0xA3, # decimal 163 + "cmd_resp": 0xA5, # decimal 165 +} + +EPOCH = Time("2000-01-01 00:00", scale="utc") + log.debug(f"padre_meddea version: {__version__}") diff --git a/padre_meddea/calibration/calibration.py b/padre_meddea/calibration/calibration.py index 426538b..d74ba7b 100644 --- a/padre_meddea/calibration/calibration.py +++ b/padre_meddea/calibration/calibration.py @@ -57,26 +57,58 @@ def process_file(filename: Path, overwrite=False) -> list: # fill in metadata primary_hdr["DATE"] = (Time.now().fits, "FITS file creation date in UTC") - for this_keyword, this_str in zip(["DATE-BEG", "DATE-END", "DATE-AVG"], ["Acquisition start time", "Acquisition end time", "Average time of acquisition"]): - primary_hdr[this_keyword] = (event_list.meta.get(this_keyword, ""), this_str) + for this_keyword, this_str in zip( + ["DATE-BEG", "DATE-END", "DATE-AVG"], + [ + "Acquisition start time", + "Acquisition end time", + "Average time of acquisition", + ], + ): + primary_hdr[this_keyword] = ( + event_list.meta.get(this_keyword, ""), + this_str, + ) primary_hdr["DATEREF"] = (primary_hdr["DATE-BEG"], "Reference date") - primary_hdr["LEVEL"] = 1 + primary_hdr["LEVEL"] = (0, "Data level of fits file") # add processing information primary_hdr["PRSTEP1"] = ("PROCESS Raw to L1", "Processing step type") - primary_hdr["PRPROC1"] = ("padre_meddea.calibration.process", "Name of procedure performing PRSTEP1") - primary_hdr["PRPVER1"] = (padre_meddea.__version__, "Version of procedure PRPROC1") - primary_hdr["PRLIB1A"] = ("padre_meddea", "Software library containing PRPROC1") + primary_hdr["PRPROC1"] = ( + "padre_meddea.calibration.process", + "Name of procedure performing PRSTEP1", + ) + primary_hdr["PRPVER1"] = ( + padre_meddea.__version__, + "Version of procedure PRPROC1", + ) + primary_hdr["PRLIB1A"] = ( + "padre_meddea", + "Software library containing PRPROC1", + ) primary_hdr["PRVER1A"] = (padre_meddea.__version__, "Version of PRLIB1A") repo = git.Repo(search_parent_directories=True) - primary_hdr["PRHSH1A"] = (repo.head.object.hexsha, "GIT commit hash for PRLIB1A") - primary_hdr["PRBRA1A"] = (repo.active_branch.name, "GIT/SVN repository branch of PRLIB1A") + primary_hdr["PRHSH1A"] = ( + repo.head.object.hexsha, + "GIT commit hash for PRLIB1A", + ) + primary_hdr["PRBRA1A"] = ( + repo.active_branch.name, + "GIT/SVN repository branch of PRLIB1A", + ) commits = list(repo.iter_commits("main", max_count=1)) - primary_hdr["PRVER1B"] = (Time(commits[0].committed_datetime).fits, "Date of last commit of PRLIB1B") + primary_hdr["PRVER1B"] = ( + Time(commits[0].committed_datetime).fits, + "Date of last commit of PRLIB1B", + ) # primary_hdr["PRLOG1"] add log information, need to do this after the fact # primary_hdr["PRENV1"] add information about processing env, need to do this after the fact + # custom keywords + primary_hdr["DATATYPE"] = ("event_list", "Description of the data") + primary_hdr["ORIGAPID"] = (0xA0, "APID(s) of the originating data") + primary_hdr["ORIGFILE"] = (file_path.name, "Originating file(s)") # add common fits keywords fits_meta = read_fits_keyword_file( diff --git a/padre_meddea/io/file_tools.py b/padre_meddea/io/file_tools.py index 1928fa7..43007ff 100644 --- a/padre_meddea/io/file_tools.py +++ b/padre_meddea/io/file_tools.py @@ -18,19 +18,13 @@ count_packets, split_by_apid, ) +import astropy.io.fits as fits import padre_meddea -from padre_meddea.util.util import has_baseline +from padre_meddea import EPOCH, APID +from padre_meddea.util.util import has_baseline, calc_time, channel_to_pixel -__all__ = ["read_file", "read_raw_file"] - -APID = { - "spectrum": 0xA2, # decimal 162 - "photon": 0xA0, # decimal 160 - "housekeeping": 0xA3, # decimal 163 - "cmd_resp": 0xA5, # decimal 165 -} -EPOCH = dt.datetime(2000, 1, 1) +__all__ = ["read_file", "read_raw_file", "read_fits"] def read_file(filename: Path): @@ -50,9 +44,9 @@ def read_file(filename: Path): -------- """ if filename.suffix == "bin": # raw binary file - result = read_raw_file(filename) + result = read_raw_file(Path(filename)) elif filename.suffix == "fits": # level 0 or above - pass + read_fits(Path(filename)) else: raise ValueError("File extension {filename.suffix} not recognized.") return result @@ -63,7 +57,7 @@ def read_raw_file(filename: Path): Read a level 0 data file. Parameters - ---------- + --------- filename : Path A file to read @@ -83,6 +77,38 @@ def read_raw_file(filename: Path): return result +def read_fits(filename: Path): + """ + Read a fits file. + """ + hdu = fits.open(filename) + if (hdu[0].header["LEVEL"] == 0) and (hdu[0].header["DATATYPE"] == "event_list"): + num_events = len(hdu["SCI"].data["seqcount"]) + ph_times = calc_time( + hdu["sci"].data["pkttimes"], + hdu["sci"].data["pktclock"], + hdu["sci"].data["clocks"], + ) + # TODO: protect in case of non-pixel channel + pixel = np.array( + [channel_to_pixel(this_chan) for this_chan in hdu["sci"].data["channel"]], + dtype=np.uint8, + ) + event_list = TimeSeries( + time=ph_times, + data={ + "atod": hdu["sci"].data["atod"], + "asic": hdu["sci"].data["asic"], + "channel": hdu["sci"].data["channel"], + "pixel": pixel, + "clocks": hdu["sci"].data["clocks"], + "pktnum": hdu["sci"].data["seqcount"], + }, + ) + event_list.sort() + return event_list + + def parse_ph_packets(filename: Path): """Given a binary file, read only the photon packets and return an event list. @@ -131,25 +157,27 @@ def parse_ph_packets(filename: Path): packet_definition = packet_definition_ph() pkt = ccsdspy.VariableLength(packet_definition) ph_data = pkt.load(packet_stream, include_primary_header=True) - packet_time = ph_data["TIME_S"] + ph_data["TIME_CLOCKS"] / 20.0e6 + # packet_time = ph_data["TIME_S"] + ph_data["TIME_CLOCKS"] / 20.0e6 if has_baseline(filename): WORDS_PER_HIT = 4 else: WORDS_PER_HIT = 3 pkt_list = Table() - pkt_list['seqcount'] = ph_data["CCSDS_SEQUENCE_COUNT"] - pkt_list['time_s'] = ph_data["TIME_S"] - pkt_list['clocks'] = ph_data["TIME_CLOCKS"] - pkt_list['livetime'] = ph_data["LIVE_TIME"] - pkt_list['inttime'] = ph_data["INTEGRATION_TIME"] + pkt_list["seqcount"] = ph_data["CCSDS_SEQUENCE_COUNT"] + pkt_list["time_s"] = ph_data["TIME_S"] + pkt_list["clocks"] = ph_data["TIME_CLOCKS"] + pkt_list["livetime"] = ph_data["LIVE_TIME"] + pkt_list["inttime"] = ph_data["INTEGRATION_TIME"] total_hits = 0 for this_packet in ph_data["PIXEL_DATA"]: total_hits += len(this_packet[0::WORDS_PER_HIT]) if (total_hits - np.floor(total_hits)) != 0: - raise ValueError(f"Got non-integer number of hits {total_hits - np.floor(total_hits)}.") + raise ValueError( + f"Got non-integer number of hits {total_hits - np.floor(total_hits)}." + ) # hit_list definition # 0, raw id number @@ -158,27 +186,31 @@ def parse_ph_packets(filename: Path): # 3, energy 12 bits # 4, packet sequence number 12 bits # 5, clock number - # 6, baseline (if exists) + # 6, baseline (if exists) otherwise 0 - hit_list = np.zeros((7, total_hits), dtype="uint16") - time_stamps = np.zeros(total_hits) + hit_list = np.zeros((9, total_hits), dtype="uint16") + time_s = np.zeros(total_hits, dtype="uint32") + time_clk = np.zeros(total_hits, dtype="uint32") # 0 time, 1 asic_num, 2 channel num, 3 hit channel i = 0 - for this_pkt_num, this_ph_data, this_time in zip( + # iterate over packets + for this_pkt_num, this_ph_data, pkt_s, pkt_clk in zip( ph_data["CCSDS_SEQUENCE_COUNT"], ph_data["PIXEL_DATA"], - packet_time, + ph_data["TIME_S"], + ph_data["TIME_CLOCKS"], ): num_hits = len(this_ph_data[0::3]) ids = this_ph_data[1::3] asic_num = (ids & 0b11100000) >> 5 channel_num = ids & 0b00011111 hit_list[0, i : i + num_hits] = this_ph_data[0::WORDS_PER_HIT] - time_stamps[i : i + num_hits] = this_ph_data[0::WORDS_PER_HIT] * 12.8e-6 + this_time hit_list[1, i : i + num_hits] = asic_num hit_list[2, i : i + num_hits] = channel_num hit_list[4, i : i + num_hits] = this_pkt_num hit_list[5, i : i + num_hits] = this_ph_data[0::WORDS_PER_HIT] + time_s[i : i + num_hits] = pkt_s + time_clk[i : i + num_hits] = pkt_clk if WORDS_PER_HIT == 4: hit_list[6, i : i + num_hits] = this_ph_data[2::WORDS_PER_HIT] # baseline hit_list[3, i : i + num_hits] = this_ph_data[3::WORDS_PER_HIT] # hit energy @@ -188,17 +220,21 @@ def parse_ph_packets(filename: Path): event_list = Table() event_list["seqcount"] = hit_list[4, :] - event_list["clock"] = hit_list[5, :] + event_list["clocks"] = hit_list[5, :] event_list["asic"] = hit_list[1, :].astype(np.uint8) event_list["channel"] = hit_list[2, :].astype(np.uint8) event_list["atod"] = hit_list[3, :] event_list["baseline"] = hit_list[6, :] # if baseline not present then all zeros - - ph_times = Time([EPOCH + dt.timedelta(seconds=this_t) for this_t in time_stamps]) - event_list.meta.update({'DATE-BEG': ph_times.min().fits}) - event_list.meta.update({'DATE-END': ph_times.max().fits}) - event_list.meta.update({'DATE-AVG': ph_times[int(len(ph_times)/2.)].fits}) - #event_list = TimeSeries( + event_list["pkttimes"] = time_s + event_list["pktclock"] = time_clk + date_beg = calc_time(time_s[0], time_clk[0]) + date_end = calc_time(time_s[-1], time_clk[-1]) + event_list.meta.update({"DATE-BEG": date_beg.fits}) + event_list.meta.update({"DATE-END": date_end.fits}) + center_index = int(len(time_s) / 2.0) + date_avg = calc_time(time_s[center_index], time_clk[center_index]) + event_list.meta.update({"DATE-AVG": date_avg.fits}) + # event_list = TimeSeries( # time=ph_times, # data={ # "atod": hit_list[3, :], @@ -208,8 +244,8 @@ def parse_ph_packets(filename: Path): # "pktnum": hit_list[4, :], # # "num": np.ones(len(hit_list[2, :])), # } - #) - #event_list.sort() + # ) + # event_list.sort() # int_time = (event_list.time.max() - event_list.time.min()).to("min") # event_list.meta.update({"int_time": int_time}) # event_list.meta.update({"avg rate": (total_hits * u.ct) / int_time}) diff --git a/padre_meddea/util/util.py b/padre_meddea/util/util.py index 8bb407e..6036b95 100644 --- a/padre_meddea/util/util.py +++ b/padre_meddea/util/util.py @@ -6,16 +6,17 @@ from datetime import datetime, timezone import time from pathlib import Path +import warnings import numpy as np -from astropy.time import Time + +from astropy.time import Time, TimeDelta +import astropy.units as u from ccsdspy.utils import split_packet_bytes, split_by_apid +from padre_meddea import EPOCH -__all__ = [ - "create_science_filename", - "has_baseline" -] +__all__ = ["create_science_filename", "has_baseline"] TIME_FORMAT = "%Y%m%dT%H%M%S" VALID_DATA_LEVELS = ["l0", "l1", "ql", "l2", "l3", "l4"] @@ -105,7 +106,48 @@ def create_science_filename( return filename + FILENAME_EXTENSION -def has_baseline(filename: Path, packet_count=10): +def calc_time(pkt_time_s, pkt_time_clk, ph_clk=0): + """ + Convert times to a Time object + """ + deltat = TimeDelta( + pkt_time_s * u.s + pkt_time_clk * 0.05 * u.us + ph_clk * 12.8 * u.us + ) + result = Time(EPOCH + deltat) + + return result + + +def channel_to_pixel(channel: int) -> int: + """ + Given a channel pixel number, return the pixel number. + """ + CHANNEL_TO_PIX = { + 26: 0, + 15: 1, + 8: 2, + 1: 3, + 29: 4, + 13: 5, + 5: 6, + 0: 7, + 30: 8, + 21: 9, + 11: 10, + 3: 11, + 31: 12, + } + + if channel in CHANNEL_TO_PIX.keys(): + return CHANNEL_TO_PIX[channel] + else: + warnings.warn( + f"Found unconnected channel, {channel}. Returning channel + 12 ={channel+12}." + ) + return channel + 12 + + +def has_baseline(filename: Path, packet_count=10) -> bool: """Given a stream of photon packets, check whether the baseline measurement is included. Baseline packets have one extra word per photon for a total of 4 words (8 bytes). @@ -116,7 +158,7 @@ def has_baseline(filename: Path, packet_count=10): ---------- packet_bytes : byte string Photon packet bytes, must be an integer number of whole packets and greaterh - + Returns ------- result : bool @@ -132,4 +174,4 @@ def has_baseline(filename: Path, packet_count=10): num_hits = np.zeros(packet_count) for i in range(packet_count): num_hits[i] = (len(packet_bytes[i]) - HEADER_BYTES) / BYTES_PER_PHOTON - return np.sum(num_hits - np.floor(num_hits)) == 1 \ No newline at end of file + return np.sum(num_hits - np.floor(num_hits)) == 1 diff --git a/pyproject.toml b/pyproject.toml index 7d934b9..4e40be7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,7 +33,7 @@ dependencies = [ 'ccsdspy==1.3.*', 'matplotlib==3.9.*', 'swxsoc @ git+https://github.com/swxsoc/swxsoc.git@main', - 'gitpython==3.1.* + 'gitpython==3.1.*', ] [project.optional-dependencies] From b7b91a00ca9495406a94bfc850383da6b9d5fa87 Mon Sep 17 00:00:00 2001 From: Steven Daniel Christe Date: Fri, 18 Oct 2024 14:55:22 -0400 Subject: [PATCH 06/21] added and re-organized documentation --- docs/user-guide/index.rst | 1 + docs/user-guide/level0.rst | 93 ++++++++++++++++++++++------------- docs/user-guide/raw.rst | 41 +++++++++++++++ padre_meddea/io/file_tools.py | 5 +- 4 files changed, 105 insertions(+), 35 deletions(-) create mode 100644 docs/user-guide/raw.rst diff --git a/docs/user-guide/index.rst b/docs/user-guide/index.rst index 9ad7126..a62d94f 100644 --- a/docs/user-guide/index.rst +++ b/docs/user-guide/index.rst @@ -12,6 +12,7 @@ For more details checkout the :ref:`reference`. Brief Tour data + raw level0 spectrum customization diff --git a/docs/user-guide/level0.rst b/docs/user-guide/level0.rst index 8ac1ddc..4023f2d 100644 --- a/docs/user-guide/level0.rst +++ b/docs/user-guide/level0.rst @@ -6,36 +6,63 @@ Level 0 Data Overview ======== -This section document the format of the level 0 binary files. -Level 0 binary files are raw telemetry or command response packets generated by the instrument. - - -Photon Packet -------------- -The packet which provides individual hit data or photons for each detector. - -+-------------+-----------+--------+------------------------------------------------------------+ -| Name | Bit Size | Type | Description | -+=============+===========+========+============================================================+ -| HEADER WORD | 16 | UINT | value is always 65131 | -+-------------+-----------+--------+------------------------------------------------------------+ -| FLAGS | 16 | UNIT | various flags, To be defined | -+-------------+-----------+--------+------------------------------------------------------------+ -| Packet Size | 16 | UINT | The size of the packet which can be used to determine the | -| | | | number of hits included in the packet | -+-------------+-----------+--------+------------------------------------------------------------+ -| Time Stamp | 48 | UINT | | -+-------------+-----------+--------+------------------------------------------------------------+ -| Checksum | 16 | UINT | For data integrity | -+-------------+-----------+--------+------------------------------------------------------------+ -| Pixel data | | | This field is repeated based on the number of hits detected| -+-------------+-----------+--------+------------------------------------------------------------+ - - -+-------------+-----------+--------+------------------------------------------------------------+ -| Pixel Data | | | | -+-------------+-----------+--------+------------------------------------------------------------+ -| Detector id | 16 | UINT | The detector id for the location of the hit | -+-------------+-----------+--------+------------------------------------------------------------+ -| Hit energy | 16 | UINT | The ADC value for the energy of the hit | -+-------------+-----------+--------+------------------------------------------------------------+ +Level 0 data are provided in the `FITS file `__ format. +For more information on how to read or write FITS file see `astropy.fits `__. +This section describes the organization the level 0 FITS files. +Level 0 fits files generally include the unconverted data from the raw binary files of ccsds packets. +The purpose of these files is to provide the raw data from the raw binary files in a more convinient form for analysis. +It also provides metadata information which summary the data in the file. + +Level 0 event files +------------------- + +This file contains the data from all events that triggered the detectors. +They consist of 3 HDUs including the primary HDU. +The primary HDU contains no data and is only used for metadata. +The two other HDUs are named `SCI` and `PKT`. +`SCI` includes the event data while `PKT` includes the packet data. +Each data packet may include one of more event therefore there is a one to many relationship between them. +In order to understand the relationship between the events and packets, each event provides the associated packet sequence number. +This sequence number can be used to lookup the packet data for that event. + +Primary HDU +*********** +No data is provided. +Stay tuned for a list of metadata + +PKT HDU +******* +The following columns are provided for each data packet. +The bits column provide the number of significant bits and not the bit length of the column itself. +The columns in the FITS file are provided in the smallest possible data type. + +======== ============================================= ==== +name bits +======== ============================================= ==== +seqcount packet sequence number, should be consecutive 12 +pkttimes the packet time in seconds since EPOCH 32 +pktclock the packet subsecond time in clocks 32 +livetime live time 16 +inttime integration time in real time 16 +flags flags 16 +======== ============================================= ==== + +SCI HDU +******* +The following columns are provided for each event or photon detected. +The bits column provide the number of significant bits and not the bit length of the column itself. +The columns in the FITS file are provided in the smallest possible data type. + +======== ============================================================================================ ==== +name description bits +======== ============================================================================================ ==== +seqcount packet sequence number 12 +clocks the clock number 16 +asic the asic number or detector id 3 +channel the asic channel which is related to the pixel 5 +atod the uncalibrated energy of the event in ADC counts 12 +baseline the baseline measurement if exists, otherwise all zeros 12 +pkttimes the packet time in seconds since EPOCH, also exists in PKT, 32 +pktclock the packet time in clocks since EPOCH, also exists in PKT 32 +======== ============================================================================================ ==== + diff --git a/docs/user-guide/raw.rst b/docs/user-guide/raw.rst new file mode 100644 index 0000000..e4a9823 --- /dev/null +++ b/docs/user-guide/raw.rst @@ -0,0 +1,41 @@ +.. _raw: + +*************** +Raw Binary Data +*************** + +Overview +======== +This section document the format of the raw binary files. +Level 0 binary files are raw telemetry or command response packets generated by the instrument. + + +Photon Packet +------------- +The packet which provides individual hit data or photons for each detector. + ++-------------+-----------+--------+------------------------------------------------------------+ +| Name | Bit Size | Type | Description | ++=============+===========+========+============================================================+ +| HEADER WORD | 16 | UINT | value is always 65131 | ++-------------+-----------+--------+------------------------------------------------------------+ +| FLAGS | 16 | UNIT | various flags, To be defined | ++-------------+-----------+--------+------------------------------------------------------------+ +| Packet Size | 16 | UINT | The size of the packet which can be used to determine the | +| | | | number of hits included in the packet | ++-------------+-----------+--------+------------------------------------------------------------+ +| Time Stamp | 48 | UINT | | ++-------------+-----------+--------+------------------------------------------------------------+ +| Checksum | 16 | UINT | For data integrity | ++-------------+-----------+--------+------------------------------------------------------------+ +| Pixel data | | | This field is repeated based on the number of hits detected| ++-------------+-----------+--------+------------------------------------------------------------+ + + ++-------------+-----------+--------+------------------------------------------------------------+ +| Pixel Data | | | | ++-------------+-----------+--------+------------------------------------------------------------+ +| Detector id | 16 | UINT | The detector id for the location of the hit | ++-------------+-----------+--------+------------------------------------------------------------+ +| Hit energy | 16 | UINT | The ADC value for the energy of the hit | ++-------------+-----------+--------+------------------------------------------------------------+ diff --git a/padre_meddea/io/file_tools.py b/padre_meddea/io/file_tools.py index 43007ff..2502524 100644 --- a/padre_meddea/io/file_tools.py +++ b/padre_meddea/io/file_tools.py @@ -165,10 +165,11 @@ def parse_ph_packets(filename: Path): pkt_list = Table() pkt_list["seqcount"] = ph_data["CCSDS_SEQUENCE_COUNT"] - pkt_list["time_s"] = ph_data["TIME_S"] - pkt_list["clocks"] = ph_data["TIME_CLOCKS"] + pkt_list["pkttimes"] = ph_data["TIME_S"] + pkt_list["pktclock"] = ph_data["TIME_CLOCKS"] pkt_list["livetime"] = ph_data["LIVE_TIME"] pkt_list["inttime"] = ph_data["INTEGRATION_TIME"] + pkt_list["flags"] = ph_data["FLAGS"] total_hits = 0 for this_packet in ph_data["PIXEL_DATA"]: From 1628df4d8390f54b617df5ebea54cf5da96f4e92 Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Sat, 19 Oct 2024 17:24:16 -0400 Subject: [PATCH 07/21] added housekeeping l1 fits files --- docs/user-guide/level0.rst | 8 ++ padre_meddea/calibration/calibration.py | 131 +++++++++++++++++------- padre_meddea/io/file_tools.py | 62 ++++++----- padre_meddea/util/util.py | 2 +- 4 files changed, 137 insertions(+), 66 deletions(-) diff --git a/docs/user-guide/level0.rst b/docs/user-guide/level0.rst index 4023f2d..107ad50 100644 --- a/docs/user-guide/level0.rst +++ b/docs/user-guide/level0.rst @@ -66,3 +66,11 @@ pkttimes the packet time in seconds since EPOCH, also exists in PKT, pktclock the packet time in clocks since EPOCH, also exists in PKT 32 ======== ============================================================================================ ==== +Level 0 spectrum files +---------------------- +Summary spectra are created for 24 pixels at a regular cadence (normally every 10 s) +Each spectrum has a total of 512 energy bins. + +Level 0 housekeeping files +-------------------------- +These files contain housekeeping data as described in the housekeeping packet. \ No newline at end of file diff --git a/padre_meddea/calibration/calibration.py b/padre_meddea/calibration/calibration.py index d74ba7b..2285f86 100644 --- a/padre_meddea/calibration/calibration.py +++ b/padre_meddea/calibration/calibration.py @@ -17,7 +17,7 @@ from padre_meddea import log from padre_meddea.io import file_tools -from padre_meddea.util.util import create_science_filename +from padre_meddea.util.util import create_science_filename, calc_time from padre_meddea.io.file_tools import read_raw_file __all__ = [ @@ -70,44 +70,14 @@ def process_file(filename: Path, overwrite=False) -> list: this_str, ) - primary_hdr["DATEREF"] = (primary_hdr["DATE-BEG"], "Reference date") primary_hdr["LEVEL"] = (0, "Data level of fits file") # add processing information - primary_hdr["PRSTEP1"] = ("PROCESS Raw to L1", "Processing step type") - primary_hdr["PRPROC1"] = ( - "padre_meddea.calibration.process", - "Name of procedure performing PRSTEP1", - ) - primary_hdr["PRPVER1"] = ( - padre_meddea.__version__, - "Version of procedure PRPROC1", - ) - primary_hdr["PRLIB1A"] = ( - "padre_meddea", - "Software library containing PRPROC1", - ) - primary_hdr["PRVER1A"] = (padre_meddea.__version__, "Version of PRLIB1A") - repo = git.Repo(search_parent_directories=True) - primary_hdr["PRHSH1A"] = ( - repo.head.object.hexsha, - "GIT commit hash for PRLIB1A", - ) - primary_hdr["PRBRA1A"] = ( - repo.active_branch.name, - "GIT/SVN repository branch of PRLIB1A", - ) - commits = list(repo.iter_commits("main", max_count=1)) - primary_hdr["PRVER1B"] = ( - Time(commits[0].committed_datetime).fits, - "Date of last commit of PRLIB1B", - ) - # primary_hdr["PRLOG1"] add log information, need to do this after the fact - # primary_hdr["PRENV1"] add information about processing env, need to do this after the fact - + primary_hdr = add_process_info_to_header(primary_hdr) + # custom keywords primary_hdr["DATATYPE"] = ("event_list", "Description of the data") - primary_hdr["ORIGAPID"] = (0xA0, "APID(s) of the originating data") + primary_hdr["ORIGAPID"] = (padre_meddea.APID["photon"], "APID(s) of the originating data") primary_hdr["ORIGFILE"] = (file_path.name, "Originating file(s)") # add common fits keywords @@ -142,15 +112,54 @@ def process_file(filename: Path, overwrite=False) -> list: hdul.writeto(path, overwrite=overwrite) # Store the output file path in a list - output_files = [path] + output_files.append(path) if parsed_data["housekeeping"] is not None: hk_data = parsed_data["housekeeping"] - hk_data.meta["INSTRUME"] = "meddea" + # send data to AWS Timestream for Grafana dashboard + record_timeseries(hk_data, "housekeeping") + hk_table = Table(hk_data) + primary_hdr = fits.Header() + # fill in metadata + primary_hdr["DATE"] = (Time.now().fits, "FITS file creation date in UTC") + primary_hdr["LEVEL"] = (0, "Data level of fits file") + primary_hdr["DATATYPE"] = ("housekeeping", "Description of the data") + primary_hdr["ORIGAPID"] = (padre_meddea.APID["housekeeping"], "APID(s) of the originating data") + primary_hdr["ORIGFILE"] = (file_path.name, "Originating file(s)") + date_beg = calc_time(hk_data['timestamp'][0]) + primary_hdr["DATEREF"] = (date_beg.fits, "Reference date") - if "CHECKSUM" in hk_data.colnames: - hk_data.remove_column("CHECKSUM") + # add processing information + primary_hdr = add_process_info_to_header(primary_hdr) + + # add common fits keywords + fits_meta = read_fits_keyword_file( + padre_meddea._data_directory / "fits_keywords_primaryhdu.csv" + ) + for row in fits_meta: + primary_hdr[row["keyword"]] = (row["value"], row["comment"]) + hk_table['seqcount'] = hk_table["CCSDS_SEQUENCE_COUNT"] + colnames_to_remove = ["CCSDS_VERSION_NUMBER", "CCSDS_PACKET_TYPE", "CCSDS_SECONDARY_FLAG", "CCSDS_SEQUENCE_FLAG", "CCSDS_APID", "CCSDS_SEQUENCE_COUNT", "CCSDS_PACKET_LENGTH", "CHECKSUM", "time"] + for this_col in colnames_to_remove: + if this_col in hk_table.colnames: + hk_table.remove_column(this_col) + + empty_primary = fits.PrimaryHDU(header=primary_hdr) + hk_hdu = fits.BinTableHDU(hk_table, name="HK") + hk_hdu.add_checksum() + hdul = fits.HDUList([empty_primary, hk_hdu]) + + path = create_science_filename( + time=date_beg, + level="l1", + descriptor="hk", + test=True, + version="0.1.0", + ) + hdul.writeto(path, overwrite=overwrite) + output_files.append(path) + + - record_timeseries(hk_data, "housekeeping") # calibrated_file = calibrate_file(data_filename) # data_plot_files = plot_file(data_filename) @@ -167,6 +176,50 @@ def raw_to_l0(filename: Path): data = file_tools.read_raw_file(filename) +def add_process_info_to_header(header: fits.Header) -> fits.Header: + """Add processing info metadata to fits header. + + Parameters + ---------- + header : fits.Header + + Returns + ------- + header : fits.Header + """ + header["PRSTEP1"] = ("PROCESS Raw to L1", "Processing step type") + header["PRPROC1"] = ( + "padre_meddea.calibration.process", + "Name of procedure performing PRSTEP1", + ) + header["PRPVER1"] = ( + padre_meddea.__version__, + "Version of procedure PRPROC1", + ) + header["PRLIB1A"] = ( + "padre_meddea", + "Software library containing PRPROC1", + ) + header["PRVER1A"] = (padre_meddea.__version__, "Version of PRLIB1A") + repo = git.Repo(search_parent_directories=True) + header["PRHSH1A"] = ( + repo.head.object.hexsha, + "GIT commit hash for PRLIB1A", + ) + header["PRBRA1A"] = ( + repo.active_branch.name, + "GIT/SVN repository branch of PRLIB1A", + ) + commits = list(repo.iter_commits("main", max_count=1)) + header["PRVER1B"] = ( + Time(commits[0].committed_datetime).fits, + "Date of last commit of PRLIB1B", + ) + # primary_hdr["PRLOG1"] add log information, need to do this after the fact + # primary_hdr["PRENV1"] add information about processing env, need to do this after the fact + return header + + def get_calibration_file(time: Time) -> Path: """ Given a time, return the appropriate calibration file. diff --git a/padre_meddea/io/file_tools.py b/padre_meddea/io/file_tools.py index 2502524..1d3480f 100644 --- a/padre_meddea/io/file_tools.py +++ b/padre_meddea/io/file_tools.py @@ -82,31 +82,41 @@ def read_fits(filename: Path): Read a fits file. """ hdu = fits.open(filename) + if (hdu[0].header["LEVEL"] == 0) and (hdu[0].header["DATATYPE"] == "event_list"): - num_events = len(hdu["SCI"].data["seqcount"]) - ph_times = calc_time( - hdu["sci"].data["pkttimes"], - hdu["sci"].data["pktclock"], - hdu["sci"].data["clocks"], - ) - # TODO: protect in case of non-pixel channel - pixel = np.array( - [channel_to_pixel(this_chan) for this_chan in hdu["sci"].data["channel"]], - dtype=np.uint8, - ) - event_list = TimeSeries( - time=ph_times, - data={ - "atod": hdu["sci"].data["atod"], - "asic": hdu["sci"].data["asic"], - "channel": hdu["sci"].data["channel"], - "pixel": pixel, - "clocks": hdu["sci"].data["clocks"], - "pktnum": hdu["sci"].data["seqcount"], - }, - ) - event_list.sort() + event_list = read_fits_l0_event_list(filename) # do I need to close the file since it is being opened again right after this? return event_list + else: + raise ValueError(F"File contents of {filename} not recogized.") + + +def read_fits_l0_event_list(filename: Path) -> TimeSeries: + """ + """ + hdu = fits.open(filename) + num_events = len(hdu["SCI"].data["seqcount"]) + ph_times = calc_time( + hdu["sci"].data["pkttimes"], + hdu["sci"].data["pktclock"], + hdu["sci"].data["clocks"], + ) + pixel = np.array( + [channel_to_pixel(this_chan) for this_chan in hdu["sci"].data["channel"]], + dtype=np.uint8, + ) + event_list = TimeSeries( + time=ph_times, + data={ + "atod": hdu["sci"].data["atod"], + "asic": hdu["sci"].data["asic"], + "channel": hdu["sci"].data["channel"], + "pixel": pixel, + "clocks": hdu["sci"].data["clocks"], + "pktnum": hdu["sci"].data["seqcount"], + }, + ) + event_list.sort() + return event_list def parse_ph_packets(filename: Path): @@ -273,9 +283,9 @@ def parse_hk_packets(filename: Path): return None packet_definition = packet_definition_hk() pkt = ccsdspy.FixedLength(packet_definition) - hk_data = pkt.load(packet_bytes) + hk_data = pkt.load(packet_bytes, include_primary_header=True) hk_timestamps = [ - dt.timedelta(seconds=int(this_t)) + EPOCH for this_t in hk_data["TIMESTAMP"] + dt.timedelta(seconds=int(this_t)) + EPOCH for this_t in hk_data["timestamp"] ] hk_data = TimeSeries(time=hk_timestamps, data=hk_data) return hk_data @@ -336,7 +346,7 @@ def packet_definition_hk(): """Return the packet definiton for the housekeeping packets.""" hk_table = ascii.read(padre_meddea._data_directory / "hk_packet_def.csv") hk_table.add_index("name") - p = [PacketField(name="TIMESTAMP", data_type="uint", bit_length=32)] + p = [PacketField(name="timestamp", data_type="uint", bit_length=32)] for this_hk in hk_table["name"]: p += [PacketField(name=this_hk, data_type="uint", bit_length=16)] p += [PacketField(name="CHECKSUM", data_type="uint", bit_length=16)] diff --git a/padre_meddea/util/util.py b/padre_meddea/util/util.py index 6036b95..b248fa4 100644 --- a/padre_meddea/util/util.py +++ b/padre_meddea/util/util.py @@ -106,7 +106,7 @@ def create_science_filename( return filename + FILENAME_EXTENSION -def calc_time(pkt_time_s, pkt_time_clk, ph_clk=0): +def calc_time(pkt_time_s, pkt_time_clk=0, ph_clk=0): """ Convert times to a Time object """ From a826db303286ae6a7340c60af5888abee6344c3a Mon Sep 17 00:00:00 2001 From: Steven Daniel Christe Date: Wed, 23 Oct 2024 14:28:27 -0400 Subject: [PATCH 08/21] cleaning up fits tools --- docs/user-guide/level0.rst | 3 +- padre_meddea/calibration/calibration.py | 63 +------------ padre_meddea/data/README.rst | 2 + padre_meddea/data/calibration/README.rst | 2 +- padre_meddea/data/fits/README.rst | 13 +++ .../data/fits_keywords_primaryhdu.csv | 11 --- padre_meddea/io/fits_tools.py | 92 +++++++++++++++++++ padre_meddea/tests/test_fits_tools.py | 41 +++++++++ 8 files changed, 152 insertions(+), 75 deletions(-) create mode 100644 padre_meddea/data/fits/README.rst delete mode 100644 padre_meddea/data/fits_keywords_primaryhdu.csv create mode 100644 padre_meddea/io/fits_tools.py create mode 100644 padre_meddea/tests/test_fits_tools.py diff --git a/docs/user-guide/level0.rst b/docs/user-guide/level0.rst index 107ad50..c4ec64d 100644 --- a/docs/user-guide/level0.rst +++ b/docs/user-guide/level0.rst @@ -73,4 +73,5 @@ Each spectrum has a total of 512 energy bins. Level 0 housekeeping files -------------------------- -These files contain housekeeping data as described in the housekeeping packet. \ No newline at end of file +These files contain housekeeping data as described in the housekeeping packet. +It also includes any register read responses that may exist during that time period. diff --git a/padre_meddea/calibration/calibration.py b/padre_meddea/calibration/calibration.py index 2285f86..229af5c 100644 --- a/padre_meddea/calibration/calibration.py +++ b/padre_meddea/calibration/calibration.py @@ -11,11 +11,10 @@ from astropy.table import Table from swxsoc.util.util import record_timeseries -import git import padre_meddea from padre_meddea import log -from padre_meddea.io import file_tools +from padre_meddea.io import file_tools, fits_tools from padre_meddea.util.util import create_science_filename, calc_time from padre_meddea.io.file_tools import read_raw_file @@ -80,13 +79,6 @@ def process_file(filename: Path, overwrite=False) -> list: primary_hdr["ORIGAPID"] = (padre_meddea.APID["photon"], "APID(s) of the originating data") primary_hdr["ORIGFILE"] = (file_path.name, "Originating file(s)") - # add common fits keywords - fits_meta = read_fits_keyword_file( - padre_meddea._data_directory / "fits_keywords_primaryhdu.csv" - ) - for row in fits_meta: - primary_hdr[row["keyword"]] = (row["value"], row["comment"]) - empty_primary = fits.PrimaryHDU(header=primary_hdr) pkt_hdu = fits.BinTableHDU(pkt_list, name="PKT") pkt_hdu.add_checksum() @@ -176,50 +168,6 @@ def raw_to_l0(filename: Path): data = file_tools.read_raw_file(filename) -def add_process_info_to_header(header: fits.Header) -> fits.Header: - """Add processing info metadata to fits header. - - Parameters - ---------- - header : fits.Header - - Returns - ------- - header : fits.Header - """ - header["PRSTEP1"] = ("PROCESS Raw to L1", "Processing step type") - header["PRPROC1"] = ( - "padre_meddea.calibration.process", - "Name of procedure performing PRSTEP1", - ) - header["PRPVER1"] = ( - padre_meddea.__version__, - "Version of procedure PRPROC1", - ) - header["PRLIB1A"] = ( - "padre_meddea", - "Software library containing PRPROC1", - ) - header["PRVER1A"] = (padre_meddea.__version__, "Version of PRLIB1A") - repo = git.Repo(search_parent_directories=True) - header["PRHSH1A"] = ( - repo.head.object.hexsha, - "GIT commit hash for PRLIB1A", - ) - header["PRBRA1A"] = ( - repo.active_branch.name, - "GIT/SVN repository branch of PRLIB1A", - ) - commits = list(repo.iter_commits("main", max_count=1)) - header["PRVER1B"] = ( - Time(commits[0].committed_datetime).fits, - "Date of last commit of PRLIB1B", - ) - # primary_hdr["PRLOG1"] add log information, need to do this after the fact - # primary_hdr["PRENV1"] add information about processing env, need to do this after the fact - return header - - def get_calibration_file(time: Time) -> Path: """ Given a time, return the appropriate calibration file. @@ -262,12 +210,3 @@ def read_calibration_file(calib_filename: Path): # if can't read the file return None - - -def read_fits_keyword_file(csv_file: Path): - """Read csv file with default fits metadata information.""" - fits_meta_table = ascii.read( - padre_meddea._data_directory / "fits_keywords_primaryhdu.csv", - format="csv", - ) - return fits_meta_table diff --git a/padre_meddea/data/README.rst b/padre_meddea/data/README.rst index ed399b9..07fe7e6 100644 --- a/padre_meddea/data/README.rst +++ b/padre_meddea/data/README.rst @@ -39,3 +39,5 @@ Stores detector constants. hk_channel_defs.csv ------------------- Stores the definitions for the values provided in housekeeping packets. + +fits_ \ No newline at end of file diff --git a/padre_meddea/data/calibration/README.rst b/padre_meddea/data/calibration/README.rst index 739027f..d54ee28 100644 --- a/padre_meddea/data/calibration/README.rst +++ b/padre_meddea/data/calibration/README.rst @@ -1,4 +1,4 @@ -Calbiration directory +Calibration directory ===================== This directory contains calibration files included with the package source diff --git a/padre_meddea/data/fits/README.rst b/padre_meddea/data/fits/README.rst new file mode 100644 index 0000000..8fde905 --- /dev/null +++ b/padre_meddea/data/fits/README.rst @@ -0,0 +1,13 @@ +FITS file data directory +======================== + +This directory contains data related to FITS file creation. + +fits_keyword_dict.csv +--------------------- +Keyword names and comments. +Used to lookup standard comment based on keyword. + +fits_keyword_primaryhdu.csv +--------------------------- +Standard keyword values. \ No newline at end of file diff --git a/padre_meddea/data/fits_keywords_primaryhdu.csv b/padre_meddea/data/fits_keywords_primaryhdu.csv deleted file mode 100644 index 74223dc..0000000 --- a/padre_meddea/data/fits_keywords_primaryhdu.csv +++ /dev/null @@ -1,11 +0,0 @@ -keyword,value,comment -AUTHOR,Steven D. Christe,Who designed the observation -CREATOR,padre_meddea,Name of software pipeline that produced the FITS file -DETECTOR,meddea,Name of the detector -INFO_URL,https://padre-meddea.readthedocs.io/en/latest/user-guide/data.html,a human-readable web page describing the data -OBSRVTRY,PADRE,Name of the observatory -TIMESYS,UTC,Time scale of the time-related keywords -TELESCOP,PADRE/MeDDEA,Telescope/Sensor name -INSTRUME,MeDDEA,Instrument name -MISSION,PADRE,Mission name -ORIGIN,NASA GSFC,File originator \ No newline at end of file diff --git a/padre_meddea/io/fits_tools.py b/padre_meddea/io/fits_tools.py new file mode 100644 index 0000000..9dbb0f9 --- /dev/null +++ b/padre_meddea/io/fits_tools.py @@ -0,0 +1,92 @@ +""" +This module provides a utilities to manage fits files reading and writing. +""" +import re +import git + +from astropy.io import ascii +import astropy.io.fits as fits +from astropy.time import Time + +import padre_meddea + +FITS_HDR0 = ascii.read(padre_meddea._data_directory / "fits" / "fits_keywords_primaryhdu.csv",format="csv") +FITS_HDR0.add_index('keyword') +FITS_HDR_KEYTOCOMMENT = ascii.read(padre_meddea._data_directory / "fits" / "fits_keywords_dict.csv",format="csv") +FITS_HDR_KEYTOCOMMENT.add_index('keyword') + + +def get_primary_header() -> fits.Header: + """Return a standard FITS file primary header.""" + header = fits.Header() + for row in FITS_HDR0: + this_comment = FITS_HDR_KEYTOCOMMENT.loc[row["keyword"]]['comment'] + header[row["keyword"]] = (row["value"], this_comment) + + return header + + +def get_std_comment(keyword: str) -> str: + """Given a keyword, return the standard comment for a header card.""" + if keyword.upper() in FITS_HDR_KEYTOCOMMENT['keyword']: + return FITS_HDR_KEYTOCOMMENT.loc[keyword]['comment'] + for this_row in FITS_HDR_KEYTOCOMMENT: + res = re.fullmatch(this_row['keyword'], keyword) + if res: + comment = this_row['comment'] + if len(res.groupdict()) > 0: # check if there was a match + for key, value in res.groupdict().items(): + comment = comment.replace(f"<{key}>", value) + return comment + + +def add_process_info_to_header(header: fits.Header, n=1) -> fits.Header: + """Add processing info metadata to a fits header. + + It adds the following SOLARNET compatible FITS cards; + PRSTEPn, PRPROCn, PRPVERn, PRLIBnA, PRVERnA, PRLIBnA, PRHSHnA, PRVERnB + + Parameters + ---------- + header : fits.Header + The fits header to add the new cards to + n : int, default 1 + The processing step number. Must be greater than or equal to 1. + + Returns + ------- + header : fits.Header + """ + if n < 1: + ValueError("Processing number, n, must be greater than or equal to 1.") + header[f"PRSTEP{n}"] = ("PROCESS Raw to L1", get_std_comment(f'PRSTEP{n}')) + header[f"PRPROC{n}"] = ( + "padre_meddea.calibration.process", + get_std_comment(f'PRPROC{n}'), + ) + header[f"PRPVER{n}"] = ( + padre_meddea.__version__, + get_std_comment(f'PRPVER{n}'), + ) + header[f"PRLIB{n}A"] = ( + "padre_meddea", + get_std_comment(f'PRLIB{n}A'), + ) + header[f"PRVER{n}A"] = (padre_meddea.__version__, get_std_comment(f'PRVER{n}A')) + repo = git.Repo(search_parent_directories=True) + header[f"PRHSH{n}A"] = ( + repo.head.object.hexsha, + get_std_comment(f'PRHSH{n}A'), + ) + header[f"PRBRA{n}A"] = ( + repo.active_branch.name, + get_std_comment(f'PRBRA{n}A'), + ) + commits = list(repo.iter_commits("main", max_count=1)) + header[f"PRVER{n}B"] = ( + Time(commits[0].committed_datetime).fits, + get_std_comment(f'PRVER{n}B'), + ) + # primary_hdr["PRLOG1"] add log information, need to do this after the fact + # primary_hdr["PRENV1"] add information about processing env, need to do this after the fact + return header \ No newline at end of file diff --git a/padre_meddea/tests/test_fits_tools.py b/padre_meddea/tests/test_fits_tools.py new file mode 100644 index 0000000..7b988fd --- /dev/null +++ b/padre_meddea/tests/test_fits_tools.py @@ -0,0 +1,41 @@ +"""Test for the fits_tools module""" +import pytest + +import astropy.io.fits as fits + +from padre_meddea.io.fits_tools import * + + +def test_comment_lookup_hdr0(): + """Test that all keywords in fits_keyword_primaryhdu are listed in fits_keyword_dict""" + hdr0_keywords = list(FITS_HDR0['keyword']) + keyword_to_comment = list(FITS_HDR_KEYTOCOMMENT['keyword']) + for this_keyword in hdr0_keywords: + assert this_keyword in keyword_to_comment + + +def test_get_primary_header(): + assert isinstance(get_primary_header(), fits.Header) + + +def test_add_process_info_to_header(): + """Test that new header cards are added.""" + header = get_primary_header() + orig_header = header.copy() + header = add_process_info_to_header(header) + # check that keywords were added + assert len(header) > len(orig_header) + orig_keywords = [this_keyword for this_keyword in orig_header] + # check that the cards that were added have content + for this_card in header.cards: + if this_card.keyword not in orig_keywords: + assert len(this_card.value) > 0 + assert len(this_card.comment) > 0 + +@pytest.mark.parametrize("test_input,expected", [("PRSTEP2", "Processing step type"), + ("PRSTEP1", "Processing step type"), + ("PRPROC3", "Name of procedure performing PRSTEP3"), + ("PRHSH5A", "GIT commit hash for PRLIB5A") + ]) +def test_get_std_comment(test_input, expected): + assert get_std_comment(test_input) == expected From 51bbde9562fd4efded72e0540c95fab86bb8acb6 Mon Sep 17 00:00:00 2001 From: Steven Daniel Christe Date: Fri, 25 Oct 2024 15:53:22 -0400 Subject: [PATCH 09/21] first working commit of fits file creation --- padre_meddea/calibration/calibration.py | 132 ++++++++++-------- padre_meddea/data/fits/fits_keywords_dict.csv | 28 ++++ .../data/fits/fits_keywords_primaryhdu.csv | 11 ++ padre_meddea/io/file_tools.py | 132 +++++++++++++----- padre_meddea/io/fits_tools.py | 65 +++++---- padre_meddea/tests/test_fits_tools.py | 20 ++- padre_meddea/util/util.py | 3 +- 7 files changed, 261 insertions(+), 130 deletions(-) create mode 100644 padre_meddea/data/fits/fits_keywords_dict.csv create mode 100644 padre_meddea/data/fits/fits_keywords_primaryhdu.csv diff --git a/padre_meddea/calibration/calibration.py b/padre_meddea/calibration/calibration.py index 229af5c..b8d8de7 100644 --- a/padre_meddea/calibration/calibration.py +++ b/padre_meddea/calibration/calibration.py @@ -18,6 +18,11 @@ from padre_meddea.util.util import create_science_filename, calc_time from padre_meddea.io.file_tools import read_raw_file +from padre_meddea.io.fits_tools import ( + add_process_info_to_header, + get_primary_header, + get_std_comment, +) __all__ = [ "process_file", @@ -51,40 +56,28 @@ def process_file(filename: Path, overwrite=False) -> list: parsed_data = read_raw_file(file_path) if parsed_data["photons"] is not None: # we have event list data event_list, pkt_list = parsed_data["photons"] + primary_hdr = get_primary_header() + primary_hdr = add_process_info_to_header(primary_hdr) + primary_hdr["LEVEL"] = (0, get_std_comment("LEVEL")) + primary_hdr["DATATYPE"] = ("event_list", get_std_comment("DATATYPE")) + primary_hdr["ORIGAPID"] = ( + padre_meddea.APID["photon"], + get_std_comment("ORIGAPID"), + ) + primary_hdr["ORIGFILE"] = (file_path.name, get_std_comment("ORIGFILE")) - primary_hdr = fits.Header() - - # fill in metadata - primary_hdr["DATE"] = (Time.now().fits, "FITS file creation date in UTC") - for this_keyword, this_str in zip( - ["DATE-BEG", "DATE-END", "DATE-AVG"], - [ - "Acquisition start time", - "Acquisition end time", - "Average time of acquisition", - ], - ): + for this_keyword in ["DATE-BEG", "DATE-END", "DATE-AVG"]: primary_hdr[this_keyword] = ( event_list.meta.get(this_keyword, ""), - this_str, + get_std_comment(this_keyword), ) - primary_hdr["LEVEL"] = (0, "Data level of fits file") - - # add processing information - primary_hdr = add_process_info_to_header(primary_hdr) - - # custom keywords - primary_hdr["DATATYPE"] = ("event_list", "Description of the data") - primary_hdr["ORIGAPID"] = (padre_meddea.APID["photon"], "APID(s) of the originating data") - primary_hdr["ORIGFILE"] = (file_path.name, "Originating file(s)") - - empty_primary = fits.PrimaryHDU(header=primary_hdr) + empty_primary_hdu = fits.PrimaryHDU(header=primary_hdr) pkt_hdu = fits.BinTableHDU(pkt_list, name="PKT") pkt_hdu.add_checksum() hit_hdu = fits.BinTableHDU(event_list, name="SCI") hit_hdu.add_checksum() - hdul = fits.HDUList([empty_primary, hit_hdu, pkt_hdu]) + hdul = fits.HDUList([empty_primary_hdu, hit_hdu, pkt_hdu]) path = create_science_filename( time=primary_hdr["DATE-BEG"], @@ -102,7 +95,6 @@ def process_file(filename: Path, overwrite=False) -> list: # Write the file, with the overwrite option controlled by the environment variable hdul.writeto(path, overwrite=overwrite) - # Store the output file path in a list output_files.append(path) if parsed_data["housekeeping"] is not None: @@ -110,35 +102,70 @@ def process_file(filename: Path, overwrite=False) -> list: # send data to AWS Timestream for Grafana dashboard record_timeseries(hk_data, "housekeeping") hk_table = Table(hk_data) - primary_hdr = fits.Header() - # fill in metadata - primary_hdr["DATE"] = (Time.now().fits, "FITS file creation date in UTC") - primary_hdr["LEVEL"] = (0, "Data level of fits file") - primary_hdr["DATATYPE"] = ("housekeeping", "Description of the data") - primary_hdr["ORIGAPID"] = (padre_meddea.APID["housekeeping"], "APID(s) of the originating data") - primary_hdr["ORIGFILE"] = (file_path.name, "Originating file(s)") - date_beg = calc_time(hk_data['timestamp'][0]) - primary_hdr["DATEREF"] = (date_beg.fits, "Reference date") - - # add processing information - primary_hdr = add_process_info_to_header(primary_hdr) - # add common fits keywords - fits_meta = read_fits_keyword_file( - padre_meddea._data_directory / "fits_keywords_primaryhdu.csv" + primary_hdr = get_primary_header() + primary_hdr = add_process_info_to_header(primary_hdr) + primary_hdr["LEVEL"] = (0, get_std_comment("LEVEL")) + primary_hdr["DATATYPE"] = ("housekeeping", get_std_comment("DATATYPE")) + primary_hdr["ORIGAPID"] = ( + padre_meddea.APID["housekeeping"], + get_std_comment("ORIGAPID"), ) - for row in fits_meta: - primary_hdr[row["keyword"]] = (row["value"], row["comment"]) - hk_table['seqcount'] = hk_table["CCSDS_SEQUENCE_COUNT"] - colnames_to_remove = ["CCSDS_VERSION_NUMBER", "CCSDS_PACKET_TYPE", "CCSDS_SECONDARY_FLAG", "CCSDS_SEQUENCE_FLAG", "CCSDS_APID", "CCSDS_SEQUENCE_COUNT", "CCSDS_PACKET_LENGTH", "CHECKSUM", "time"] + primary_hdr["ORIGFILE"] = (file_path.name, get_std_comment("ORIGFILE")) + + date_beg = calc_time(hk_data["timestamp"][0]) + primary_hdr["DATEREF"] = (date_beg.fits, get_std_comment("DATEREF")) + + hk_table["seqcount"] = hk_table["CCSDS_SEQUENCE_COUNT"] + colnames_to_remove = [ + "CCSDS_VERSION_NUMBER", + "CCSDS_PACKET_TYPE", + "CCSDS_SECONDARY_FLAG", + "CCSDS_SEQUENCE_FLAG", + "CCSDS_APID", + "CCSDS_SEQUENCE_COUNT", + "CCSDS_PACKET_LENGTH", + "CHECKSUM", + "time", + ] for this_col in colnames_to_remove: if this_col in hk_table.colnames: hk_table.remove_column(this_col) - empty_primary = fits.PrimaryHDU(header=primary_hdr) - hk_hdu = fits.BinTableHDU(hk_table, name="HK") + empty_primary_hdu = fits.PrimaryHDU(header=primary_hdr) + hk_hdu = fits.BinTableHDU(data=hk_table, name="HK") hk_hdu.add_checksum() - hdul = fits.HDUList([empty_primary, hk_hdu]) + + # add command response data if it exists + if parsed_data["cmd_resp"] is not None: + data_ts = parsed_data["cmd_resp"] + this_header = fits.Header() + this_header["DATEREF"] = ( + data_ts.time[0].fits, + get_std_comment("DATEREF"), + ) + record_timeseries(data_ts, "housekeeping") + data_table = Table(data_ts) + colnames_to_remove = [ + "CCSDS_VERSION_NUMBER", + "CCSDS_PACKET_TYPE", + "CCSDS_SECONDARY_FLAG", + "CCSDS_SEQUENCE_FLAG", + "CCSDS_APID", + "CCSDS_SEQUENCE_COUNT", + "CCSDS_PACKET_LENGTH", + "CHECKSUM", + "time", + ] + for this_col in colnames_to_remove: + if this_col in hk_table.colnames: + data_table.remove_column(this_col) + cmd_hdu = fits.BinTableHDU(data=data_table, name="READ") + cmd_hdu.add_checksum() + else: # if None still end an empty Binary Table + this_header = fits.Header() + cmd_hdu = fits.BinTableHDU(data=None, header=this_header, name="READ") + hdul = fits.HDUList([empty_primary_hdu, hk_hdu, cmd_hdu]) path = create_science_filename( time=date_beg, @@ -149,13 +176,8 @@ def process_file(filename: Path, overwrite=False) -> list: ) hdul.writeto(path, overwrite=overwrite) output_files.append(path) - - - - - # calibrated_file = calibrate_file(data_filename) - # data_plot_files = plot_file(data_filename) - # calib_plot_files = plot_file(calibrated_file) + if parsed_data["spectra"] is not None: + spec_data = parsed_data["spectra"] # add other tasks below return output_files diff --git a/padre_meddea/data/fits/fits_keywords_dict.csv b/padre_meddea/data/fits/fits_keywords_dict.csv new file mode 100644 index 0000000..60892bb --- /dev/null +++ b/padre_meddea/data/fits/fits_keywords_dict.csv @@ -0,0 +1,28 @@ +keyword,comment +AUTHOR,Who designed the observation +CREATOR,Name of software pipeline that produced the FITS file +DETECTOR,Name of the detector +INFO_URL,a human-readable web page describing the data +OBSRVTRY,Name of the observatory +TIMESYS,Time scale of the time-related keywords +TELESCOP,Telescope/Sensor name +INSTRUME,Instrument name +MISSION,Mission name +ORIGIN,File originator +DATE-BEG,Acquisition start time +DATE-END,Acquisition end time +DATE-AVG,Average time of acquisition +LEVEL,Data level of fits file +DATE,File creation date in UTC +DATATYPE,Description of the data +ORIGAPID,APID(s) of the originating data +ORIGFILE,Originating file(s) +DATEREF,Reference date +PRSTEP(?P[1-9]),Processing step type +PRPROC(?P[1-9]),Name of procedure performing PRSTEP +PRPVER(?P[1-9]),Version of procedure PRPROC +PRLIB(?P[1-9])A,Software library containing PRPROC +PRVER(?P[1-9])A,Version of PRLIB1A +PRHSH(?P[1-9])A,GIT commit hash for PRLIBA +PRBRA(?P[1-9])A,GIT/SVN repository branch of PRLIBA +PRVER(?P[1-9])B,Date of last commit of PRLIBB diff --git a/padre_meddea/data/fits/fits_keywords_primaryhdu.csv b/padre_meddea/data/fits/fits_keywords_primaryhdu.csv new file mode 100644 index 0000000..50ff863 --- /dev/null +++ b/padre_meddea/data/fits/fits_keywords_primaryhdu.csv @@ -0,0 +1,11 @@ +keyword,value +AUTHOR,Steven D. Christe +CREATOR,padre_meddea +DETECTOR,meddea +INFO_URL,https://padre-meddea.readthedocs.io/en/latest/user-guide/data.html +OBSRVTRY,PADRE +TIMESYS,UTC +TELESCOP,PADRE/MeDDEA +INSTRUME,MeDDEA +MISSION,PADRE +ORIGIN,NASA GSFC \ No newline at end of file diff --git a/padre_meddea/io/file_tools.py b/padre_meddea/io/file_tools.py index 1d3480f..2839701 100644 --- a/padre_meddea/io/file_tools.py +++ b/padre_meddea/io/file_tools.py @@ -82,17 +82,21 @@ def read_fits(filename: Path): Read a fits file. """ hdu = fits.open(filename) - + if (hdu[0].header["LEVEL"] == 0) and (hdu[0].header["DATATYPE"] == "event_list"): - event_list = read_fits_l0_event_list(filename) # do I need to close the file since it is being opened again right after this? + event_list = read_fits_l0_event_list( + filename + ) # do I need to close the file since it is being opened again right after this? return event_list + if (hdu[0].header["LEVEL"] == 0) and (hdu[0].header["DATATYPE"] == "housekeeping"): + hk_list = read_fits_l0_housekeeping(filename) + return hk_list else: - raise ValueError(F"File contents of {filename} not recogized.") + raise ValueError(f"File contents of {filename} not recogized.") def read_fits_l0_event_list(filename: Path) -> TimeSeries: - """ - """ + """ """ hdu = fits.open(filename) num_events = len(hdu["SCI"].data["seqcount"]) ph_times = calc_time( @@ -100,6 +104,7 @@ def read_fits_l0_event_list(filename: Path) -> TimeSeries: hdu["sci"].data["pktclock"], hdu["sci"].data["clocks"], ) + # add the pixel conversions pixel = np.array( [channel_to_pixel(this_chan) for this_chan in hdu["sci"].data["channel"]], dtype=np.uint8, @@ -119,6 +124,19 @@ def read_fits_l0_event_list(filename: Path) -> TimeSeries: return event_list +def read_fits_l0_housekeeping(filename: Path) -> TimeSeries: + """Read a level 0 housekeeping file + """ + hdu = fits.open(filename) + colnames = [this_col.name for this_col in hdu['HK'].data.columns] + times = calc_time(hdu["HK"].data["timestamp"]) + hk_list = TimeSeries( + time = times, + data = {key: hdu['hk'].data[key] for key in colnames} + ) + return hk_list + + def parse_ph_packets(filename: Path): """Given a binary file, read only the photon packets and return an event list. @@ -127,25 +145,25 @@ def parse_ph_packets(filename: Path): Photon packet format is (words are 16 bits) described below. - == ============================================================= - # Description - == ============================================================= - 0 CCSDS header 1 (0x00A0) - 1 CCSDS header 2 (0b11 and sequence count) - 2 CCSDS header 3 payload size (remaining packet size - 1 octet) - 3 time_stamp_s 1 - 4 time_stamp_s 2 - 5 time_stamp_clocks 1 - 6 time_stamp_clocks 2 - 7 integration time in clock counts - 8 live time in clock counts - 9 drop counter ([15] Int.Time Overflow, [14:12] decimation level, [11:0] # dropped photons) - 10 checksum - 11 start of pixel data - - pixel time step in clock count - - pixel_location (ASIC # bits[7:5], pixel num bits[4:0]) - - pixel_data 12 bit ADC count - == ============================================================= + == ============================================================= + # Description + == ============================================================= + 0 CCSDS header 1 (0x00A0) + 1 CCSDS header 2 (0b11 and sequence count) + 2 CCSDS header 3 payload size (remaining packet size - 1 octet) + 3 time_stamp_s 1 + 4 time_stamp_s 2 + 5 time_stamp_clocks 1 + 6 time_stamp_clocks 2 + 7 integration time in clock counts + 8 live time in clock counts + 9 drop counter ([15] Int.Time Overflow, [14:12] decimation level, [11:0] # dropped photons) + 10 checksum + 11 start of pixel data + - pixel time step in clock count + - pixel_location (ASIC # bits[7:5], pixel num bits[4:0]) + - pixel_data 12 bit ADC count + == ============================================================= Parameters ---------- @@ -211,8 +229,8 @@ def parse_ph_packets(filename: Path): ph_data["TIME_S"], ph_data["TIME_CLOCKS"], ): - num_hits = len(this_ph_data[0::3]) - ids = this_ph_data[1::3] + num_hits = len(this_ph_data[0::WORDS_PER_HIT]) + ids = this_ph_data[1::WORDS_PER_HIT] asic_num = (ids & 0b11100000) >> 5 channel_num = ids & 0b00011111 hit_list[0, i : i + num_hits] = this_ph_data[0::WORDS_PER_HIT] @@ -311,10 +329,8 @@ def parse_spectrum_packets(filename: Path): return None packet_definition = packet_definition_hist2() pkt = ccsdspy.FixedLength(packet_definition) - data = pkt.load(packet_bytes) - timestamps = [ - dt.timedelta(seconds=int(this_t)) + EPOCH for this_t in data["TIMESTAMPS"] - ] + data = pkt.load(packet_bytes, include_primary_header=True) + timestamps = calc_time(data["TIMESTAMPS"], data["TIMESTAMPCLOCK"]) num_packets = len(timestamps) h = data["HISTOGRAM_DATA"].reshape((num_packets, 24, 513)) histogram_data = h[:, :, 1:] # remove the id field @@ -325,6 +341,23 @@ def parse_spectrum_packets(filename: Path): def parse_cmd_response_packets(filename: Path): """Given a raw binary file, read only the command response packets and return the data. + The packet is defined as follows + + == ============================================================= + # Description + == ============================================================= + 0 CCSDS header 1 (0x00A5) + 1 CCSDS header 2 (0b11 and sequence count) + 2 CCSDS header 3 payload size (remaining packet size - 1 octet) + 3 time_stamp_s 1 + 4 time_stamp_s 2 + 5 time_stamp_clocks 1 + 6 time_stamp_clocks 2 + 7 register address + 8 register value + 9 checksum + == ============================================================= + Parameters ---------- filename : Path @@ -333,13 +366,26 @@ def parse_cmd_response_packets(filename: Path): Returns ------- cmd_resp_list : astropy.time.TimeSeries or None - A list of command responses.""" - return None - - -def packet_definition_cmd_resp(): - """Return the packet definiton for a command response packet.""" - pass + A list of register read responses. + """ + with open(filename, "rb") as mixed_file: + stream_by_apid = split_by_apid(mixed_file) + packet_bytes = stream_by_apid.get(APID["cmd_resp"], None) + if packet_bytes is None: + return None + packet_definition = packet_definition_cmd_response() + pkt = ccsdspy.FixedLength(packet_definition, include_primary_header=True) + data = pkt.load(packet_bytes) + timestamps = calc_time(data["TIMESTAMPS"], data["TIMESTAMPCLOCK"]) + data = { + "time_s": data["TIMESTAMPS"], + "time_clock": data["TIMESTAMPCLOCK"], + "address": data["ADDRESS"], + "value": data["VALUE"], + "seqcount": data["CCSDS_SEQUENCE_COUNT"], + } + ts = TimeSeries(time=timestamps, data=data) + return ts def packet_definition_hk(): @@ -430,6 +476,18 @@ def packet_definition_ph(): return p +def packet_definition_cmd_response(): + """Return the packet definiton for the register read response""" + p = [ + PacketField(name="TIMESTAMPS", data_type="uint", bit_length=32), + PacketField(name="TIMESTAMPCLOCK", data_type="uint", bit_length=32), + PacketField(name="ADDR", data_type="uint", bit_length=16), + PacketField(name="VALUE", data_type="uint", bit_length=16), + PacketField(name="CHECKSUM", data_type="uint", bit_length=16), + ] + return p + + def inspect_raw_file(filename: Path): """Given a raw binary file of packets, provide some high level summary information.""" with open(filename, "rb") as mixed_file: diff --git a/padre_meddea/io/fits_tools.py b/padre_meddea/io/fits_tools.py index 9dbb0f9..40c914c 100644 --- a/padre_meddea/io/fits_tools.py +++ b/padre_meddea/io/fits_tools.py @@ -1,6 +1,7 @@ """ This module provides a utilities to manage fits files reading and writing. """ + import re import git @@ -10,48 +11,52 @@ import padre_meddea -FITS_HDR0 = ascii.read(padre_meddea._data_directory / "fits" / "fits_keywords_primaryhdu.csv",format="csv") -FITS_HDR0.add_index('keyword') -FITS_HDR_KEYTOCOMMENT = ascii.read(padre_meddea._data_directory / "fits" / "fits_keywords_dict.csv",format="csv") -FITS_HDR_KEYTOCOMMENT.add_index('keyword') - - -def get_primary_header() -> fits.Header: - """Return a standard FITS file primary header.""" - header = fits.Header() - for row in FITS_HDR0: - this_comment = FITS_HDR_KEYTOCOMMENT.loc[row["keyword"]]['comment'] - header[row["keyword"]] = (row["value"], this_comment) - - return header +FITS_HDR0 = ascii.read( + padre_meddea._data_directory / "fits" / "fits_keywords_primaryhdu.csv", format="csv" +) +FITS_HDR0.add_index("keyword") +FITS_HDR_KEYTOCOMMENT = ascii.read( + padre_meddea._data_directory / "fits" / "fits_keywords_dict.csv", format="csv" +) +FITS_HDR_KEYTOCOMMENT.add_index("keyword") def get_std_comment(keyword: str) -> str: """Given a keyword, return the standard comment for a header card.""" - if keyword.upper() in FITS_HDR_KEYTOCOMMENT['keyword']: - return FITS_HDR_KEYTOCOMMENT.loc[keyword]['comment'] + if keyword.upper() in FITS_HDR_KEYTOCOMMENT["keyword"]: + return FITS_HDR_KEYTOCOMMENT.loc[keyword]["comment"] for this_row in FITS_HDR_KEYTOCOMMENT: - res = re.fullmatch(this_row['keyword'], keyword) + res = re.fullmatch(this_row["keyword"], keyword) if res: - comment = this_row['comment'] + comment = this_row["comment"] if len(res.groupdict()) > 0: # check if there was a match for key, value in res.groupdict().items(): comment = comment.replace(f"<{key}>", value) return comment +def get_primary_header() -> fits.Header: + """Return a standard FITS file primary header.""" + header = fits.Header() + header["DATE"] = (Time.now().fits, get_std_comment("DATE")) + for row in FITS_HDR0: + this_comment = get_std_comment(row["keyword"]) + header[row["keyword"]] = (row["value"], this_comment) + return header + + def add_process_info_to_header(header: fits.Header, n=1) -> fits.Header: """Add processing info metadata to a fits header. - + It adds the following SOLARNET compatible FITS cards; PRSTEPn, PRPROCn, PRPVERn, PRLIBnA, PRVERnA, PRLIBnA, PRHSHnA, PRVERnB - + Parameters ---------- header : fits.Header The fits header to add the new cards to n : int, default 1 - The processing step number. Must be greater than or equal to 1. + The processing step number. Must be >= 1 and <= 9. Returns ------- @@ -59,34 +64,34 @@ def add_process_info_to_header(header: fits.Header, n=1) -> fits.Header: """ if n < 1: ValueError("Processing number, n, must be greater than or equal to 1.") - header[f"PRSTEP{n}"] = ("PROCESS Raw to L1", get_std_comment(f'PRSTEP{n}')) + header[f"PRSTEP{n}"] = ("PROCESS Raw to L1", get_std_comment(f"PRSTEP{n}")) header[f"PRPROC{n}"] = ( "padre_meddea.calibration.process", - get_std_comment(f'PRPROC{n}'), + get_std_comment(f"PRPROC{n}"), ) header[f"PRPVER{n}"] = ( padre_meddea.__version__, - get_std_comment(f'PRPVER{n}'), + get_std_comment(f"PRPVER{n}"), ) header[f"PRLIB{n}A"] = ( "padre_meddea", - get_std_comment(f'PRLIB{n}A'), + get_std_comment(f"PRLIB{n}A"), ) - header[f"PRVER{n}A"] = (padre_meddea.__version__, get_std_comment(f'PRVER{n}A')) + header[f"PRVER{n}A"] = (padre_meddea.__version__, get_std_comment(f"PRVER{n}A")) repo = git.Repo(search_parent_directories=True) header[f"PRHSH{n}A"] = ( repo.head.object.hexsha, - get_std_comment(f'PRHSH{n}A'), + get_std_comment(f"PRHSH{n}A"), ) header[f"PRBRA{n}A"] = ( repo.active_branch.name, - get_std_comment(f'PRBRA{n}A'), + get_std_comment(f"PRBRA{n}A"), ) commits = list(repo.iter_commits("main", max_count=1)) header[f"PRVER{n}B"] = ( Time(commits[0].committed_datetime).fits, - get_std_comment(f'PRVER{n}B'), + get_std_comment(f"PRVER{n}B"), ) # primary_hdr["PRLOG1"] add log information, need to do this after the fact # primary_hdr["PRENV1"] add information about processing env, need to do this after the fact - return header \ No newline at end of file + return header diff --git a/padre_meddea/tests/test_fits_tools.py b/padre_meddea/tests/test_fits_tools.py index 7b988fd..e1252be 100644 --- a/padre_meddea/tests/test_fits_tools.py +++ b/padre_meddea/tests/test_fits_tools.py @@ -1,4 +1,5 @@ """Test for the fits_tools module""" + import pytest import astropy.io.fits as fits @@ -8,8 +9,8 @@ def test_comment_lookup_hdr0(): """Test that all keywords in fits_keyword_primaryhdu are listed in fits_keyword_dict""" - hdr0_keywords = list(FITS_HDR0['keyword']) - keyword_to_comment = list(FITS_HDR_KEYTOCOMMENT['keyword']) + hdr0_keywords = list(FITS_HDR0["keyword"]) + keyword_to_comment = list(FITS_HDR_KEYTOCOMMENT["keyword"]) for this_keyword in hdr0_keywords: assert this_keyword in keyword_to_comment @@ -32,10 +33,15 @@ def test_add_process_info_to_header(): assert len(this_card.value) > 0 assert len(this_card.comment) > 0 -@pytest.mark.parametrize("test_input,expected", [("PRSTEP2", "Processing step type"), - ("PRSTEP1", "Processing step type"), - ("PRPROC3", "Name of procedure performing PRSTEP3"), - ("PRHSH5A", "GIT commit hash for PRLIB5A") - ]) + +@pytest.mark.parametrize( + "test_input,expected", + [ + ("PRSTEP2", "Processing step type"), + ("PRSTEP1", "Processing step type"), + ("PRPROC3", "Name of procedure performing PRSTEP3"), + ("PRHSH5A", "GIT commit hash for PRLIB5A"), + ], +) def test_get_std_comment(test_input, expected): assert get_std_comment(test_input) == expected diff --git a/padre_meddea/util/util.py b/padre_meddea/util/util.py index b248fa4..d9db002 100644 --- a/padre_meddea/util/util.py +++ b/padre_meddea/util/util.py @@ -174,4 +174,5 @@ def has_baseline(filename: Path, packet_count=10) -> bool: num_hits = np.zeros(packet_count) for i in range(packet_count): num_hits[i] = (len(packet_bytes[i]) - HEADER_BYTES) / BYTES_PER_PHOTON - return np.sum(num_hits - np.floor(num_hits)) == 1 + # check if there is any remainder for non integer number of hits + return np.sum(num_hits - np.floor(num_hits)) == 0 From 5727a108654d076608206b359ce91e5d6d3cbd52 Mon Sep 17 00:00:00 2001 From: Steven Daniel Christe Date: Wed, 13 Nov 2024 16:01:47 -0500 Subject: [PATCH 10/21] bug fix for cmd_resp --- padre_meddea/__init__.py | 2 +- padre_meddea/io/file_tools.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/padre_meddea/__init__.py b/padre_meddea/__init__.py index 90ef540..0b0ffd0 100644 --- a/padre_meddea/__init__.py +++ b/padre_meddea/__init__.py @@ -59,7 +59,7 @@ "spectrum": 0xA2, # decimal 162 "photon": 0xA0, # decimal 160 "housekeeping": 0xA3, # decimal 163 - "cmd_resp": 0xA5, # decimal 165 + "cmd_resp": 0x99, # decimal 153 } EPOCH = Time("2000-01-01 00:00", scale="utc") diff --git a/padre_meddea/io/file_tools.py b/padre_meddea/io/file_tools.py index 2839701..fbeb197 100644 --- a/padre_meddea/io/file_tools.py +++ b/padre_meddea/io/file_tools.py @@ -374,13 +374,13 @@ def parse_cmd_response_packets(filename: Path): if packet_bytes is None: return None packet_definition = packet_definition_cmd_response() - pkt = ccsdspy.FixedLength(packet_definition, include_primary_header=True) - data = pkt.load(packet_bytes) + pkt = ccsdspy.FixedLength(packet_definition) + data = pkt.load(packet_bytes, include_primary_header=True) timestamps = calc_time(data["TIMESTAMPS"], data["TIMESTAMPCLOCK"]) data = { "time_s": data["TIMESTAMPS"], "time_clock": data["TIMESTAMPCLOCK"], - "address": data["ADDRESS"], + "address": data["ADDR"], "value": data["VALUE"], "seqcount": data["CCSDS_SEQUENCE_COUNT"], } From 564f1da538d20d87868fef7d4f718950e15e854c Mon Sep 17 00:00:00 2001 From: Steven Daniel Christe Date: Thu, 14 Nov 2024 09:23:22 -0500 Subject: [PATCH 11/21] Update util.py --- padre_meddea/util/util.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/padre_meddea/util/util.py b/padre_meddea/util/util.py index d9db002..a09e1b5 100644 --- a/padre_meddea/util/util.py +++ b/padre_meddea/util/util.py @@ -176,3 +176,8 @@ def has_baseline(filename: Path, packet_count=10) -> bool: num_hits[i] = (len(packet_bytes[i]) - HEADER_BYTES) / BYTES_PER_PHOTON # check if there is any remainder for non integer number of hits return np.sum(num_hits - np.floor(num_hits)) == 0 + + +def is_consecutive(arr: np.array): + """Return True if the array is all consecutive integers or has not missing numbers.""" + return np.all(np.diff(arr) == 1) \ No newline at end of file From be8715797183fc89734bb64764fcb9d87a8ea08a Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Sun, 17 Nov 2024 14:07:29 -0500 Subject: [PATCH 12/21] bug fix and protecting has_baseline --- padre_meddea/tests/test_file_tools.py | 5 +++-- padre_meddea/util/util.py | 16 ++++++++++------ 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/padre_meddea/tests/test_file_tools.py b/padre_meddea/tests/test_file_tools.py index 1088b38..1300b1c 100644 --- a/padre_meddea/tests/test_file_tools.py +++ b/padre_meddea/tests/test_file_tools.py @@ -27,11 +27,12 @@ def test_read_file_bad_file(): def test_read_ph_file(): - ph_list = parse_ph_packets(ph_packet_file) + ph_list, pkt_list = parse_ph_packets(ph_packet_file) # check that there are the right number of events # assert len(ph_list) == 760 # check that there are in fact 4 packets - assert len(np.unique(ph_list["packet_num"])) == NUM_PACKETS + assert len(np.unique(ph_list["seqcount"])) == NUM_PACKETS + assert len(pkt_list) == NUM_PACKETS def test_read_hk_file(): diff --git a/padre_meddea/util/util.py b/padre_meddea/util/util.py index a09e1b5..08991e4 100644 --- a/padre_meddea/util/util.py +++ b/padre_meddea/util/util.py @@ -14,7 +14,7 @@ import astropy.units as u from ccsdspy.utils import split_packet_bytes, split_by_apid -from padre_meddea import EPOCH +from padre_meddea import EPOCH, APID __all__ = ["create_science_filename", "has_baseline"] @@ -169,11 +169,15 @@ def has_baseline(filename: Path, packet_count=10) -> bool: with open(filename, "rb") as mixed_file: stream_by_apid = split_by_apid(mixed_file) - packet_stream = stream_by_apid[160] - packet_bytes = split_packet_bytes(packet_stream) - num_hits = np.zeros(packet_count) - for i in range(packet_count): - num_hits[i] = (len(packet_bytes[i]) - HEADER_BYTES) / BYTES_PER_PHOTON + if APID['photon'] in stream_by_apid.keys(): # only applicable to photon packets + packet_stream = stream_by_apid[APID['photon']] + packet_bytes = split_packet_bytes(packet_stream) + packet_count = min(len(packet_bytes), packet_count) # in case we have fewer than packet_count in the file + num_hits = np.zeros(packet_count) + for i in range(packet_count): + num_hits[i] = (len(packet_bytes[i]) - HEADER_BYTES) / BYTES_PER_PHOTON + else: + raise ValueError("Only works on photon packets.") # check if there is any remainder for non integer number of hits return np.sum(num_hits - np.floor(num_hits)) == 0 From 9ed516fad3684e665bcbe3d323fdad335d670ae5 Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Sun, 17 Nov 2024 15:01:20 -0500 Subject: [PATCH 13/21] fixed swxsoc config and added tests --- padre_meddea/__init__.py | 9 +++ padre_meddea/data/test/README.rst | 7 ++ padre_meddea/io/file_tools.py | 8 +- padre_meddea/tests/test_util.py | 100 +++++++++++++++++++++++++ padre_meddea/util/util.py | 118 +++++------------------------- 5 files changed, 138 insertions(+), 104 deletions(-) create mode 100644 padre_meddea/data/test/README.rst create mode 100644 padre_meddea/tests/test_util.py diff --git a/padre_meddea/__init__.py b/padre_meddea/__init__.py index b1ad58a..58360a4 100644 --- a/padre_meddea/__init__.py +++ b/padre_meddea/__init__.py @@ -23,6 +23,15 @@ # Load user configuration config = swxsoc_config +config["mission"]["mission_name"] = "padre" +config["mission"]["file_extension"] = ".fits" +config["mission"]["inst_names"] = ["meddea", "sharp"] +config["mission"]["inst_fullnames"] = [ + "Measuring Directivity to Determine Electron Anisotropy", + "sharp", +] +config["mission"]["inst_to_shortname"] = {"meddea": "meddea", "sharp": "sharp"} + log = swxsoc_log diff --git a/padre_meddea/data/test/README.rst b/padre_meddea/data/test/README.rst new file mode 100644 index 0000000..7f51d9c --- /dev/null +++ b/padre_meddea/data/test/README.rst @@ -0,0 +1,7 @@ +Test file directory +=================== + +apid[###]_4packets.bin +---------------------- +Contains 4 packets for each APID from telemetered packets. +Note that the photon packet file (160) does not contain baseline measurements. \ No newline at end of file diff --git a/padre_meddea/io/file_tools.py b/padre_meddea/io/file_tools.py index fbeb197..9a130af 100644 --- a/padre_meddea/io/file_tools.py +++ b/padre_meddea/io/file_tools.py @@ -125,14 +125,12 @@ def read_fits_l0_event_list(filename: Path) -> TimeSeries: def read_fits_l0_housekeeping(filename: Path) -> TimeSeries: - """Read a level 0 housekeeping file - """ + """Read a level 0 housekeeping file""" hdu = fits.open(filename) - colnames = [this_col.name for this_col in hdu['HK'].data.columns] + colnames = [this_col.name for this_col in hdu["HK"].data.columns] times = calc_time(hdu["HK"].data["timestamp"]) hk_list = TimeSeries( - time = times, - data = {key: hdu['hk'].data[key] for key in colnames} + time=times, data={key: hdu["hk"].data[key] for key in colnames} ) return hk_list diff --git a/padre_meddea/tests/test_util.py b/padre_meddea/tests/test_util.py new file mode 100644 index 0000000..1c663d8 --- /dev/null +++ b/padre_meddea/tests/test_util.py @@ -0,0 +1,100 @@ +import pytest + +import numpy as np + +import padre_meddea +from padre_meddea import EPOCH +import padre_meddea.util.util as util + +from astropy.time import TimeDelta +import astropy.units as u + +TIME = "2024-04-06T12:06:21" +TIME_FORMATTED = "20240406T120621" + + +# fmt: off +@pytest.mark.parametrize("instrument,time,level,version,result", [ + ("meddea", TIME, "l1", "1.2.3", f"padre_meddea_l1_{TIME_FORMATTED}_v1.2.3.fits"), + ("meddea", TIME, "l2", "2.4.5", f"padre_meddea_l2_{TIME_FORMATTED}_v2.4.5.fits"), + ("sharp", TIME, "l2", "1.3.5", f"padre_sharp_l2_{TIME_FORMATTED}_v1.3.5.fits"), + ("sharp", TIME, "l3", "2.4.5", f"padre_sharp_l3_{TIME_FORMATTED}_v2.4.5.fits"), +] +) +def test_science_filename_output_a(instrument, time, level, version, result): + """Test simple cases with expected output. + Since we are using the swxsoc create_science_filename, we are testing whether we did the config correctly in __init__.py""" + assert ( + util.create_science_filename(instrument, time, level=level, version=version) + == result + ) +# fmt: on + + +def test_is_consecutive(): + """Test if consecutive""" + assert util.is_consecutive(np.arange(10)) + assert util.is_consecutive(range(100)) + assert util.is_consecutive(np.arange(10, 100, 1)) + + +def test_is_not_consecutive(): + assert not util.is_consecutive([0, 2, 3, 4, 5]) + assert not util.is_consecutive( + np.concatenate((np.arange(1, 10), np.arange(11, 20))) + ) + + +@pytest.mark.parametrize( + "input,output", + [ + (26, 0), + (15, 1), + (8, 2), + (1, 3), + (29, 4), + (13, 5), + (5, 6), + (0, 7), + (30, 8), + (21, 9), + (11, 10), + (3, 11), + (31, 12), + (14, 14 + 12), # unconnected channel gets 12 added to it + ], +) +def test_channel_to_pix(input, output): + assert util.channel_to_pixel(input) == output + + +def test_has_baseline(): + assert not util.has_baseline( + padre_meddea._test_files_directory / "apid160_4packets.bin" + ) + + +def test_has_baseline_error(): + with pytest.raises(ValueError): + util.has_baseline(padre_meddea._test_files_directory / "apid162_4packets.bin") + + +@pytest.mark.parametrize( + "pkt_time_s,pkt_time_clk,ph_clk,output", + [ + (0, 0, 0, EPOCH), + (1, 0, 0, EPOCH + TimeDelta(1 * u.s)), + (10, 0, 0, EPOCH + TimeDelta(10 * u.s)), + (0, 1, 0, EPOCH + TimeDelta(0.05 * u.microsecond)), + (0, 0, 1, EPOCH + TimeDelta(12.8 * u.microsecond)), + ( + 5, + 5, + 5, + EPOCH + + TimeDelta(5 * u.s + 5 * 0.05 * u.microsecond + 5 * 12.8 * u.microsecond), + ), + ], +) +def test_calc_time(pkt_time_s, pkt_time_clk, ph_clk, output): + assert util.calc_time(pkt_time_s, pkt_time_clk, ph_clk) == output diff --git a/padre_meddea/util/util.py b/padre_meddea/util/util.py index 08991e4..e0845de 100644 --- a/padre_meddea/util/util.py +++ b/padre_meddea/util/util.py @@ -2,119 +2,37 @@ This module provides general utility functions. """ -import os -from datetime import datetime, timezone -import time from pathlib import Path import warnings import numpy as np - from astropy.time import Time, TimeDelta import astropy.units as u from ccsdspy.utils import split_packet_bytes, split_by_apid -from padre_meddea import EPOCH, APID - -__all__ = ["create_science_filename", "has_baseline"] - -TIME_FORMAT = "%Y%m%dT%H%M%S" -VALID_DATA_LEVELS = ["l0", "l1", "ql", "l2", "l3", "l4"] -FILENAME_EXTENSION = ".fits" - - -def create_science_filename( - time: str, - level: str, - version: str, - mode: str = "", - descriptor: str = "", - test: bool = False, -): - """Return a compliant filename. The format is defined as - - {mission}_{inst}_{mode}_{level}{test}_{descriptor}_{time}_v{version}.cdf - - This format is only appropriate for data level >= 1. - - Parameters - ---------- - instrument : `str` - The instrument name. Must be one of the following "eea", "nemesis", "merit", "spani" - time : `str` (in isot format) or ~astropy.time - The time - level : `str` - The data level. Must be one of the following "l0", "l1", "l2", "l3", "l4", "ql" - version : `str` - The file version which must be given as X.Y.Z - descriptor : `str` - An optional file descriptor. - mode : `str` - An optional instrument mode. - test : bool - Selects whether the file is a test file. - - Returns - ------- - filename : `str` - A CDF file name including the given parameters that matches the mission's file naming conventions - - Raises - ------ - ValueError: If the instrument is not recognized as one of the mission's instruments - ValueError: If the data level is not recognized as one of the mission's valid data levels - ValueError: If the data version does not match the mission's data version formatting conventions - ValueError: If the data product descriptor or instrument mode do not match the mission's formatting conventions - """ - test_str = "" - - if isinstance(time, str): - time_str = Time(time, format="isot").strftime(TIME_FORMAT) - else: - time_str = time.strftime(TIME_FORMAT) - - if level not in VALID_DATA_LEVELS[1:]: - raise ValueError( - f"Level, {level}, is not recognized. Must be one of {VALID_DATA_LEVELS[1:]}." - ) - # check that version is in the right format with three parts - if len(version.split(".")) != 3: - raise ValueError( - f"Version, {version}, is not formatted correctly. Should be X.Y.Z" - ) - # check that version has integers in each part - for item in version.split("."): - try: - int_value = int(item) - except ValueError: - raise ValueError(f"Version, {version}, is not all integers.") - - if test is True: - test_str = "test" - - # the parse_science_filename function depends on _ not being present elsewhere - if ("_" in mode) or ("_" in descriptor): - raise ValueError( - "The underscore symbol _ is not allowed in mode or descriptor." - ) +from swxsoc.util import create_science_filename - filename = ( - f"padre_meddea_{mode}_{level}{test_str}_{descriptor}_{time_str}_v{version}" - ) - filename = filename.replace("__", "_") # reformat if mode or descriptor not given +from padre_meddea import EPOCH, APID - return filename + FILENAME_EXTENSION +__all__ = [ + "create_science_filename", + "calc_time", + "has_baseline", + "is_consecutive", + "channel_to_pixel", +] -def calc_time(pkt_time_s, pkt_time_clk=0, ph_clk=0): +def calc_time(pkt_time_s, pkt_time_clk=0, ph_clk=0) -> Time: """ Convert times to a Time object """ deltat = TimeDelta( - pkt_time_s * u.s + pkt_time_clk * 0.05 * u.us + ph_clk * 12.8 * u.us + pkt_time_s * u.s + + pkt_time_clk * 0.05 * u.microsecond + + ph_clk * 12.8 * u.microsecond ) result = Time(EPOCH + deltat) - return result @@ -169,10 +87,12 @@ def has_baseline(filename: Path, packet_count=10) -> bool: with open(filename, "rb") as mixed_file: stream_by_apid = split_by_apid(mixed_file) - if APID['photon'] in stream_by_apid.keys(): # only applicable to photon packets - packet_stream = stream_by_apid[APID['photon']] + if APID["photon"] in stream_by_apid.keys(): # only applicable to photon packets + packet_stream = stream_by_apid[APID["photon"]] packet_bytes = split_packet_bytes(packet_stream) - packet_count = min(len(packet_bytes), packet_count) # in case we have fewer than packet_count in the file + packet_count = min( + len(packet_bytes), packet_count + ) # in case we have fewer than packet_count in the file num_hits = np.zeros(packet_count) for i in range(packet_count): num_hits[i] = (len(packet_bytes[i]) - HEADER_BYTES) / BYTES_PER_PHOTON @@ -184,4 +104,4 @@ def has_baseline(filename: Path, packet_count=10) -> bool: def is_consecutive(arr: np.array): """Return True if the array is all consecutive integers or has not missing numbers.""" - return np.all(np.diff(arr) == 1) \ No newline at end of file + return np.all(np.diff(arr) == 1) From 39d5545fe21ef7e66ab49673f64cb7df282b4319 Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Mon, 18 Nov 2024 07:45:07 -0500 Subject: [PATCH 14/21] Update fits_tools.py --- padre_meddea/io/fits_tools.py | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/padre_meddea/io/fits_tools.py b/padre_meddea/io/fits_tools.py index 40c914c..4eead03 100644 --- a/padre_meddea/io/fits_tools.py +++ b/padre_meddea/io/fits_tools.py @@ -4,6 +4,7 @@ import re import git +from git import InvalidGitRepositoryError from astropy.io import ascii import astropy.io.fits as fits @@ -78,20 +79,23 @@ def add_process_info_to_header(header: fits.Header, n=1) -> fits.Header: get_std_comment(f"PRLIB{n}A"), ) header[f"PRVER{n}A"] = (padre_meddea.__version__, get_std_comment(f"PRVER{n}A")) - repo = git.Repo(search_parent_directories=True) - header[f"PRHSH{n}A"] = ( - repo.head.object.hexsha, - get_std_comment(f"PRHSH{n}A"), - ) - header[f"PRBRA{n}A"] = ( - repo.active_branch.name, - get_std_comment(f"PRBRA{n}A"), - ) - commits = list(repo.iter_commits("main", max_count=1)) - header[f"PRVER{n}B"] = ( - Time(commits[0].committed_datetime).fits, - get_std_comment(f"PRVER{n}B"), - ) + try: + repo = git.Repo(padre_meddea.__file__, search_parent_directories=True) + header[f"PRHSH{n}A"] = ( + repo.head.object.hexsha, + get_std_comment(f"PRHSH{n}A"), + ) + header[f"PRBRA{n}A"] = ( + repo.active_branch.name, + get_std_comment(f"PRBRA{n}A"), + ) + commits = list(repo.iter_commits("main", max_count=1)) + header[f"PRVER{n}B"] = ( + Time(commits[0].committed_datetime).fits, + get_std_comment(f"PRVER{n}B"), + ) + except InvalidGitRepositoryError: + pass # primary_hdr["PRLOG1"] add log information, need to do this after the fact # primary_hdr["PRENV1"] add information about processing env, need to do this after the fact return header From 1fe918fd948222ca3500e4a6c4a621ec631cb6d4 Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Mon, 18 Nov 2024 07:53:01 -0500 Subject: [PATCH 15/21] Update fits_tools.py --- padre_meddea/io/fits_tools.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/padre_meddea/io/fits_tools.py b/padre_meddea/io/fits_tools.py index 4eead03..97ed69f 100644 --- a/padre_meddea/io/fits_tools.py +++ b/padre_meddea/io/fits_tools.py @@ -85,10 +85,10 @@ def add_process_info_to_header(header: fits.Header, n=1) -> fits.Header: repo.head.object.hexsha, get_std_comment(f"PRHSH{n}A"), ) - header[f"PRBRA{n}A"] = ( - repo.active_branch.name, - get_std_comment(f"PRBRA{n}A"), - ) + #header[f"PRBRA{n}A"] = ( + # repo.active_branch.name, + # get_std_comment(f"PRBRA{n}A"), + #) commits = list(repo.iter_commits("main", max_count=1)) header[f"PRVER{n}B"] = ( Time(commits[0].committed_datetime).fits, From 39865f83ce34d54ad070365ab8ebff7cb1d851fa Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Mon, 18 Nov 2024 08:13:08 -0500 Subject: [PATCH 16/21] Update fits_tools.py --- padre_meddea/io/fits_tools.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/padre_meddea/io/fits_tools.py b/padre_meddea/io/fits_tools.py index 97ed69f..473d593 100644 --- a/padre_meddea/io/fits_tools.py +++ b/padre_meddea/io/fits_tools.py @@ -85,15 +85,15 @@ def add_process_info_to_header(header: fits.Header, n=1) -> fits.Header: repo.head.object.hexsha, get_std_comment(f"PRHSH{n}A"), ) - #header[f"PRBRA{n}A"] = ( + # header[f"PRBRA{n}A"] = ( # repo.active_branch.name, # get_std_comment(f"PRBRA{n}A"), - #) - commits = list(repo.iter_commits("main", max_count=1)) - header[f"PRVER{n}B"] = ( - Time(commits[0].committed_datetime).fits, - get_std_comment(f"PRVER{n}B"), - ) + # ) + # commits = list(repo.iter_commits("main", max_count=1)) + # header[f"PRVER{n}B"] = ( + # Time(commits[0].committed_datetime).fits, + # get_std_comment(f"PRVER{n}B"), + # ) except InvalidGitRepositoryError: pass # primary_hdr["PRLOG1"] add log information, need to do this after the fact From 08501b46e586d45db4f0e17b3b1a6a3f9e9ba1aa Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Mon, 18 Nov 2024 08:15:18 -0500 Subject: [PATCH 17/21] Update calibration.py --- padre_meddea/calibration/calibration.py | 1 + 1 file changed, 1 insertion(+) diff --git a/padre_meddea/calibration/calibration.py b/padre_meddea/calibration/calibration.py index b8d8de7..fa89261 100644 --- a/padre_meddea/calibration/calibration.py +++ b/padre_meddea/calibration/calibration.py @@ -80,6 +80,7 @@ def process_file(filename: Path, overwrite=False) -> list: hdul = fits.HDUList([empty_primary_hdu, hit_hdu, pkt_hdu]) path = create_science_filename( + 'meddea', time=primary_hdr["DATE-BEG"], level="l1", descriptor="eventlist", From 5e056738d115ef3d662ff371fbe30b13ce2221e0 Mon Sep 17 00:00:00 2001 From: Steven Daniel Christe Date: Mon, 18 Nov 2024 12:35:08 -0500 Subject: [PATCH 18/21] removed unneeded config --- padre_meddea/__init__.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/padre_meddea/__init__.py b/padre_meddea/__init__.py index 58360a4..4a14847 100644 --- a/padre_meddea/__init__.py +++ b/padre_meddea/__init__.py @@ -21,18 +21,6 @@ print_config, ) -# Load user configuration -config = swxsoc_config -config["mission"]["mission_name"] = "padre" -config["mission"]["file_extension"] = ".fits" -config["mission"]["inst_names"] = ["meddea", "sharp"] -config["mission"]["inst_fullnames"] = [ - "Measuring Directivity to Determine Electron Anisotropy", - "sharp", -] -config["mission"]["inst_to_shortname"] = {"meddea": "meddea", "sharp": "sharp"} - - log = swxsoc_log # Then you can be explicit to control what ends up in the namespace, From 4400eddf886b46bd75ce51368bc04be4650d59db Mon Sep 17 00:00:00 2001 From: Steven Daniel Christe Date: Mon, 18 Nov 2024 12:42:15 -0500 Subject: [PATCH 19/21] made gitpython module optional --- padre_meddea/io/fits_tools.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/padre_meddea/io/fits_tools.py b/padre_meddea/io/fits_tools.py index 473d593..5d71f6c 100644 --- a/padre_meddea/io/fits_tools.py +++ b/padre_meddea/io/fits_tools.py @@ -3,8 +3,6 @@ """ import re -import git -from git import InvalidGitRepositoryError from astropy.io import ascii import astropy.io.fits as fits @@ -80,6 +78,9 @@ def add_process_info_to_header(header: fits.Header, n=1) -> fits.Header: ) header[f"PRVER{n}A"] = (padre_meddea.__version__, get_std_comment(f"PRVER{n}A")) try: + import git + from git import InvalidGitRepositoryError + repo = git.Repo(padre_meddea.__file__, search_parent_directories=True) header[f"PRHSH{n}A"] = ( repo.head.object.hexsha, @@ -94,6 +95,8 @@ def add_process_info_to_header(header: fits.Header, n=1) -> fits.Header: # Time(commits[0].committed_datetime).fits, # get_std_comment(f"PRVER{n}B"), # ) + except ModuleNotFoundError: + pass except InvalidGitRepositoryError: pass # primary_hdr["PRLOG1"] add log information, need to do this after the fact From c6b1fb858f7520af3a62c3f0b095539177616aa8 Mon Sep 17 00:00:00 2001 From: Steven Daniel Christe Date: Mon, 18 Nov 2024 12:46:57 -0500 Subject: [PATCH 20/21] black --- padre_meddea/calibration/calibration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/padre_meddea/calibration/calibration.py b/padre_meddea/calibration/calibration.py index fa89261..77145fe 100644 --- a/padre_meddea/calibration/calibration.py +++ b/padre_meddea/calibration/calibration.py @@ -80,7 +80,7 @@ def process_file(filename: Path, overwrite=False) -> list: hdul = fits.HDUList([empty_primary_hdu, hit_hdu, pkt_hdu]) path = create_science_filename( - 'meddea', + "meddea", time=primary_hdr["DATE-BEG"], level="l1", descriptor="eventlist", From 0bb8dc1d5c3c4e59521517bf537c1bf07e0d5751 Mon Sep 17 00:00:00 2001 From: Steven Daniel Christe Date: Mon, 18 Nov 2024 15:21:36 -0500 Subject: [PATCH 21/21] Update __init__.py --- padre_meddea/__init__.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/padre_meddea/__init__.py b/padre_meddea/__init__.py index 4a14847..b1ad58a 100644 --- a/padre_meddea/__init__.py +++ b/padre_meddea/__init__.py @@ -21,6 +21,9 @@ print_config, ) +# Load user configuration +config = swxsoc_config + log = swxsoc_log # Then you can be explicit to control what ends up in the namespace,