From fd55941b1e7eab68e153af617f61dec53d7c9090 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Tue, 29 Oct 2024 13:27:43 -0400 Subject: [PATCH] Add a get_aca_images method for mica l0 (#301) * Add a get_aca_images method for mica l0 * Handle instances when slot_data is empty / has no rows * Test more about empty table * Improve docstring * Implement full compliance with maude_decom API * Remove pyrightconfig.json * Ruff fixes * Fix docstring formatting * Exclude new IMG_VCDUCTR and IMG_TIME from col comparison * Remove IMG_TIME as it was cut * Add IMG_VCDUCTR column to match maude --------- Co-authored-by: Tom Aldcroft --- mica/archive/aca_l0.py | 161 +++++++++++++++++++++++++++++- mica/archive/tests/test_aca_l0.py | 71 ++++++++++++- 2 files changed, 228 insertions(+), 4 deletions(-) diff --git a/mica/archive/aca_l0.py b/mica/archive/aca_l0.py index 71b05015..4b74f6e4 100644 --- a/mica/archive/aca_l0.py +++ b/mica/archive/aca_l0.py @@ -17,9 +17,10 @@ import Ska.File import ska_dbi import tables -from astropy.table import Table +from astropy.table import Table, vstack from Chandra.Time import DateTime from chandra_aca.aca_image import ACAImage +from cxotime import CxoTimeLike from numpy import ma from mica.common import MICA_ARCHIVE, MissingDataError @@ -61,6 +62,8 @@ ("MNF", ">i4"), ("END_INTEG_TIME", ">f8"), ("INTEG", ">f4"), + ("PIXTLM", " np.ndarray: + """Extract UINT8's from a 2-D array of uint8 bit values into a new array. + + Input array is Nx8 where N is the number of rows and 8 is the number of bits. + It must be in little-endian order with the most significant bit in the rightmost + position of the array. + """ + # This returns a shape of (N, 1) so we need to squeeze it. + out = np.packbits(arr[:, bit_start:bit_stop], axis=1, bitorder="little") + return out.squeeze() + + +def get_aca_images( + start: CxoTimeLike, + stop: CxoTimeLike, +): + """ + Get ACA images from mica L0 archive in same format as maude_decom.get_aca_images. + + This gets a table of images from the mica L0 file archive that includes the same + output columns as ``chandra_aca.maude_decom.get_aca_images()``. The images are 8x8 + and are centered in the 8x8 image, with masks indicating empty pixels (for 4x4 or + 6x6 data). + + The units of the image data are DN, so you need to multiply by 5.0 / INTEG to get + e-/sec. + + This function returns additional columns that are not in the maude images: + - 'HD3TLM*': uint8, raw 8-bit header-3 telemetry values + - 'IMGSIZE': int32, 4, 6, or 8 + + Parameters + ---------- + start : CxoTimeLike + Start time of the time range. + stop : CxoTimeLike + Stop time of the time range. + + Returns + ------- + astropy.table.Table + ACA images and accompanying data. + """ + slot_data_list = [] + for slot in range(8): + slot_data_raw = get_slot_data( + start, + stop, + slot=slot, + centered_8x8=True, + ) + # Convert from numpy structured array to astropy table, keeping only good data + ok = slot_data_raw["QUALITY"] == 0 + slot_data = Table(slot_data_raw[ok]) + + # Add slot number if there are any rows (if statement not needed after + # https://github.com/astropy/astropy/pull/17102). + # if len(slot_data) > 0: + # slot_data["IMGNUM"] = slot + + for name in ["IMGFID", "IMGFUNC", "IMGNUM"]: + slot_data[name] = slot_data[f"{name}1"] + for ii in (1, 2, 3, 4): + name_ii = f"{name}{ii}" + if name_ii in slot_data.colnames: + del slot_data[name_ii] + + slot_data_list.append(slot_data) + + # Combine all slot data into one output table and sort by time and slot + out = vstack(slot_data_list) + out.sort(keys=["TIME", "IMGNUM"]) + + # Rename and rejigger colums to match chandra_aca.maude_decum.get_aca_images + is_6x6 = out["IMGSIZE"] == 6 + is_4x4 = out["IMGSIZE"] == 4 + for rc in ["ROW", "COL"]: + # IMGROW/COL0 are row/col of lower/left of 4x4, 6x6, or 8x8 image + # Make these new columns: + # IMGROW/COL_A1: row/col of pixel A1 of 4x4, 6x6, or 8x8 image + # IMGROW/COL0_8x8: row/col of lower/left of 8x8 image with smaller img centered + out[f"IMG{rc}_A1"] = out[f"IMG{rc}0"] + out[f"IMG{rc}_A1"][is_6x6] += 1 + out[f"IMG{rc}0_8X8"] = out[f"IMG{rc}0"] + out[f"IMG{rc}0_8X8"][is_6x6] -= 1 + out[f"IMG{rc}0_8X8"][is_4x4] -= 2 + + out["IMGTYPE"] = np.full(shape=len(out), fill_value=4) # 8x8 + out["IMGTYPE"][is_6x6] = 1 # 6x6 + out["IMGTYPE"][is_4x4] = 0 # 4x4 + + out.rename_column("IMGRAW", "IMG") + out.rename_column("BGDTYP", "AABGDTYP") + out.rename_column("PIXTLM", "AAPIXTLM") + del out["FILENAME"] + del out["QUALITY"] + + # maude_decom.get_aca_images() provides both VCDUCTR (VCDU for each sub-image) and + # IMG_VCDUCTR (VCDU of first sub-image). For combined images, these are identical. + out["VCDUCTR"] = out["MJF"] * 128 + out["MNF"] + out["IMG_VCDUCTR"] = out["VCDUCTR"] + + # Split uint8 GLBSTAT and IMGSTAT into individual bits. The stat_bit_names are + # in MSB order so we need to reverse them below. + for stat_name, stat_bit_names in ( + ("GLBSTAT", GLBSTAT_NAMES), + ("IMGSTAT", IMGSTAT_NAMES), + ): + for bit, name in zip(range(8), reversed(stat_bit_names)): + out[name] = (out[stat_name] & (1 << bit)) > 0 + + # Extract fields from the mica L0 8-bit COMMCNT value which is really three MSIDs. + # We need to use .data attribute of MaskedColumn to get a numpy masked array. The + # unpackbits function fails on MaskedColumn because it looks for a _mask attribute. + bits = np.unpackbits(out["COMMCNT"].data).reshape(-1, 8) + bits_rev = bits[:, ::-1] # MSB is on the right (index=7) + out["COMMCNT"] = _extract_bits_as_uint8(bits_rev, 2, 8) + out["COMMCNT_CHECKSUM_FAIL"] = bits_rev[:, 0] > 0 + out["COMMCNT_SYNTAX_ERROR"] = bits_rev[:, 1] > 0 + + # Extract fields from the mica L0 8-bit COMMPROG value which is really two MSIDs + bits = np.unpackbits(out["COMMPROG"].data).reshape(-1, 8) + bits_rev = bits[:, ::-1] + out["COMMPROG"] = _extract_bits_as_uint8(bits_rev, 2, 8) + out["COMMPROG_REPEAT"] = _extract_bits_as_uint8(bits_rev, 0, 2) + + return out + + def get_slot_data( start, stop, @@ -218,6 +373,10 @@ def get_slot_data( continue if fname in chunk.dtype.names: all_rows[fname][idx0:idx1] = chunk[fname] + elif fname == "PIXTLM": + all_rows[fname][idx0:idx1] = "ORIG" + elif fname == "BGDTYP": + all_rows[fname][idx0:idx1] = "FLAT" if "IMGSIZE" in columns: all_rows["IMGSIZE"][idx0:idx1] = f_imgsize if "FILENAME" in columns: diff --git a/mica/archive/tests/test_aca_l0.py b/mica/archive/tests/test_aca_l0.py index 0f6f5186..731a894b 100644 --- a/mica/archive/tests/test_aca_l0.py +++ b/mica/archive/tests/test_aca_l0.py @@ -3,10 +3,13 @@ import os +import astropy.units as u import numpy as np import pytest from astropy.table import Table -from Ska.Numpy import interpolate +from chandra_aca import maude_decom +from cxotime import CxoTime +from ska_numpy import interpolate from mica.archive import aca_l0, asp_l1 @@ -32,12 +35,73 @@ def test_l0_images_meta(): } +@pytest.mark.skipif(not has_l0_2012_archive, reason="Test requires 2012 L0 archive") +def test_get_aca_images(): + """ + Confirm mica and maude images sources agree for a small time range. + """ + # Includes 4x4, 6x6, 8x8 images, commanding, NPNT, NMAN, and fids + start = CxoTime("2012:270:02:25:00") + stop = CxoTime("2012:270:02:46:00") + imgs_maude = maude_decom.get_aca_images(start, stop) + imgs_maude.sort(keys=["TIME", "IMGNUM"]) + imgs_mica = aca_l0.get_aca_images(start, stop) + + assert len(imgs_maude) == len(imgs_mica) + + # Test data set includes commanding, searches, and all data types + assert set(imgs_mica["IMGSIZE"]) == {4, 6, 8} + assert set(imgs_mica["COMMPROG"]) == {0, 1, 3, 5, 7, 9, 11, 13, 15, 19, 21, 22} + assert set(imgs_mica["COMMCNT"]) == {0, 24, 13, 25} + assert set(imgs_mica["COMMPROG_REPEAT"]) == {0, 1} + + for colname in imgs_maude.colnames: + if imgs_maude[colname].dtype.kind == "f": + assert np.allclose( + imgs_mica[colname], imgs_maude[colname], rtol=0, atol=1e-3 + ) + else: + assert np.all(imgs_mica[colname] == imgs_maude[colname]) + + assert set(imgs_mica.colnames) - set(imgs_maude.colnames) == { + "HD3TLM63", + "HD3TLM77", + "HD3TLM66", + "HD3TLM74", + "HD3TLM64", + "HD3TLM76", + "IMGSIZE", + "HD3TLM65", + "HD3TLM72", + "HD3TLM73", + "HD3TLM75", + "HD3TLM62", + "HD3TLM67", + } + + +@pytest.mark.skipif(not has_l0_2012_archive, reason="Test requires 2012 L0 archive") +def test_get_aca_images_empty(): + """ + Confirm that get_aca_images returns a zero-length table when no images are found + """ + start = CxoTime("2012:270") + imgs_mica = aca_l0.get_aca_images(start, start) + imgs_maude = maude_decom.get_aca_images(start, start + 10 * u.s) + + # zero-length table with the required columns to match maude_decom + assert len(imgs_mica) == 0 + assert type(imgs_mica) is Table + assert set(imgs_maude.colnames).issubset(imgs_mica.colnames) + + has_l0_2007_archive = os.path.exists(os.path.join(aca_l0.CONFIG["data_root"], "2007")) has_asp_l1 = os.path.exists(os.path.join(asp_l1.CONFIG["data_root"])) @pytest.mark.skipif( - "not has_l0_2007_archive or not has_asp_l1", reason="Test requires 2007 L0 archive" + "not has_l0_2007_archive or not has_asp_l1", + reason="Test requires 2007 L0 archive and L1 aspect archive", ) def test_get_l0_images(): """ @@ -87,7 +151,8 @@ def test_get_l0_images(): @pytest.mark.skipif( - "not has_l0_2007_archive or not has_asp_l1", reason="Test requires 2007 L0 archive" + "not has_l0_2007_archive or not has_asp_l1", + reason="Test requires 2007 L0 archive and L1 aspect archive", ) def test_get_slot_data_8x8(): """