From c83669f05da9cfc1803f173e2b979bac7d495bae Mon Sep 17 00:00:00 2001 From: Steven Christe Date: Mon, 18 Nov 2024 15:33:14 -0500 Subject: [PATCH] Adds FITS file creation tools * improvement to photon fits files * add test for baseline * baseline aware file processing * adding more meta data to file files and new organization * adding metadata and reading fits files * added and re-organized documentation * added housekeeping l1 fits files * cleaning up fits tools * first working commit of fits file creation --- docs/user-guide/index.rst | 1 + docs/user-guide/level0.rst | 102 ++++-- docs/user-guide/raw.rst | 41 +++ padre_meddea/__init__.py | 14 +- padre_meddea/calibration/calibration.py | 146 ++++++-- padre_meddea/data/README.rst | 2 + padre_meddea/data/calibration/README.rst | 2 +- padre_meddea/data/fits/README.rst | 13 + padre_meddea/data/fits/fits_keywords_dict.csv | 28 ++ .../data/fits/fits_keywords_primaryhdu.csv | 11 + .../data/fits_keywords_primaryhdu.csv | 9 - padre_meddea/data/test/README.rst | 7 + padre_meddea/io/file_tools.py | 316 +++++++++++++----- padre_meddea/io/fits_tools.py | 104 ++++++ padre_meddea/tests/test_file_tools.py | 5 +- padre_meddea/tests/test_fits_tools.py | 47 +++ padre_meddea/tests/test_util.py | 100 ++++++ padre_meddea/util/util.py | 166 ++++----- pyproject.toml | 1 + 19 files changed, 862 insertions(+), 253 deletions(-) create mode 100644 docs/user-guide/raw.rst create mode 100644 padre_meddea/data/fits/README.rst create mode 100644 padre_meddea/data/fits/fits_keywords_dict.csv create mode 100644 padre_meddea/data/fits/fits_keywords_primaryhdu.csv delete mode 100644 padre_meddea/data/fits_keywords_primaryhdu.csv create mode 100644 padre_meddea/data/test/README.rst create mode 100644 padre_meddea/io/fits_tools.py create mode 100644 padre_meddea/tests/test_fits_tools.py create mode 100644 padre_meddea/tests/test_util.py 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..c4ec64d 100644 --- a/docs/user-guide/level0.rst +++ b/docs/user-guide/level0.rst @@ -6,36 +6,72 @@ 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 +======== ============================================================================================ ==== + +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. +It also includes any register read responses that may exist during that time period. 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/__init__.py b/padre_meddea/__init__.py index 7ecdc94..b1ad58a 100644 --- a/padre_meddea/__init__.py +++ b/padre_meddea/__init__.py @@ -2,6 +2,8 @@ import os from pathlib import Path +from astropy.time import Time + try: from ._version import version as __version__ from ._version import version_tuple @@ -31,9 +33,6 @@ _data_directory = _package_directory / "data" _test_files_directory = _package_directory / "data" / "test" -MISSION_NAME = "PADRE" -INSTRUMENT_NAME = "MeDDEA" - # the ratio of detector area for large pixels versus small pixels RATIO_TOTAL_LARGE_TO_SMALL_PIX = 0.947 @@ -63,4 +62,13 @@ 10.73, ] +APID = { + "spectrum": 0xA2, # decimal 162 + "photon": 0xA0, # decimal 160 + "housekeeping": 0xA3, # decimal 163 + "cmd_resp": 0x99, # decimal 153 +} + +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 a2cc766..77145fe 100644 --- a/padre_meddea/calibration/calibration.py +++ b/padre_meddea/calibration/calibration.py @@ -14,10 +14,15 @@ 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 +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", @@ -45,25 +50,39 @@ 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 - ph_list = parsed_data["photons"] - hdu = fits.PrimaryHDU(data=None) - hdu.header["DATE"] = (Time.now().fits, "FITS file creation date in UTC") - fits_meta = read_fits_keyword_file( - padre_meddea._data_directory / "fits_keywords_primaryhdu.csv" + 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"), ) - for row in fits_meta: - hdu.header[row["keyword"]] = (row["value"], row["comment"]) - bin_hdu = fits.BinTableHDU(data=Table(ph_list)) - hdul = fits.HDUList([hdu, bin_hdu]) + primary_hdr["ORIGFILE"] = (file_path.name, get_std_comment("ORIGFILE")) + + for this_keyword in ["DATE-BEG", "DATE-END", "DATE-AVG"]: + primary_hdr[this_keyword] = ( + event_list.meta.get(this_keyword, ""), + get_std_comment(this_keyword), + ) + + 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_hdu, hit_hdu, pkt_hdu]) path = create_science_filename( "meddea", - ph_list["time"][0].fits, - "l1", + time=primary_hdr["DATE-BEG"], + level="l1", descriptor="eventlist", test=True, version="0.1.0", @@ -77,21 +96,89 @@ 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 = [path] + output_files.append(path) if parsed_data["housekeeping"] is not None: hk_data = parsed_data["housekeeping"] - hk_data.meta["INSTRUME"] = "meddea" - - if "CHECKSUM" in hk_data.colnames: - hk_data.remove_column("CHECKSUM") - + # send data to AWS Timestream for Grafana dashboard record_timeseries(hk_data, "housekeeping") + hk_table = Table(hk_data) + + 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"), + ) + 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_hdu = fits.PrimaryHDU(header=primary_hdr) + hk_hdu = fits.BinTableHDU(data=hk_table, name="HK") + hk_hdu.add_checksum() + + # 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]) - # calibrated_file = calibrate_file(data_filename) - # data_plot_files = plot_file(data_filename) - # calib_plot_files = plot_file(calibrated_file) + 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) + if parsed_data["spectra"] is not None: + spec_data = parsed_data["spectra"] # add other tasks below return output_files @@ -146,12 +233,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/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/data/fits_keywords_primaryhdu.csv b/padre_meddea/data/fits_keywords_primaryhdu.csv deleted file mode 100644 index 1860e92..0000000 --- a/padre_meddea/data/fits_keywords_primaryhdu.csv +++ /dev/null @@ -1,9 +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 -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/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 7c48ea8..9a130af 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 @@ -16,18 +18,13 @@ count_packets, split_by_apid, ) +import astropy.io.fits as fits import padre_meddea +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): @@ -47,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 @@ -60,7 +57,7 @@ def read_raw_file(filename: Path): Read a level 0 data file. Parameters - ---------- + --------- filename : Path A file to read @@ -80,6 +77,64 @@ 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"): + 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.") + + +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"], + ) + # add the pixel conversions + 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 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. @@ -88,25 +143,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 ---------- @@ -122,67 +177,106 @@ 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 + # 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["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"]: - total_hits += len(this_packet[0::3]) - - hit_list = np.zeros((6, total_hits), dtype="uint16") - time_stamps = np.zeros(total_hits) + 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 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) otherwise 0 + + 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] + 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::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] 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] + 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 + 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_num": hit_list[1, :], - "asic_channel": hit_list[2, :], - "clock": hit_list[5, :], - "packet_num": 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["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 + 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, :], + # "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,14 +294,14 @@ 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, 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 @@ -228,15 +322,13 @@ 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) - 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 @@ -247,6 +339,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 @@ -255,20 +364,33 @@ 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) + 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["ADDR"], + "value": data["VALUE"], + "seqcount": data["CCSDS_SEQUENCE_COUNT"], + } + ts = TimeSeries(time=timestamps, data=data) + return ts 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)] @@ -352,6 +474,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 new file mode 100644 index 0000000..5d71f6c --- /dev/null +++ b/padre_meddea/io/fits_tools.py @@ -0,0 +1,104 @@ +""" +This module provides a utilities to manage fits files reading and writing. +""" + +import re + +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_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 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 >= 1 and <= 9. + + 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")) + 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, + 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 ModuleNotFoundError: + pass + 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 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/tests/test_fits_tools.py b/padre_meddea/tests/test_fits_tools.py new file mode 100644 index 0000000..e1252be --- /dev/null +++ b/padre_meddea/tests/test_fits_tools.py @@ -0,0 +1,47 @@ +"""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 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 2752c03..e0845de 100644 --- a/padre_meddea/util/util.py +++ b/padre_meddea/util/util.py @@ -2,100 +2,106 @@ 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 +from astropy.time import Time, TimeDelta +import astropy.units as u +from ccsdspy.utils import split_packet_bytes, split_by_apid + +from swxsoc.util import create_science_filename + +from padre_meddea import EPOCH, APID __all__ = [ "create_science_filename", + "calc_time", + "has_baseline", + "is_consecutive", + "channel_to_pixel", ] -TIME_FORMAT = "%Y%m%dT%H%M%S" -VALID_DATA_LEVELS = ["l0", "l1", "ql", "l2", "l3", "l4"] -FILENAME_EXTENSION = ".fits" +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.microsecond + + ph_clk * 12.8 * u.microsecond + ) + 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 create_science_filename( - instrument: str, - 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 +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). - This format is only appropriate for data level >= 1. + 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 ---------- - 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. + packet_bytes : byte string + Photon packet bytes, must be an integer number of whole packets and greaterh 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." - ) + result : bool - filename = ( - f"padre_meddea_{mode}_{level}{test_str}_{descriptor}_{time_str}_v{version}" - ) - filename = filename.replace("__", "_") # reformat if mode or descriptor not given - - return filename + FILENAME_EXTENSION + """ + 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) + 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 + + +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) diff --git a/pyproject.toml b/pyproject.toml index 4d7262d..4e40be7 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]