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 <taldcroft@gmail.com>

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", "<U4"),
+    ("BGDTYP", "<U4"),
     ("GLBSTAT", "|u1"),
     ("COMMCNT", "|u1"),
     ("COMMPROG", "|u1"),
@@ -118,6 +121,27 @@
     cda_table="cda_aca0.h5",
 )
 
+# Names of bits in GLBSTAT and IMGSTAT fields, in order from most to least significant
+GLBSTAT_NAMES = [
+    "HIGH_BGD",
+    "RAM_FAIL",
+    "ROM_FAIL",
+    "POWER_FAIL",
+    "CAL_FAIL",
+    "COMM_CHECKSUM_FAIL",
+    "RESET",
+    "SYNTAX_ERROR",
+]
+
+IMGSTAT_NAMES = [
+    "SAT_PIXEL",
+    "DEF_PIXEL",
+    "QUAD_BOUND",
+    "COMMON_COL",
+    "MULTI_STAR",
+    "ION_RAD",
+]
+
 
 def get_options():
     parser = argparse.ArgumentParser(
@@ -144,6 +168,137 @@ def get_options():
     return opt
 
 
+def _extract_bits_as_uint8(
+    arr: np.ndarray, bit_start: int, bit_stop: int
+) -> 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:

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():
     """