From a7230c6ff25971edef6a75011980e3f118478766 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 28 Aug 2024 13:16:23 -0400 Subject: [PATCH 01/28] Add some routines to get background subtracted images --- chandra_aca/dark_subtract.py | 286 +++++++++++++++++++++++ chandra_aca/tests/test_dark_subtract.py | 295 ++++++++++++++++++++++++ 2 files changed, 581 insertions(+) create mode 100644 chandra_aca/dark_subtract.py create mode 100644 chandra_aca/tests/test_dark_subtract.py diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py new file mode 100644 index 00000000..926731c5 --- /dev/null +++ b/chandra_aca/dark_subtract.py @@ -0,0 +1,286 @@ +from functools import lru_cache + +import numba +import numpy as np +from astropy.table import Table, vstack +from cheta import fetch_sci +from cxotime import CxoTime, CxoTimeLike +from mica.archive import aca_l0 +from mica.archive.aca_dark import get_dark_cal_props + +from chandra_aca import maude_decom +from chandra_aca.dark_model import dark_temp_scale_img + + +@lru_cache(maxsize=2) +def get_mica_images(start: CxoTimeLike, stop: CxoTimeLike): + """ + Get ACA images from MICA for a given time range. + + Parameters + ---------- + start : CxoTimeLike + Start time of the time range. + stop : CxoTimeLike + Stop time of the time range. + + Returns + ------- + dict + dict with a astropy.table of ACA image data for each slot. + + """ + slot_data = {} + for slot in range(8): + slot_data[slot] = aca_l0.get_slot_data( + start, + stop, + imgsize=[4, 6, 8], + slot=slot, + columns=[ + "TIME", + "IMGRAW", + "IMGSIZE", + "IMGROW0", + "IMGCOL0", + ], + centered_8x8=True, + ) + slot_data[slot] = Table(slot_data[slot]) + + # Rename and rejigger colums to match maude + IMGROW0_8X8 = slot_data[slot]["IMGROW0"].copy() + IMGCOL0_8X8 = slot_data[slot]["IMGCOL0"].copy() + IMGROW0_8X8[slot_data[slot]["IMGSIZE"] == 6] -= 1 + IMGCOL0_8X8[slot_data[slot]["IMGSIZE"] == 6] -= 1 + IMGROW0_8X8[slot_data[slot]["IMGSIZE"] == 4] -= 2 + IMGCOL0_8X8[slot_data[slot]["IMGSIZE"] == 4] -= 2 + slot_data[slot]["IMGROW0_8X8"] = IMGROW0_8X8 + slot_data[slot]["IMGCOL0_8X8"] = IMGCOL0_8X8 + slot_data[slot].rename_column("IMGRAW", "IMG") + + return slot_data + + +@lru_cache(maxsize=2) +def get_maude_images(start: CxoTimeLike, stop: CxoTimeLike, fetch_limit_hours=100): + """ + Get ACA images from MAUDE for a given time range. + + Parameters + ---------- + start : CxoTimeLike + Start time of the time range. + stop : CxoTimeLike + Stop time of the time range. + fetch_limit_hours : int + Maximum number of hours to fetch from MAUDE in one go. + + Returns + ------- + dict + dict with a astropy.table of ACA image data for each slot. + """ + tstart = CxoTime(start).secs + tstop = CxoTime(stop).secs + + if tstop - tstart > fetch_limit_hours * 60 * 60: + raise ValueError("Time range too large for maude fetch, limit is 100 hours") + + slot_chunks = [] + # Break maude fetches into max 3 hour chunks required by maude_decom fetch + for mstart in np.arange(tstart, tstop, 60 * 60 * 3): + mstop = np.min([tstop, mstart + 60 * 60 * 3]) + imgs = maude_decom.get_aca_images(mstart, mstop) + slot_chunks.append(imgs) + slot_chunks = vstack(slot_chunks) + + slot_data = {} + for slot in range(8): + slot_data[slot] = slot_chunks[slot_chunks["IMGNUM"] == slot] + + return slot_data + + +def get_dcsub_aca_images( + start=None, + stop=None, + aca_images=None, + t_ccd=None, + t_ccd_times=None, + dc_img=None, + dc_tccd=None, + source="maude", +): + """ + Get dark current subtracted ACA images for a given time range. + + Parameters + ---------- + start : str, optional + Start time of the time range. Default is None. + stop : str, optional + Stop time of the time range. Default is None. + aca_images : dict, optional + Dictionary of ACA images. Should be keyed by slot. Default is None. + t_ccd : array-like, optional + Array of CCD temperatures. Default is None. + t_ccd_times : array-like, optional + Array of CCD temperature times. Default is None. + dc_img : array-like, optional + Dark current image. Default is None. + dc_tccd : float, optional + Dark current CCD temperature. Default is None. + source : str, optional + Source of the ACA images. Can be 'maude' or 'mica'. Default is 'maude'. + + Returns + ------- + dict + Dictionary of dark current subtracted ACA images for each slot. + + Raises + ------ + ValueError + If the dark calibration does not have t_ccd or ccd_temp. + If aca_image_table is not defined or source is not 'maude' or 'mica'. + + """ + if dc_img is not None: + assert dc_img.shape == (1024, 1024) + + if dc_img is None: + dc = get_dark_cal_props(start, "nearest", include_image=True, aca_image=False) + dc_img = dc["image"] + if "t_ccd" in dc: + dc_tccd = dc["t_ccd"] + elif "ccd_temp" in dc: + dc_tccd = dc["ccd_temp"] + else: + raise ValueError("Dark cal must have t_ccd or ccd_temp") + + if t_ccd is None: + if source == "maude": + with fetch_sci.data_source("maude allow_subset=False"): + dat = fetch_sci.Msid("aacccdpt", start, stop) + else: + with fetch_sci.data_source("cxc"): + dat = fetch_sci.Msid("aacccdpt", start, stop) + t_ccd = dat.vals + t_ccd_times = dat.times + + if aca_images is not None: + # If supplied a dictionary of aca images, use that. Should be keyed by slot. + # This is for testing. + assert isinstance(aca_images, dict) + # For each slot, require aca_image have columns IMG, IMGROW0_8X8, IMGCOL0_8X8, TIME + for slot in aca_images: + assert "IMG" in aca_images[slot].colnames + assert "IMGROW0_8X8" in aca_images[slot].colnames + assert "IMGCOL0_8X8" in aca_images[slot].colnames + assert "TIME" in aca_images[slot].colnames + # require that IMGROW0_8X8 and IMGCOL0_8X8 are masked columns + assert hasattr(aca_images[slot]["IMGROW0_8X8"], "mask") + assert hasattr(aca_images[slot]["IMGCOL0_8X8"], "mask") + imgs = aca_images + elif source == "maude": + imgs = get_maude_images(start, stop) + elif source == "mica": + imgs = get_mica_images(start, stop) + else: + raise ValueError( + "aca_image_table must be defined or source must be maude or mica" + ) + + imgs_bgsub = {} + dark_imgs = {} + for slot in range(8): + if slot not in imgs: + continue + dark_imgs[slot] = get_dark_current_imgs( + imgs[slot], dc_img, dc_tccd, t_ccd, t_ccd_times + ) + img_col = "IMG" + imgs_bgsub[slot] = imgs[slot][img_col] - dark_imgs[slot] + imgs_bgsub[slot].clip(0, None) + + return imgs_bgsub + + +def get_dark_current_imgs(img_table, dc_img, dc_tccd, t_ccd, t_ccd_times): + """ + Get the scaled dark current values for a table of ACA images. + + Parameters + ---------- + img_table : astry.table.Table + A table containing information about the images. + dc_img : 1024x1024 array + The dark calibration image. + dc_tccd : float + The reference temperature of the dark calibration image. + t_ccd : array + The temperature values for the time range of the img_table. + t_ccd_times : type + The corresponding times for the temperature values. + + Returns + ------- + dark : np.array + An array containing the temperature scaled dark current for each ACA image + in the img_table in e-/s. + + """ + assert dc_img.shape == (1024, 1024) + assert len(t_ccd) == len(t_ccd_times) + + dark_raw = get_dark_backgrounds( + dc_img, + img_table["IMGROW0_8X8"].data.filled(1025), + img_table["IMGCOL0_8X8"].data.filled(1025), + ) + dt_ccd = np.interp(img_table["TIME"], t_ccd_times, t_ccd) + + dark = np.zeros((len(dark_raw), 8, 8), dtype=np.float64) + + # Scale the dark current at each dark cal 8x8 image to the ccd temperature and e-/s + for i, (eight, t_ccd_i) in enumerate(zip(dark_raw, dt_ccd, strict=True)): + img_sc = dark_temp_scale_img(eight, dc_tccd, t_ccd_i) + # Note that this just uses standard integration time of 1.696 sec + dark[i] = img_sc * 1.696 / 5 + + return dark + + +def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): + """ + Get the dark background for a stack/table of ACA image. + + Parameters + ---------- + raw_dark_img : np.array + The dark calibration image. + imgrow0 : np.array + The row of the ACA image. + imgcol0 : np.array + The column of the ACA image. + size : int, optional + The size of the ACA image. Default is 8. + + Returns + ------- + dark_img : np.array + The dark backgrounds for the ACA images. + """ + + @numba.jit(nopython=True) + def staggered_aca_slice(array_in, array_out, row, col): + for i in np.arange(len(row)): + if row[i] + size < 1024 and col[i] + size < 1024: + array_out[i] = array_in[row[i] : row[i] + size, col[i] : col[i] + size] + + dark_img = np.zeros([len(imgrow0), size, size], dtype=np.float64) + staggered_aca_slice( + raw_dark_img.astype(float), dark_img, 512 + imgrow0, 512 + imgcol0 + ) + return dark_img diff --git a/chandra_aca/tests/test_dark_subtract.py b/chandra_aca/tests/test_dark_subtract.py new file mode 100644 index 00000000..3470f442 --- /dev/null +++ b/chandra_aca/tests/test_dark_subtract.py @@ -0,0 +1,295 @@ +from pathlib import Path + +import mica.common +import numpy as np +import pytest +from astropy.table import Table +from mica.archive.aca_dark import get_dark_cal_props +from numpy import ma + +from chandra_aca.dark_subtract import ( + get_dark_backgrounds, + get_dark_current_imgs, + get_dcsub_aca_images, + get_maude_images, + get_mica_images, +) + +HAS_ACA0_ARCHIVE = (Path(mica.common.MICA_ARCHIVE) / "aca0").exists() + + +# Make a test fixture for the mock dark current image +@pytest.fixture +def mock_dc(): + """ + Make a mock dark current image with some values in the corner. + """ + mock_dc = np.zeros((1024, 1024)) + # put some ints in the corner + for i in range(16): + for j in range(16): + mock_dc[i, j] = i + j + return mock_dc + + +# Make a test fixture for the mock 8x8 data +@pytest.fixture +def mock_imgs(): + """ + Make a mock 8x8 image table. + + The intent is that these image coordinates should overlap with the mock dark current. + """ + imgs = {} + img_table = Table() + img_table["TIME"] = np.arange(8) + img_table["IMGROW0_8X8"] = ma.arange(8) - 512 + img_table["IMGCOL0_8X8"] = ma.arange(8) - 512 + img_table["IMG"] = np.ones((8, 8, 8)) * 16 * 1.696 / 5 + imgs[0] = img_table + return imgs + + +# Test fixture array of dark current images in DN +@pytest.fixture +def dc_imgs_dn(): + """ + Save expected results of extracted dark current images in this array. + """ + dc_imgs_dn = np.array( + [ + [ + [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0], + [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0], + [2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0], + [3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], + [4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0], + [5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0], + [6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0], + [7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0], + ], + [ + [2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0], + [3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], + [4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0], + [5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0], + [6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0], + [7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0], + [8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0], + [9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0], + ], + [ + [4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0], + [5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0], + [6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0], + [7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0], + [8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0], + [9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0], + [10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0], + [11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0], + ], + [ + [6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0], + [7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0], + [8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0], + [9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0], + [10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0], + [11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0], + [12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0], + [13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0], + ], + [ + [8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0], + [9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0], + [10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0], + [11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0], + [12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0], + [13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0], + [14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0], + [15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0], + ], + [ + [10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0], + [11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0], + [12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0], + [13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0], + [14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0], + [15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0], + [16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0], + [17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0], + ], + [ + [12.0, 13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0], + [13.0, 14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0], + [14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0], + [15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0], + [16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0], + [17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0], + [18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0], + [19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0], + ], + [ + [14.0, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0], + [15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0], + [16.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0], + [17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0], + [18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0], + [19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0], + [20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0], + [21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0], + ], + ] + ) + return dc_imgs_dn + + +def test_dcsub_aca_images(mock_dc, mock_imgs, dc_imgs_dn): + """ + Confirm that the pattern of background subtraction is correct. + """ + + imgs_bgsub = get_dcsub_aca_images( + aca_images=mock_imgs, + dc_img=mock_dc, + dc_tccd=-10, + t_ccd=np.repeat(-10, 8), + t_ccd_times=mock_imgs[0]["TIME"], + ) + assert imgs_bgsub[0].shape == (8, 8, 8) + # Note that the mock unsubtracted data is originally 16 * 1.696 / 5 + exp = 16 - dc_imgs_dn + assert np.allclose(imgs_bgsub[0] * 5 / 1.696, exp, atol=1e-6, rtol=0) + + +def test_get_dark_images(mock_dc, mock_imgs, dc_imgs_dn): + """ + Confirm the pattern of dark current images matches the reference. + The temperature scaling is set here to be a no-op. + """ + dc_imgs = get_dark_current_imgs( + mock_imgs[0], + dc_img=mock_dc, + dc_tccd=-10, + t_ccd=np.repeat(-10, 8), + t_ccd_times=mock_imgs[0]["TIME"], + ) + assert dc_imgs.shape == (8, 8, 8) + dc_imgs_raw = dc_imgs * 5 / 1.696 + assert np.allclose(dc_imgs_raw, dc_imgs_dn, atol=1e-6, rtol=0) + + +def test_get_dark_backgrounds(mock_dc, mock_imgs, dc_imgs_dn): + """ + Confirm the pattern of dark current matches the reference. + """ + dc_imgs_raw = get_dark_backgrounds( + mock_dc, + mock_imgs[0]["IMGROW0_8X8"].data.filled(1025), + mock_imgs[0]["IMGCOL0_8X8"].data.filled(1025), + ) + assert dc_imgs_raw.shape == (8, 8, 8) + assert np.allclose(dc_imgs_raw, dc_imgs_dn, atol=1e-6, rtol=0) + + +@pytest.mark.skipif(not HAS_ACA0_ARCHIVE, reason="No ACA0 archive") +def test_get_images(): + """ + Confirm mica and maude images sources agree for a small time range. + """ + tstart = 746484519.429 + imgs_maude = get_maude_images(tstart, tstart + 30) + imgs_mica = get_mica_images(tstart, tstart + 30) + for slot in range(8): + assert len(imgs_maude[slot]) == len(imgs_mica[slot]) + for col in ["IMG", "IMGROW0_8X8", "IMGCOL0_8X8", "TIME"]: + assert np.allclose( + imgs_maude[slot][col], imgs_mica[slot][col], atol=0.1, rtol=0 + ) + + +def test_dc_consistent(): + """ + Confirm extracted dark current is reasonably consistent with the images + edges of a 6x6. + """ + tstart = 746484519.429 + imgs_maude = get_maude_images(tstart, tstart + 30) + # This slot and this time show an overall decent correlation of dark current + # and edge pixel values in the 6x6 image - but don't make a great test. + slot = 4 + img = imgs_maude[slot]["IMG"][0] + + dc = get_dark_cal_props(tstart, "nearest", include_image=True, aca_image=False) + dc_imgs = get_dark_current_imgs( + imgs_maude[slot], + dc_img=dc["image"], + dc_tccd=dc["t_ccd"], + t_ccd=np.repeat(-9.9, len(imgs_maude[slot])), + t_ccd_times=imgs_maude[slot]["TIME"], + ) + edge_mask_6x6 = np.zeros((8, 8), dtype=bool) + edge_mask_6x6[1, 2:6] = True + edge_mask_6x6[6, 2:6] = True + edge_mask_6x6[2:6, 1] = True + edge_mask_6x6[2:6, 6] = True + + # Here just show that a lot of pixels are within 50 e-/s of the dark current + # In the 6x6, star spoiling is a large effect. + assert ( + np.percentile(img[edge_mask_6x6].filled(0) - dc_imgs[0][edge_mask_6x6], 50) < 50 + ) + + +def test_dcsub_aca_images_maude(): + """ + Test that dark current background subtracted images match + references saved in the test. + """ + # This one is just a regression test + tstart = 471139130.93277466 # dwell tstart for obsid 15600 - random obsid + tstop = tstart + 20 + imgs_bgsub = get_dcsub_aca_images(tstart, tstop, source="maude") + + exp0 = np.array( + [ + [7, 8, 22, 49, 71, 75, 23, 17], + [22, 31, 33, 75, 125, 83, 54, 24], + [28, 50, 166, 367, 491, 311, 144, 45], + [48, 152, 748, 4124, 5668, 1548, 435, 94], + [50, 202, 966, 5895, 8722, 2425, 615, 141], + [29, 88, 331, 999, 1581, 631, 382, 103], + [16, 36, 108, 279, 302, 229, 112, 62], + [14, 14, 33, 69, 71, 49, 33, 25], + ] + ) + np.allclose(exp0, imgs_bgsub[0][0].filled(0), atol=1, rtol=0) + + # slot 4 is 6x6 and masked + exp4 = np.array( + [ + [0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 35, 41, 51, 63, 0, 0], + [0, 84, 301, 572, 498, 124, 56, 0], + [0, 311, 1982, 4991, 4314, 475, 116, 0], + [0, 773, 4551, 3157, 4647, 687, 100, 0], + [0, 617, 2586, 3120, 2617, 296, 55, 0], + [0, 0, 425, 1372, 344, 70, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0], + ] + ) + np.allclose(exp4, imgs_bgsub[4][0].filled(0), atol=1, rtol=0) + + +@pytest.mark.skipif(not HAS_ACA0_ARCHIVE, reason="No ACA0 archive") +def test_dcsub_aca_images_mica_maude(): + """ + Confirm background subtracted mica/maude images match for small time range. + """ + tstart = 471139130.93277466 # dwell tstart for obsid 15600 - random obsid + tstop = tstart + 30 + imgs_bgsub = get_dcsub_aca_images(tstart, tstop, source="mica") + + imgs_bgsub_maude = get_dcsub_aca_images(tstart, tstop, source="maude") + + for slot in range(8): + assert np.allclose(imgs_bgsub[slot], imgs_bgsub_maude[slot], atol=1e-3, rtol=0) From ab2a8f15103e4f8445f7d2e63607f56473bc3efb Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Tue, 17 Sep 2024 11:11:49 -0400 Subject: [PATCH 02/28] Break into more subroutines --- chandra_aca/dark_subtract.py | 319 +++++++++++++++++------- chandra_aca/tests/test_dark_subtract.py | 54 ++++ 2 files changed, 283 insertions(+), 90 deletions(-) diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index 926731c5..025b0508 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -12,60 +12,40 @@ from chandra_aca.dark_model import dark_temp_scale_img -@lru_cache(maxsize=2) -def get_mica_images(start: CxoTimeLike, stop: CxoTimeLike): +def get_dark_data(time: CxoTimeLike): """ - Get ACA images from MICA for a given time range. + Retrieve the dark calibration image and temperature for a given time. + + This is a light wrapper around mica.archive.aca_dark.get_dark_cal_props. Parameters ---------- - start : CxoTimeLike - Start time of the time range. - stop : CxoTimeLike - Stop time of the time range. + time : CxoTimeLike + Reference time for the dark calibration image. Returns ------- - dict - dict with a astropy.table of ACA image data for each slot. - + dc_img : np.array + Dark calibration image. + dc_tccd : float + Dark calibration CCD temperature. """ - slot_data = {} - for slot in range(8): - slot_data[slot] = aca_l0.get_slot_data( - start, - stop, - imgsize=[4, 6, 8], - slot=slot, - columns=[ - "TIME", - "IMGRAW", - "IMGSIZE", - "IMGROW0", - "IMGCOL0", - ], - centered_8x8=True, - ) - slot_data[slot] = Table(slot_data[slot]) - - # Rename and rejigger colums to match maude - IMGROW0_8X8 = slot_data[slot]["IMGROW0"].copy() - IMGCOL0_8X8 = slot_data[slot]["IMGCOL0"].copy() - IMGROW0_8X8[slot_data[slot]["IMGSIZE"] == 6] -= 1 - IMGCOL0_8X8[slot_data[slot]["IMGSIZE"] == 6] -= 1 - IMGROW0_8X8[slot_data[slot]["IMGSIZE"] == 4] -= 2 - IMGCOL0_8X8[slot_data[slot]["IMGSIZE"] == 4] -= 2 - slot_data[slot]["IMGROW0_8X8"] = IMGROW0_8X8 - slot_data[slot]["IMGCOL0_8X8"] = IMGCOL0_8X8 - slot_data[slot].rename_column("IMGRAW", "IMG") - - return slot_data + dc = get_dark_cal_props(time, "nearest", include_image=True, aca_image=False) + dc_img = dc["image"] + if "t_ccd" in dc: + dc_tccd = dc["t_ccd"] + elif "ccd_temp" in dc: + dc_tccd = dc["ccd_temp"] + else: + raise ValueError("Dark cal must have t_ccd or ccd_temp") + return dc_img, dc_tccd -@lru_cache(maxsize=2) -def get_maude_images(start: CxoTimeLike, stop: CxoTimeLike, fetch_limit_hours=100): +def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): """ - Get ACA images from MAUDE for a given time range. + Get the CCD temperature for a given time range. + + This is a light wrapper around cheta.fetch_sci.Msid("aacccdpt", start, stop). Parameters ---------- @@ -73,33 +53,26 @@ def get_maude_images(start: CxoTimeLike, stop: CxoTimeLike, fetch_limit_hours=10 Start time of the time range. stop : CxoTimeLike Stop time of the time range. - fetch_limit_hours : int - Maximum number of hours to fetch from MAUDE in one go. + source : str, optional + Source of the CCD temperature data. If 'maude', override the fetch_sci data source + to 'maude allow_subset=False'. Else, use the cheta default. Returns ------- - dict - dict with a astropy.table of ACA image data for each slot. + t_ccd : np.array + CCD temperature values. + t_ccd_times : np.array + Times corresponding to the CCD temperature values. """ - tstart = CxoTime(start).secs - tstop = CxoTime(stop).secs - - if tstop - tstart > fetch_limit_hours * 60 * 60: - raise ValueError("Time range too large for maude fetch, limit is 100 hours") - - slot_chunks = [] - # Break maude fetches into max 3 hour chunks required by maude_decom fetch - for mstart in np.arange(tstart, tstop, 60 * 60 * 3): - mstop = np.min([tstop, mstart + 60 * 60 * 3]) - imgs = maude_decom.get_aca_images(mstart, mstop) - slot_chunks.append(imgs) - slot_chunks = vstack(slot_chunks) - - slot_data = {} - for slot in range(8): - slot_data[slot] = slot_chunks[slot_chunks["IMGNUM"] == slot] - - return slot_data + if source == "maude": + # Override the data_source to be explicit about maude source. + with fetch_sci.data_source("maude allow_subset=False"): + dat = fetch_sci.Msid("aacccdpt", start, stop) + else: + dat = fetch_sci.Msid("aacccdpt", start, stop) + t_ccd = dat.vals + t_ccd_times = dat.times + return t_ccd, t_ccd_times def get_dcsub_aca_images( @@ -111,6 +84,7 @@ def get_dcsub_aca_images( dc_img=None, dc_tccd=None, source="maude", + full_table=False, ): """ Get dark current subtracted ACA images for a given time range. @@ -122,7 +96,8 @@ def get_dcsub_aca_images( stop : str, optional Stop time of the time range. Default is None. aca_images : dict, optional - Dictionary of ACA images. Should be keyed by slot. Default is None. + Dictionary of ACA image tables. Should be keyed by slot. Default is None. + This is intended for testing. t_ccd : array-like, optional Array of CCD temperatures. Default is None. t_ccd_times : array-like, optional @@ -133,11 +108,15 @@ def get_dcsub_aca_images( Dark current CCD temperature. Default is None. source : str, optional Source of the ACA images. Can be 'maude' or 'mica'. Default is 'maude'. + full_table : bool, optional + If True, return the full table of ACA images with dark and bgsub columns. Returns ------- dict - Dictionary of dark current subtracted ACA images for each slot. + If full_table - dictionary of aca slot data keyed by slot including raw and + dark current subtracted images. + If not full_table - dictionary of dark current subtracted ACA images keyed by slot. Raises ------ @@ -150,24 +129,10 @@ def get_dcsub_aca_images( assert dc_img.shape == (1024, 1024) if dc_img is None: - dc = get_dark_cal_props(start, "nearest", include_image=True, aca_image=False) - dc_img = dc["image"] - if "t_ccd" in dc: - dc_tccd = dc["t_ccd"] - elif "ccd_temp" in dc: - dc_tccd = dc["ccd_temp"] - else: - raise ValueError("Dark cal must have t_ccd or ccd_temp") + dc_img, dc_tccd = get_dark_data(start) if t_ccd is None: - if source == "maude": - with fetch_sci.data_source("maude allow_subset=False"): - dat = fetch_sci.Msid("aacccdpt", start, stop) - else: - with fetch_sci.data_source("cxc"): - dat = fetch_sci.Msid("aacccdpt", start, stop) - t_ccd = dat.vals - t_ccd_times = dat.times + t_ccd, t_ccd_times = get_tccd_data(start, stop, source=source) if aca_images is not None: # If supplied a dictionary of aca images, use that. Should be keyed by slot. @@ -182,7 +147,7 @@ def get_dcsub_aca_images( # require that IMGROW0_8X8 and IMGCOL0_8X8 are masked columns assert hasattr(aca_images[slot]["IMGROW0_8X8"], "mask") assert hasattr(aca_images[slot]["IMGCOL0_8X8"], "mask") - imgs = aca_images + imgs = aca_images.copy() elif source == "maude": imgs = get_maude_images(start, stop) elif source == "mica": @@ -203,14 +168,22 @@ def get_dcsub_aca_images( img_col = "IMG" imgs_bgsub[slot] = imgs[slot][img_col] - dark_imgs[slot] imgs_bgsub[slot].clip(0, None) + if full_table: + imgs[slot]["IMGBGSUB"] = imgs_bgsub[slot] + imgs[slot]["IMGDARK"] = dark_imgs[slot] + imgs[slot]["AACCCDPT"] = np.interp(imgs[slot]["TIME"], t_ccd_times, t_ccd) - return imgs_bgsub + if not full_table: + return imgs_bgsub + else: + return imgs def get_dark_current_imgs(img_table, dc_img, dc_tccd, t_ccd, t_ccd_times): """ Get the scaled dark current values for a table of ACA images. + Parameters ---------- img_table : astry.table.Table @@ -220,7 +193,7 @@ def get_dark_current_imgs(img_table, dc_img, dc_tccd, t_ccd, t_ccd_times): dc_tccd : float The reference temperature of the dark calibration image. t_ccd : array - The temperature values for the time range of the img_table. + The temperature values covering the time range of the img_table. t_ccd_times : type The corresponding times for the temperature values. @@ -244,10 +217,13 @@ def get_dark_current_imgs(img_table, dc_img, dc_tccd, t_ccd, t_ccd_times): dark = np.zeros((len(dark_raw), 8, 8), dtype=np.float64) # Scale the dark current at each dark cal 8x8 image to the ccd temperature and e-/s - for i, (eight, t_ccd_i) in enumerate(zip(dark_raw, dt_ccd, strict=True)): + for i, (eight, t_ccd_i, img) in enumerate( + zip(dark_raw, dt_ccd, img_table, strict=True) + ): img_sc = dark_temp_scale_img(eight, dc_tccd, t_ccd_i) - # Note that this just uses standard integration time of 1.696 sec - dark[i] = img_sc * 1.696 / 5 + # Default to using 1.696 integ time if not present in the image table + integ = img["INTEG"] if "INTEG" in img.colnames else 1.696 + dark[i] = img_sc * integ / 5 return dark @@ -273,6 +249,7 @@ def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): The dark backgrounds for the ACA images. """ + # Borrowed from the agasc code @numba.jit(nopython=True) def staggered_aca_slice(array_in, array_out, row, col): for i in np.arange(len(row)): @@ -284,3 +261,165 @@ def staggered_aca_slice(array_in, array_out, row, col): raw_dark_img.astype(float), dark_img, 512 + imgrow0, 512 + imgcol0 ) return dark_img + + +@lru_cache(maxsize=2) +def get_mica_images(start: CxoTimeLike, stop: CxoTimeLike, cols=None): + """ + Get ACA images from mica for a given time range. + + These mica images are centered in the 8x8 images via the centered_8x8=True option + to aca_l0.get_slot_data to match the maude images. + + Parameters + ---------- + start : CxoTimeLike + Start time of the time range. + stop : CxoTimeLike + Stop time of the time range. + + Returns + ------- + dict + dict with a astropy.table of ACA image data for each slot. + + """ + if cols is None: + cols = [ + "TIME", + "QUALITY", + "IMGFUNC1", + "IMGRAW", + "IMGSIZE", + "IMGROW0", + "IMGCOL0", + "BGDAVG", + "INTEG", + ] + else: + for col in ["TIME", "IMGRAW", "IMGROW0", "IMGCOL0"]: + if col not in cols: + cols.append(col) + + slot_data = {} + for slot in range(8): + slot_data[slot] = aca_l0.get_slot_data( + start, + stop, + imgsize=[4, 6, 8], + slot=slot, + columns=cols, + centered_8x8=True, + ) + slot_data[slot] = Table(slot_data[slot]) + + # Rename and rejigger colums to match maude + IMGROW0_8X8 = slot_data[slot]["IMGROW0"].copy() + IMGCOL0_8X8 = slot_data[slot]["IMGCOL0"].copy() + IMGROW0_8X8[slot_data[slot]["IMGSIZE"] == 6] -= 1 + IMGCOL0_8X8[slot_data[slot]["IMGSIZE"] == 6] -= 1 + IMGROW0_8X8[slot_data[slot]["IMGSIZE"] == 4] -= 2 + IMGCOL0_8X8[slot_data[slot]["IMGSIZE"] == 4] -= 2 + slot_data[slot]["IMGROW0_8X8"] = IMGROW0_8X8 + slot_data[slot]["IMGCOL0_8X8"] = IMGCOL0_8X8 + slot_data[slot].rename_column("IMGRAW", "IMG") + + return slot_data + + +@lru_cache(maxsize=2) +def get_maude_images(start: CxoTimeLike, stop: CxoTimeLike, fetch_limit_hours=100): + """ + Get ACA images from MAUDE for a given time range. + + Parameters + ---------- + start : CxoTimeLike + Start time of the time range. + stop : CxoTimeLike + Stop time of the time range. + fetch_limit_hours : int + Maximum number of hours to fetch from MAUDE in one go. + + Returns + ------- + dict + dict with a astropy.table of ACA image data for each slot. + """ + tstart = CxoTime(start).secs + tstop = CxoTime(stop).secs + + if tstop - tstart > fetch_limit_hours * 60 * 60: + raise ValueError("Time range too large for maude fetch, limit is 100 hours") + + slot_chunks = [] + # Break maude fetches into max 3 hour chunks required by maude_decom fetch + for mstart in np.arange(tstart, tstop, 60 * 60 * 3): + mstop = np.min([tstop, mstart + 60 * 60 * 3]) + imgs = maude_decom.get_aca_images(mstart, mstop) + slot_chunks.append(imgs) + slot_chunks = vstack(slot_chunks) + + slot_data = {} + for slot in range(8): + slot_data[slot] = slot_chunks[slot_chunks["IMGNUM"] == slot] + + return slot_data + + +def get_aca_images( + start=None, + stop=None, + aca_images=None, + t_ccd=None, + t_ccd_times=None, + dc_img=None, + dc_tccd=None, + source="maude", +): + """ + Get aca images for a given time range. + + + Parameters + ---------- + start : str, optional + Start time of the time range. Default is None. + stop : str, optional + Stop time of the time range. Default is None. + aca_images : dict, optional + Dictionary of ACA images. Should be keyed by slot. Default is None. + t_ccd : array-like, optional + Array of CCD temperatures. Default is None. + t_ccd_times : array-like, optional + Array of CCD temperature times. Default is None. + dc_img : array-like, optional + Dark current image. Default is None. + dc_tccd : float, optional + Dark current CCD temperature. Default is None. + source : str, optional + Source of the ACA images. Can be 'maude' or 'mica'. Default is 'maude'. + + Returns + ------- + dict + Dictionary of dark current subtracted ACA images for each slot. + + Raises + ------ + ValueError + If the dark calibration does not have t_ccd or ccd_temp. + If aca_image_table is not defined or source is not 'maude' or 'mica'. + + """ + return get_dcsub_aca_images( + start=start, + stop=stop, + aca_images=aca_images, + t_ccd=t_ccd, + t_ccd_times=t_ccd_times, + dc_img=dc_img, + dc_tccd=dc_tccd, + source=source, + full_table=True, + ) diff --git a/chandra_aca/tests/test_dark_subtract.py b/chandra_aca/tests/test_dark_subtract.py index 3470f442..2be24821 100644 --- a/chandra_aca/tests/test_dark_subtract.py +++ b/chandra_aca/tests/test_dark_subtract.py @@ -6,6 +6,9 @@ from astropy.table import Table from mica.archive.aca_dark import get_dark_cal_props from numpy import ma +from cxotime import CxoTime +import cheta.fetch + from chandra_aca.dark_subtract import ( get_dark_backgrounds, @@ -13,6 +16,9 @@ get_dcsub_aca_images, get_maude_images, get_mica_images, + get_aca_images, + get_dark_data, + get_tccd_data, ) HAS_ACA0_ARCHIVE = (Path(mica.common.MICA_ARCHIVE) / "aca0").exists() @@ -161,6 +167,54 @@ def test_dcsub_aca_images(mock_dc, mock_imgs, dc_imgs_dn): assert np.allclose(imgs_bgsub[0] * 5 / 1.696, exp, atol=1e-6, rtol=0) +def test_get_dark_data_recent(): + date = "2021:001:00:00:00.000" + dark_data_img, dark_data_tccd = get_dark_data(date) + assert dark_data_img.shape == (1024, 1024) + assert dark_data_tccd == -8.56 + + +def test_get_dark_data_old(): + # Confirm this works for a really old one too + date = CxoTime("2001:001:00:00:00.000").secs + dark_data_img, dark_data_tccd = get_dark_data(date) + assert dark_data_img.shape == (1024, 1024) + assert dark_data_tccd == -10.3720703125 + + +def test_get_tccd_data(): + start = "2021:001:00:00:00.000" + stop = "2021:001:00:00:30.000" + t_ccd_maude, t_ccd_times_maude = get_tccd_data(start, stop, source="maude") + assert np.isclose(t_ccd_maude[0], -6.55137778) + + # Confirm the t_ccd values are the same for maude as default fetch data source + # but only bother if cxc is in the sources. Technically this should be a skipif. + if "cxc" in cheta.fetch.data_source.sources(): + t_ccd, t_ccd_times = get_tccd_data(start, stop) + assert np.allclose(t_ccd_maude, t_ccd) + assert np.allclose(t_ccd_times_maude, t_ccd_times) + + +def test_get_aca_images(mock_dc, mock_imgs, dc_imgs_dn): + """ + Confirm the pattern of dark current images matches the reference. + """ + dat = get_aca_images( + aca_images=mock_imgs, + dc_img=mock_dc, + dc_tccd=-10, + t_ccd=np.repeat(-10, 8), + t_ccd_times=mock_imgs[0]["TIME"], + ) + for col in ["IMG", "IMGROW0_8X8", "IMGCOL0_8X8", "TIME", "IMGBGSUB"]: + assert col in dat[0].colnames + assert np.allclose(dat[0][col], mock_imgs[0][col], atol=1e-6, rtol=0) + assert dat[0]["IMGBGSUB"].shape == (8, 8, 8) + exp = 16 - dc_imgs_dn + assert np.allclose(dat[0]["IMGBGSUB"] * 5 / 1.696, exp, atol=1e-6, rtol=0) + + def test_get_dark_images(mock_dc, mock_imgs, dc_imgs_dn): """ Confirm the pattern of dark current images matches the reference. From ee6086f1ee994c4ea0c8848573421a1fc5861e40 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Tue, 17 Sep 2024 12:03:41 -0400 Subject: [PATCH 03/28] Add retry --- chandra_aca/dark_subtract.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index 025b0508..b715e193 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -2,11 +2,13 @@ import numba import numpy as np +import requests.exceptions from astropy.table import Table, vstack from cheta import fetch_sci from cxotime import CxoTime, CxoTimeLike from mica.archive import aca_l0 from mica.archive.aca_dark import get_dark_cal_props +from ska_helpers import retry from chandra_aca import maude_decom from chandra_aca.dark_model import dark_temp_scale_img @@ -41,6 +43,7 @@ def get_dark_data(time: CxoTimeLike): return dc_img, dc_tccd +@retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): """ Get the CCD temperature for a given time range. @@ -75,6 +78,7 @@ def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): return t_ccd, t_ccd_times +@retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) def get_dcsub_aca_images( start=None, stop=None, @@ -183,7 +187,6 @@ def get_dark_current_imgs(img_table, dc_img, dc_tccd, t_ccd, t_ccd_times): """ Get the scaled dark current values for a table of ACA images. - Parameters ---------- img_table : astry.table.Table @@ -327,6 +330,7 @@ def get_mica_images(start: CxoTimeLike, stop: CxoTimeLike, cols=None): return slot_data +@retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) @lru_cache(maxsize=2) def get_maude_images(start: CxoTimeLike, stop: CxoTimeLike, fetch_limit_hours=100): """ @@ -356,6 +360,8 @@ def get_maude_images(start: CxoTimeLike, stop: CxoTimeLike, fetch_limit_hours=10 # Break maude fetches into max 3 hour chunks required by maude_decom fetch for mstart in np.arange(tstart, tstop, 60 * 60 * 3): mstop = np.min([tstop, mstart + 60 * 60 * 3]) + # Add a retry here if necessary + imgs = maude_decom.get_aca_images(mstart, mstop) slot_chunks.append(imgs) slot_chunks = vstack(slot_chunks) @@ -367,6 +373,7 @@ def get_maude_images(start: CxoTimeLike, stop: CxoTimeLike, fetch_limit_hours=10 return slot_data +@retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) def get_aca_images( start=None, stop=None, @@ -380,7 +387,6 @@ def get_aca_images( """ Get aca images for a given time range. - Parameters ---------- start : str, optional From 7a3290ab00766855ba96a54afa2ee8e550ab470d Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Fri, 20 Sep 2024 16:16:45 -0400 Subject: [PATCH 04/28] First round simplifications --- chandra_aca/dark_subtract.py | 257 +++++------------------- chandra_aca/tests/test_dark_subtract.py | 120 ++++------- 2 files changed, 89 insertions(+), 288 deletions(-) diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index b715e193..d1395858 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -1,12 +1,9 @@ -from functools import lru_cache - import numba import numpy as np import requests.exceptions -from astropy.table import Table, vstack +from astropy.table import Table from cheta import fetch_sci -from cxotime import CxoTime, CxoTimeLike -from mica.archive import aca_l0 +from cxotime import CxoTimeLike from mica.archive.aca_dark import get_dark_cal_props from ska_helpers import retry @@ -14,35 +11,6 @@ from chandra_aca.dark_model import dark_temp_scale_img -def get_dark_data(time: CxoTimeLike): - """ - Retrieve the dark calibration image and temperature for a given time. - - This is a light wrapper around mica.archive.aca_dark.get_dark_cal_props. - - Parameters - ---------- - time : CxoTimeLike - Reference time for the dark calibration image. - - Returns - ------- - dc_img : np.array - Dark calibration image. - dc_tccd : float - Dark calibration CCD temperature. - """ - dc = get_dark_cal_props(time, "nearest", include_image=True, aca_image=False) - dc_img = dc["image"] - if "t_ccd" in dc: - dc_tccd = dc["t_ccd"] - elif "ccd_temp" in dc: - dc_tccd = dc["ccd_temp"] - else: - raise ValueError("Dark cal must have t_ccd or ccd_temp") - return dc_img, dc_tccd - - @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): """ @@ -79,10 +47,10 @@ def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) -def get_dcsub_aca_images( +def _get_dcsub_aca_images( start=None, stop=None, - aca_images=None, + aca_image_table=None, t_ccd=None, t_ccd_times=None, dc_img=None, @@ -99,7 +67,7 @@ def get_dcsub_aca_images( Start time of the time range. Default is None. stop : str, optional Stop time of the time range. Default is None. - aca_images : dict, optional + aca_image_table : dict, optional Dictionary of ACA image tables. Should be keyed by slot. Default is None. This is intended for testing. t_ccd : array-like, optional @@ -129,58 +97,65 @@ def get_dcsub_aca_images( If aca_image_table is not defined or source is not 'maude' or 'mica'. """ - if dc_img is not None: - assert dc_img.shape == (1024, 1024) - if dc_img is None: - dc_img, dc_tccd = get_dark_data(start) - - if t_ccd is None: - t_ccd, t_ccd_times = get_tccd_data(start, stop, source=source) - - if aca_images is not None: - # If supplied a dictionary of aca images, use that. Should be keyed by slot. - # This is for testing. - assert isinstance(aca_images, dict) - # For each slot, require aca_image have columns IMG, IMGROW0_8X8, IMGCOL0_8X8, TIME - for slot in aca_images: - assert "IMG" in aca_images[slot].colnames - assert "IMGROW0_8X8" in aca_images[slot].colnames - assert "IMGCOL0_8X8" in aca_images[slot].colnames - assert "TIME" in aca_images[slot].colnames - # require that IMGROW0_8X8 and IMGCOL0_8X8 are masked columns - assert hasattr(aca_images[slot]["IMGROW0_8X8"], "mask") - assert hasattr(aca_images[slot]["IMGCOL0_8X8"], "mask") - imgs = aca_images.copy() + if aca_image_table is not None: + # If supplied an image table, use that, but confirm it has the right pieces. + assert type(aca_image_table) is Table + # Require aca_image_table have columns IMG, IMGROW0_8X8, IMGCOL0_8X8, TIME + assert "IMG" in aca_image_table.colnames + assert "IMGROW0_8X8" in aca_image_table.colnames + assert "IMGCOL0_8X8" in aca_image_table.colnames + assert "TIME" in aca_image_table.colnames + # require that IMGROW0_8X8 and IMGCOL0_8X8 are masked columns + assert hasattr(aca_image_table["IMGROW0_8X8"], "mask") + assert hasattr(aca_image_table["IMGCOL0_8X8"], "mask") + # Make a copy so we don't modify the input table + img_table = aca_image_table.copy() elif source == "maude": - imgs = get_maude_images(start, stop) + if start is None or stop is None: + raise ValueError("start and stop must be defined for maude source") + img_table = maude_decom.get_aca_images(start, stop) elif source == "mica": - imgs = get_mica_images(start, stop) + raise NotImplementedError("mica source not yet implemented") else: raise ValueError( "aca_image_table must be defined or source must be maude or mica" ) - imgs_bgsub = {} - dark_imgs = {} - for slot in range(8): - if slot not in imgs: - continue - dark_imgs[slot] = get_dark_current_imgs( - imgs[slot], dc_img, dc_tccd, t_ccd, t_ccd_times + if dc_img is not None: + assert dc_img.shape == (1024, 1024) + + if start is None: + start = img_table["TIME"].min() + if stop is None: + stop = img_table["TIME"].max() + + if dc_img is None: + dark_data = get_dark_cal_props( + start, select="nearest", include_image=True, aca_image=False ) - img_col = "IMG" - imgs_bgsub[slot] = imgs[slot][img_col] - dark_imgs[slot] - imgs_bgsub[slot].clip(0, None) - if full_table: - imgs[slot]["IMGBGSUB"] = imgs_bgsub[slot] - imgs[slot]["IMGDARK"] = dark_imgs[slot] - imgs[slot]["AACCCDPT"] = np.interp(imgs[slot]["TIME"], t_ccd_times, t_ccd) + dc_img = dark_data["image"] + dc_tccd = dark_data["ccd_temp"] + + if t_ccd is None: + t_ccd, t_ccd_times = get_tccd_data(start, stop, source=source) + + imgs_bgsub = {} + imgs_dark = {} + + imgs_dark = get_dark_current_imgs(img_table, dc_img, dc_tccd, t_ccd, t_ccd_times) + imgs_bgsub = img_table["IMG"] - imgs_dark + imgs_bgsub.clip(0, None) + + if full_table: + img_table["IMGBGSUB"] = imgs_bgsub + img_table["IMGDARK"] = imgs_dark + img_table["AACCCDPT"] = np.interp(img_table["TIME"], t_ccd_times, t_ccd) if not full_table: return imgs_bgsub else: - return imgs + return img_table def get_dark_current_imgs(img_table, dc_img, dc_tccd, t_ccd, t_ccd_times): @@ -266,122 +241,11 @@ def staggered_aca_slice(array_in, array_out, row, col): return dark_img -@lru_cache(maxsize=2) -def get_mica_images(start: CxoTimeLike, stop: CxoTimeLike, cols=None): - """ - Get ACA images from mica for a given time range. - - These mica images are centered in the 8x8 images via the centered_8x8=True option - to aca_l0.get_slot_data to match the maude images. - - Parameters - ---------- - start : CxoTimeLike - Start time of the time range. - stop : CxoTimeLike - Stop time of the time range. - - Returns - ------- - dict - dict with a astropy.table of ACA image data for each slot. - - """ - if cols is None: - cols = [ - "TIME", - "QUALITY", - "IMGFUNC1", - "IMGRAW", - "IMGSIZE", - "IMGROW0", - "IMGCOL0", - "BGDAVG", - "INTEG", - ] - else: - for col in ["TIME", "IMGRAW", "IMGROW0", "IMGCOL0"]: - if col not in cols: - cols.append(col) - - slot_data = {} - for slot in range(8): - slot_data[slot] = aca_l0.get_slot_data( - start, - stop, - imgsize=[4, 6, 8], - slot=slot, - columns=cols, - centered_8x8=True, - ) - slot_data[slot] = Table(slot_data[slot]) - - # Rename and rejigger colums to match maude - IMGROW0_8X8 = slot_data[slot]["IMGROW0"].copy() - IMGCOL0_8X8 = slot_data[slot]["IMGCOL0"].copy() - IMGROW0_8X8[slot_data[slot]["IMGSIZE"] == 6] -= 1 - IMGCOL0_8X8[slot_data[slot]["IMGSIZE"] == 6] -= 1 - IMGROW0_8X8[slot_data[slot]["IMGSIZE"] == 4] -= 2 - IMGCOL0_8X8[slot_data[slot]["IMGSIZE"] == 4] -= 2 - slot_data[slot]["IMGROW0_8X8"] = IMGROW0_8X8 - slot_data[slot]["IMGCOL0_8X8"] = IMGCOL0_8X8 - slot_data[slot].rename_column("IMGRAW", "IMG") - - return slot_data - - -@retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) -@lru_cache(maxsize=2) -def get_maude_images(start: CxoTimeLike, stop: CxoTimeLike, fetch_limit_hours=100): - """ - Get ACA images from MAUDE for a given time range. - - Parameters - ---------- - start : CxoTimeLike - Start time of the time range. - stop : CxoTimeLike - Stop time of the time range. - fetch_limit_hours : int - Maximum number of hours to fetch from MAUDE in one go. - - Returns - ------- - dict - dict with a astropy.table of ACA image data for each slot. - """ - tstart = CxoTime(start).secs - tstop = CxoTime(stop).secs - - if tstop - tstart > fetch_limit_hours * 60 * 60: - raise ValueError("Time range too large for maude fetch, limit is 100 hours") - - slot_chunks = [] - # Break maude fetches into max 3 hour chunks required by maude_decom fetch - for mstart in np.arange(tstart, tstop, 60 * 60 * 3): - mstop = np.min([tstop, mstart + 60 * 60 * 3]) - # Add a retry here if necessary - - imgs = maude_decom.get_aca_images(mstart, mstop) - slot_chunks.append(imgs) - slot_chunks = vstack(slot_chunks) - - slot_data = {} - for slot in range(8): - slot_data[slot] = slot_chunks[slot_chunks["IMGNUM"] == slot] - - return slot_data - - @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) def get_aca_images( start=None, stop=None, - aca_images=None, - t_ccd=None, - t_ccd_times=None, - dc_img=None, - dc_tccd=None, + aca_image_table=None, source="maude", ): """ @@ -395,14 +259,6 @@ def get_aca_images( Stop time of the time range. Default is None. aca_images : dict, optional Dictionary of ACA images. Should be keyed by slot. Default is None. - t_ccd : array-like, optional - Array of CCD temperatures. Default is None. - t_ccd_times : array-like, optional - Array of CCD temperature times. Default is None. - dc_img : array-like, optional - Dark current image. Default is None. - dc_tccd : float, optional - Dark current CCD temperature. Default is None. source : str, optional Source of the ACA images. Can be 'maude' or 'mica'. Default is 'maude'. @@ -414,18 +270,13 @@ def get_aca_images( Raises ------ ValueError - If the dark calibration does not have t_ccd or ccd_temp. If aca_image_table is not defined or source is not 'maude' or 'mica'. """ - return get_dcsub_aca_images( + return _get_dcsub_aca_images( start=start, stop=stop, - aca_images=aca_images, - t_ccd=t_ccd, - t_ccd_times=t_ccd_times, - dc_img=dc_img, - dc_tccd=dc_tccd, + aca_image_table=aca_image_table, source=source, full_table=True, ) diff --git a/chandra_aca/tests/test_dark_subtract.py b/chandra_aca/tests/test_dark_subtract.py index 2be24821..ef5d57d7 100644 --- a/chandra_aca/tests/test_dark_subtract.py +++ b/chandra_aca/tests/test_dark_subtract.py @@ -1,23 +1,18 @@ from pathlib import Path +import cheta.fetch import mica.common import numpy as np import pytest from astropy.table import Table from mica.archive.aca_dark import get_dark_cal_props from numpy import ma -from cxotime import CxoTime -import cheta.fetch - +from chandra_aca import maude_decom from chandra_aca.dark_subtract import ( + _get_dcsub_aca_images, get_dark_backgrounds, get_dark_current_imgs, - get_dcsub_aca_images, - get_maude_images, - get_mica_images, - get_aca_images, - get_dark_data, get_tccd_data, ) @@ -40,20 +35,19 @@ def mock_dc(): # Make a test fixture for the mock 8x8 data @pytest.fixture -def mock_imgs(): +def mock_img_table(): """ Make a mock 8x8 image table. The intent is that these image coordinates should overlap with the mock dark current. """ - imgs = {} img_table = Table() img_table["TIME"] = np.arange(8) img_table["IMGROW0_8X8"] = ma.arange(8) - 512 img_table["IMGCOL0_8X8"] = ma.arange(8) - 512 img_table["IMG"] = np.ones((8, 8, 8)) * 16 * 1.696 / 5 - imgs[0] = img_table - return imgs + img_table["IMGNUM"] = 0 + return img_table # Test fixture array of dark current images in DN @@ -149,37 +143,22 @@ def dc_imgs_dn(): return dc_imgs_dn -def test_dcsub_aca_images(mock_dc, mock_imgs, dc_imgs_dn): +def test_dcsub_aca_images(mock_dc, mock_img_table, dc_imgs_dn): """ Confirm that the pattern of background subtraction is correct. """ - imgs_bgsub = get_dcsub_aca_images( - aca_images=mock_imgs, + imgs_bgsub = _get_dcsub_aca_images( + aca_image_table=mock_img_table, dc_img=mock_dc, dc_tccd=-10, t_ccd=np.repeat(-10, 8), - t_ccd_times=mock_imgs[0]["TIME"], + t_ccd_times=mock_img_table["TIME"], ) - assert imgs_bgsub[0].shape == (8, 8, 8) + assert imgs_bgsub.shape == (8, 8, 8) # Note that the mock unsubtracted data is originally 16 * 1.696 / 5 exp = 16 - dc_imgs_dn - assert np.allclose(imgs_bgsub[0] * 5 / 1.696, exp, atol=1e-6, rtol=0) - - -def test_get_dark_data_recent(): - date = "2021:001:00:00:00.000" - dark_data_img, dark_data_tccd = get_dark_data(date) - assert dark_data_img.shape == (1024, 1024) - assert dark_data_tccd == -8.56 - - -def test_get_dark_data_old(): - # Confirm this works for a really old one too - date = CxoTime("2001:001:00:00:00.000").secs - dark_data_img, dark_data_tccd = get_dark_data(date) - assert dark_data_img.shape == (1024, 1024) - assert dark_data_tccd == -10.3720703125 + assert np.allclose(imgs_bgsub * 5 / 1.696, exp, atol=1e-6, rtol=0) def test_get_tccd_data(): @@ -196,90 +175,76 @@ def test_get_tccd_data(): assert np.allclose(t_ccd_times_maude, t_ccd_times) -def test_get_aca_images(mock_dc, mock_imgs, dc_imgs_dn): +def test_get_aca_images(mock_dc, mock_img_table, dc_imgs_dn): """ Confirm the pattern of dark current images matches the reference. """ - dat = get_aca_images( - aca_images=mock_imgs, + dat = _get_dcsub_aca_images( + aca_image_table=mock_img_table, dc_img=mock_dc, dc_tccd=-10, t_ccd=np.repeat(-10, 8), - t_ccd_times=mock_imgs[0]["TIME"], + t_ccd_times=mock_img_table["TIME"], + full_table=True, ) - for col in ["IMG", "IMGROW0_8X8", "IMGCOL0_8X8", "TIME", "IMGBGSUB"]: - assert col in dat[0].colnames - assert np.allclose(dat[0][col], mock_imgs[0][col], atol=1e-6, rtol=0) - assert dat[0]["IMGBGSUB"].shape == (8, 8, 8) + for col in ["IMG", "IMGROW0_8X8", "IMGCOL0_8X8", "TIME"]: + assert col in dat.colnames + assert np.allclose(dat[col], mock_img_table[col], atol=1e-6, rtol=0) + assert dat["IMGBGSUB"].shape == (8, 8, 8) exp = 16 - dc_imgs_dn - assert np.allclose(dat[0]["IMGBGSUB"] * 5 / 1.696, exp, atol=1e-6, rtol=0) + assert np.allclose(dat["IMGBGSUB"] * 5 / 1.696, exp, atol=1e-6, rtol=0) -def test_get_dark_images(mock_dc, mock_imgs, dc_imgs_dn): +def test_get_dark_images(mock_dc, mock_img_table, dc_imgs_dn): """ Confirm the pattern of dark current images matches the reference. The temperature scaling is set here to be a no-op. """ dc_imgs = get_dark_current_imgs( - mock_imgs[0], + img_table=mock_img_table, dc_img=mock_dc, dc_tccd=-10, t_ccd=np.repeat(-10, 8), - t_ccd_times=mock_imgs[0]["TIME"], + t_ccd_times=mock_img_table["TIME"], ) assert dc_imgs.shape == (8, 8, 8) dc_imgs_raw = dc_imgs * 5 / 1.696 assert np.allclose(dc_imgs_raw, dc_imgs_dn, atol=1e-6, rtol=0) -def test_get_dark_backgrounds(mock_dc, mock_imgs, dc_imgs_dn): +def test_get_dark_backgrounds(mock_dc, mock_img_table, dc_imgs_dn): """ Confirm the pattern of dark current matches the reference. """ dc_imgs_raw = get_dark_backgrounds( mock_dc, - mock_imgs[0]["IMGROW0_8X8"].data.filled(1025), - mock_imgs[0]["IMGCOL0_8X8"].data.filled(1025), + mock_img_table["IMGROW0_8X8"].data.filled(1025), + mock_img_table["IMGCOL0_8X8"].data.filled(1025), ) assert dc_imgs_raw.shape == (8, 8, 8) assert np.allclose(dc_imgs_raw, dc_imgs_dn, atol=1e-6, rtol=0) -@pytest.mark.skipif(not HAS_ACA0_ARCHIVE, reason="No ACA0 archive") -def test_get_images(): - """ - Confirm mica and maude images sources agree for a small time range. - """ - tstart = 746484519.429 - imgs_maude = get_maude_images(tstart, tstart + 30) - imgs_mica = get_mica_images(tstart, tstart + 30) - for slot in range(8): - assert len(imgs_maude[slot]) == len(imgs_mica[slot]) - for col in ["IMG", "IMGROW0_8X8", "IMGCOL0_8X8", "TIME"]: - assert np.allclose( - imgs_maude[slot][col], imgs_mica[slot][col], atol=0.1, rtol=0 - ) - - def test_dc_consistent(): """ Confirm extracted dark current is reasonably consistent with the images edges of a 6x6. """ tstart = 746484519.429 - imgs_maude = get_maude_images(tstart, tstart + 30) + img_table = maude_decom.get_aca_images(tstart, tstart + 30) # This slot and this time show an overall decent correlation of dark current # and edge pixel values in the 6x6 image - but don't make a great test. slot = 4 - img = imgs_maude[slot]["IMG"][0] + img_table_slot = img_table[img_table["IMGNUM"] == slot] + img = img_table_slot[0]["IMG"] dc = get_dark_cal_props(tstart, "nearest", include_image=True, aca_image=False) dc_imgs = get_dark_current_imgs( - imgs_maude[slot], + img_table_slot, dc_img=dc["image"], dc_tccd=dc["t_ccd"], - t_ccd=np.repeat(-9.9, len(imgs_maude[slot])), - t_ccd_times=imgs_maude[slot]["TIME"], + t_ccd=np.repeat(-9.9, len(img_table_slot)), + t_ccd_times=img_table_slot["TIME"], ) edge_mask_6x6 = np.zeros((8, 8), dtype=bool) edge_mask_6x6[1, 2:6] = True @@ -302,7 +267,7 @@ def test_dcsub_aca_images_maude(): # This one is just a regression test tstart = 471139130.93277466 # dwell tstart for obsid 15600 - random obsid tstop = tstart + 20 - imgs_bgsub = get_dcsub_aca_images(tstart, tstop, source="maude") + imgs_bgsub = _get_dcsub_aca_images(tstart, tstop, source="maude") exp0 = np.array( [ @@ -332,18 +297,3 @@ def test_dcsub_aca_images_maude(): ] ) np.allclose(exp4, imgs_bgsub[4][0].filled(0), atol=1, rtol=0) - - -@pytest.mark.skipif(not HAS_ACA0_ARCHIVE, reason="No ACA0 archive") -def test_dcsub_aca_images_mica_maude(): - """ - Confirm background subtracted mica/maude images match for small time range. - """ - tstart = 471139130.93277466 # dwell tstart for obsid 15600 - random obsid - tstop = tstart + 30 - imgs_bgsub = get_dcsub_aca_images(tstart, tstop, source="mica") - - imgs_bgsub_maude = get_dcsub_aca_images(tstart, tstop, source="maude") - - for slot in range(8): - assert np.allclose(imgs_bgsub[slot], imgs_bgsub_maude[slot], atol=1e-3, rtol=0) From 820211a2b5013ea83195d46a6dbfeba8fbb841db Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Tue, 24 Sep 2024 17:20:16 -0400 Subject: [PATCH 05/28] Second round simplifications --- chandra_aca/dark_subtract.py | 174 ++++++++++++------------ chandra_aca/tests/test_dark_subtract.py | 31 +++-- 2 files changed, 102 insertions(+), 103 deletions(-) diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index d1395858..b6663180 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -4,6 +4,7 @@ from astropy.table import Table from cheta import fetch_sci from cxotime import CxoTimeLike +from mica.archive import aca_l0 from mica.archive.aca_dark import get_dark_cal_props from ska_helpers import retry @@ -16,7 +17,8 @@ def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): """ Get the CCD temperature for a given time range. - This is a light wrapper around cheta.fetch_sci.Msid("aacccdpt", start, stop). + This is a light wrapper around cheta.fetch_sci.Msid("aacccdpt", start, stop) + to handle an option to use the maude data source in an explicit way. Parameters ---------- @@ -30,7 +32,7 @@ def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): Returns ------- - t_ccd : np.array + t_ccd_vals : np.array CCD temperature values. t_ccd_times : np.array Times corresponding to the CCD temperature values. @@ -41,60 +43,52 @@ def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): dat = fetch_sci.Msid("aacccdpt", start, stop) else: dat = fetch_sci.Msid("aacccdpt", start, stop) - t_ccd = dat.vals - t_ccd_times = dat.times - return t_ccd, t_ccd_times + return dat.vals, dat.times @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) def _get_dcsub_aca_images( - start=None, - stop=None, + start: CxoTimeLike = None, + stop: CxoTimeLike = None, aca_image_table=None, - t_ccd=None, + t_ccd_vals=None, t_ccd_times=None, - dc_img=None, - dc_tccd=None, + img_dark=None, + tccd_dark=None, source="maude", - full_table=False, ): """ Get dark current subtracted ACA images for a given time range. + This private method has expanded options for unit testing. + Users should use get_dcsub_aca_images instead. + Parameters ---------- - start : str, optional + start : CxoTimeLike, optional Start time of the time range. Default is None. - stop : str, optional + stop : CxoTimeLike, optional Stop time of the time range. Default is None. - aca_image_table : dict, optional - Dictionary of ACA image tables. Should be keyed by slot. Default is None. - This is intended for testing. - t_ccd : array-like, optional + aca_image_table : astropy table of ACA images, optional + + t_ccd_vals : array-like, optional Array of CCD temperatures. Default is None. t_ccd_times : array-like, optional Array of CCD temperature times. Default is None. - dc_img : array-like, optional + img_dark : array-like, optional Dark current image. Default is None. - dc_tccd : float, optional + tccd_dark : float, optional Dark current CCD temperature. Default is None. source : str, optional Source of the ACA images. Can be 'maude' or 'mica'. Default is 'maude'. - full_table : bool, optional - If True, return the full table of ACA images with dark and bgsub columns. Returns ------- - dict - If full_table - dictionary of aca slot data keyed by slot including raw and - dark current subtracted images. - If not full_table - dictionary of dark current subtracted ACA images keyed by slot. - - Raises - ------ - ValueError - If the dark calibration does not have t_ccd or ccd_temp. - If aca_image_table is not defined or source is not 'maude' or 'mica'. + astropy.table.Table + Table of aca image data including raw and dark current subtracted images. + Raw images are in column 'IMG', dark current subtracted images are in column 'IMGBGSUB', + dark current cut-out images are in column 'IMGDARK', and interpolated + CCD temperature values are in column 'AACCCDPT'. """ @@ -111,99 +105,104 @@ def _get_dcsub_aca_images( assert hasattr(aca_image_table["IMGCOL0_8X8"], "mask") # Make a copy so we don't modify the input table img_table = aca_image_table.copy() + + # And it doesn't make sense to also define start and stop so raise errors if so + if start is not None: + raise ValueError("Do not define start if aca_image_table is defined") + if stop is not None: + raise ValueError("Do not define stop if aca_image_table is defined") elif source == "maude": if start is None or stop is None: raise ValueError("start and stop must be defined for maude source") img_table = maude_decom.get_aca_images(start, stop) elif source == "mica": - raise NotImplementedError("mica source not yet implemented") + img_table = aca_l0.get_aca_images(start, stop) else: raise ValueError( "aca_image_table must be defined or source must be maude or mica" ) - if dc_img is not None: - assert dc_img.shape == (1024, 1024) + if start is None and aca_image_table is None: + raise ValueError("start must be defined if aca_image_table is not defined") + if stop is None and aca_image_table is None: + raise ValueError("stop must be defined if aca_image_table is not defined") if start is None: start = img_table["TIME"].min() if stop is None: stop = img_table["TIME"].max() - if dc_img is None: + if img_dark is None: dark_data = get_dark_cal_props( start, select="nearest", include_image=True, aca_image=False ) - dc_img = dark_data["image"] - dc_tccd = dark_data["ccd_temp"] + img_dark = dark_data["image"] + tccd_dark = dark_data["ccd_temp"] - if t_ccd is None: - t_ccd, t_ccd_times = get_tccd_data(start, stop, source=source) + assert img_dark.shape == (1024, 1024) - imgs_bgsub = {} - imgs_dark = {} + if t_ccd_vals is None: + t_ccd_vals, t_ccd_times = get_tccd_data(start, stop, source=source) - imgs_dark = get_dark_current_imgs(img_table, dc_img, dc_tccd, t_ccd, t_ccd_times) + imgs_dark = get_dark_current_imgs( + img_table, img_dark, tccd_dark, t_ccd_vals, t_ccd_times + ) imgs_bgsub = img_table["IMG"] - imgs_dark imgs_bgsub.clip(0, None) - if full_table: - img_table["IMGBGSUB"] = imgs_bgsub - img_table["IMGDARK"] = imgs_dark - img_table["AACCCDPT"] = np.interp(img_table["TIME"], t_ccd_times, t_ccd) + img_table["IMGBGSUB"] = imgs_bgsub + img_table["IMGDARK"] = imgs_dark + img_table["AACCCDPT"] = np.interp(img_table["TIME"], t_ccd_times, t_ccd_vals) - if not full_table: - return imgs_bgsub - else: - return img_table + return img_table -def get_dark_current_imgs(img_table, dc_img, dc_tccd, t_ccd, t_ccd_times): +def get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccd_vals, t_ccd_times): """ Get the scaled dark current values for a table of ACA images. Parameters ---------- - img_table : astry.table.Table - A table containing information about the images. - dc_img : 1024x1024 array + img_table : astropy.table.Table + A table containing the ACA images and expected metadata. + img_dark : 1024x1024 array The dark calibration image. - dc_tccd : float + tccd_dark : float The reference temperature of the dark calibration image. - t_ccd : array - The temperature values covering the time range of the img_table. - t_ccd_times : type + t_ccd_vals : array + The cheta temperature values covering the time range of the img_table. + t_ccd_times : array The corresponding times for the temperature values. Returns ------- - dark : np.array + imgs_dark : np.array An array containing the temperature scaled dark current for each ACA image in the img_table in e-/s. """ - assert dc_img.shape == (1024, 1024) - assert len(t_ccd) == len(t_ccd_times) + assert img_dark.shape == (1024, 1024) + assert len(t_ccd_vals) == len(t_ccd_times) - dark_raw = get_dark_backgrounds( - dc_img, + imgs_dark_unscaled = get_dark_backgrounds( + img_dark, img_table["IMGROW0_8X8"].data.filled(1025), img_table["IMGCOL0_8X8"].data.filled(1025), ) - dt_ccd = np.interp(img_table["TIME"], t_ccd_times, t_ccd) + dt_ccd = np.interp(img_table["TIME"], t_ccd_times, t_ccd_vals) - dark = np.zeros((len(dark_raw), 8, 8), dtype=np.float64) + imgs_dark = np.zeros((len(img_table), 8, 8), dtype=np.float64) # Scale the dark current at each dark cal 8x8 image to the ccd temperature and e-/s - for i, (eight, t_ccd_i, img) in enumerate( - zip(dark_raw, dt_ccd, img_table, strict=True) + for i, (img_dark_unscaled, t_ccd_i, img) in enumerate( + zip(imgs_dark_unscaled, dt_ccd, img_table, strict=True) ): - img_sc = dark_temp_scale_img(eight, dc_tccd, t_ccd_i) + img_scaled = dark_temp_scale_img(img_dark_unscaled, tccd_dark, t_ccd_i) # Default to using 1.696 integ time if not present in the image table integ = img["INTEG"] if "INTEG" in img.colnames else 1.696 - dark[i] = img_sc * integ / 5 + imgs_dark[i] = img_scaled * integ / 5 - return dark + return imgs_dark def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): @@ -223,8 +222,10 @@ def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): Returns ------- - dark_img : np.array - The dark backgrounds for the ACA images. + imgs_dark : np.array + The dark backgrounds for the ACA images sampled from raw_dark_img. + This will have the same length as imgrow0 and imgcol0, with shape + (len(imgrow0), size, size). These pixels have not been scaled. """ # Borrowed from the agasc code @@ -234,17 +235,17 @@ def staggered_aca_slice(array_in, array_out, row, col): if row[i] + size < 1024 and col[i] + size < 1024: array_out[i] = array_in[row[i] : row[i] + size, col[i] : col[i] + size] - dark_img = np.zeros([len(imgrow0), size, size], dtype=np.float64) + imgs_dark = np.zeros([len(imgrow0), size, size], dtype=np.float64) staggered_aca_slice( - raw_dark_img.astype(float), dark_img, 512 + imgrow0, 512 + imgcol0 + raw_dark_img.astype(float), imgs_dark, 512 + imgrow0, 512 + imgcol0 ) - return dark_img + return imgs_dark @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) -def get_aca_images( - start=None, - stop=None, +def get_dcsub_aca_images( + start: CxoTimeLike = None, + stop: CxoTimeLike = None, aca_image_table=None, source="maude", ): @@ -257,20 +258,18 @@ def get_aca_images( Start time of the time range. Default is None. stop : str, optional Stop time of the time range. Default is None. - aca_images : dict, optional - Dictionary of ACA images. Should be keyed by slot. Default is None. + aca_image_table : astropy.table.Table + Table including ACA images. source : str, optional Source of the ACA images. Can be 'maude' or 'mica'. Default is 'maude'. Returns ------- - dict - Dictionary of dark current subtracted ACA images for each slot. - - Raises - ------ - ValueError - If aca_image_table is not defined or source is not 'maude' or 'mica'. + astropy.table.Table + Table of aca image data including raw and dark current subtracted images. + Raw images are in column 'IMG', dark current subtracted images are in column 'IMGBGSUB', + dark current cut-out images are in column 'IMGDARK', and interpolated + CCD temperature values are in column 'AACCCDPT'. """ return _get_dcsub_aca_images( @@ -278,5 +277,4 @@ def get_aca_images( stop=stop, aca_image_table=aca_image_table, source=source, - full_table=True, ) diff --git a/chandra_aca/tests/test_dark_subtract.py b/chandra_aca/tests/test_dark_subtract.py index ef5d57d7..144096b1 100644 --- a/chandra_aca/tests/test_dark_subtract.py +++ b/chandra_aca/tests/test_dark_subtract.py @@ -148,13 +148,14 @@ def test_dcsub_aca_images(mock_dc, mock_img_table, dc_imgs_dn): Confirm that the pattern of background subtraction is correct. """ - imgs_bgsub = _get_dcsub_aca_images( + table_bgsub = _get_dcsub_aca_images( aca_image_table=mock_img_table, - dc_img=mock_dc, - dc_tccd=-10, - t_ccd=np.repeat(-10, 8), + img_dark=mock_dc, + tccd_dark=-10, + t_ccd_vals=np.repeat(-10, 8), t_ccd_times=mock_img_table["TIME"], ) + imgs_bgsub = table_bgsub["IMGBGSUB"] assert imgs_bgsub.shape == (8, 8, 8) # Note that the mock unsubtracted data is originally 16 * 1.696 / 5 exp = 16 - dc_imgs_dn @@ -181,11 +182,10 @@ def test_get_aca_images(mock_dc, mock_img_table, dc_imgs_dn): """ dat = _get_dcsub_aca_images( aca_image_table=mock_img_table, - dc_img=mock_dc, - dc_tccd=-10, - t_ccd=np.repeat(-10, 8), + img_dark=mock_dc, + tccd_dark=-10, + t_ccd_vals=np.repeat(-10, 8), t_ccd_times=mock_img_table["TIME"], - full_table=True, ) for col in ["IMG", "IMGROW0_8X8", "IMGCOL0_8X8", "TIME"]: assert col in dat.colnames @@ -202,9 +202,9 @@ def test_get_dark_images(mock_dc, mock_img_table, dc_imgs_dn): """ dc_imgs = get_dark_current_imgs( img_table=mock_img_table, - dc_img=mock_dc, - dc_tccd=-10, - t_ccd=np.repeat(-10, 8), + img_dark=mock_dc, + tccd_dark=-10, + t_ccd_vals=np.repeat(-10, 8), t_ccd_times=mock_img_table["TIME"], ) assert dc_imgs.shape == (8, 8, 8) @@ -241,9 +241,9 @@ def test_dc_consistent(): dc = get_dark_cal_props(tstart, "nearest", include_image=True, aca_image=False) dc_imgs = get_dark_current_imgs( img_table_slot, - dc_img=dc["image"], - dc_tccd=dc["t_ccd"], - t_ccd=np.repeat(-9.9, len(img_table_slot)), + img_dark=dc["image"], + tccd_dark=dc["t_ccd"], + t_ccd_vals=np.repeat(-9.9, len(img_table_slot)), t_ccd_times=img_table_slot["TIME"], ) edge_mask_6x6 = np.zeros((8, 8), dtype=bool) @@ -267,7 +267,8 @@ def test_dcsub_aca_images_maude(): # This one is just a regression test tstart = 471139130.93277466 # dwell tstart for obsid 15600 - random obsid tstop = tstart + 20 - imgs_bgsub = _get_dcsub_aca_images(tstart, tstop, source="maude") + table_bgsub = _get_dcsub_aca_images(tstart, tstop, source="maude") + imgs_bgsub = table_bgsub["IMGBGSUB"] exp0 = np.array( [ From f18ea88386583a079eb3754c8047967ffd20555d Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Tue, 24 Sep 2024 20:07:02 -0400 Subject: [PATCH 06/28] Improve docstring and start/stop errors --- chandra_aca/dark_subtract.py | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index b6663180..006e412c 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -92,10 +92,18 @@ def _get_dcsub_aca_images( """ + if start is None and aca_image_table is None: + raise ValueError("start must be defined if aca_image_table is not defined") + if stop is None and aca_image_table is None: + raise ValueError("stop must be defined if aca_image_table is not defined") + if aca_image_table is not None: # If supplied an image table, use that, but confirm it has the right pieces. assert type(aca_image_table) is Table - # Require aca_image_table have columns IMG, IMGROW0_8X8, IMGCOL0_8X8, TIME + + # Require aca_image_table have columns IMG, IMGROW0_8X8, IMGCOL0_8X8, TIME. + # If this is a table of mica data that doesn't have those columns, see + # mica.archive.aca_l0.get_aca_images for how to get the data or use that method. assert "IMG" in aca_image_table.colnames assert "IMGROW0_8X8" in aca_image_table.colnames assert "IMGCOL0_8X8" in aca_image_table.colnames @@ -112,8 +120,6 @@ def _get_dcsub_aca_images( if stop is not None: raise ValueError("Do not define stop if aca_image_table is defined") elif source == "maude": - if start is None or stop is None: - raise ValueError("start and stop must be defined for maude source") img_table = maude_decom.get_aca_images(start, stop) elif source == "mica": img_table = aca_l0.get_aca_images(start, stop) @@ -122,11 +128,6 @@ def _get_dcsub_aca_images( "aca_image_table must be defined or source must be maude or mica" ) - if start is None and aca_image_table is None: - raise ValueError("start must be defined if aca_image_table is not defined") - if stop is None and aca_image_table is None: - raise ValueError("stop must be defined if aca_image_table is not defined") - if start is None: start = img_table["TIME"].min() if stop is None: @@ -252,11 +253,19 @@ def get_dcsub_aca_images( """ Get aca images for a given time range. + One must supply either aca_image_table or start and stop times. If start and stop + times are supplied, an aca image table will be fetched from the maude or mica data + and a table will be returned that includes dark current background subtracted images + in one of the astropy table columns. If an aca_image_table is supplied, the images + within that table will be used as the raw images, and a new table that includes + dark current background subtracted images will be returned. + + Parameters ---------- - start : str, optional + start : CxoTimeLike, optional Start time of the time range. Default is None. - stop : str, optional + stop : CxoTimeLike, optional Stop time of the time range. Default is None. aca_image_table : astropy.table.Table Table including ACA images. From 99361f1dcfb0caee8413315e3b54578d87a47051 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 21 Oct 2024 11:34:20 -0400 Subject: [PATCH 07/28] Use default padding and smoothing in get_tccd_data --- chandra_aca/dark_subtract.py | 57 ++++++++++++++++++++----- chandra_aca/tests/test_dark_subtract.py | 9 ++-- 2 files changed, 52 insertions(+), 14 deletions(-) diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index 006e412c..262a1080 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -1,19 +1,21 @@ import numba import numpy as np import requests.exceptions +import scipy.signal from astropy.table import Table from cheta import fetch_sci -from cxotime import CxoTimeLike +from cxotime import CxoTime, CxoTimeLike from mica.archive import aca_l0 from mica.archive.aca_dark import get_dark_cal_props from ska_helpers import retry +from ska_numpy import smooth from chandra_aca import maude_decom from chandra_aca.dark_model import dark_temp_scale_img @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) -def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): +def get_tccd_data(times: CxoTimeLike, pad=600, smooth_window=30, source="maude"): """ Get the CCD temperature for a given time range. @@ -22,13 +24,17 @@ def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): Parameters ---------- - start : CxoTimeLike - Start time of the time range. - stop : CxoTimeLike - Stop time of the time range. + times: + The times for which to get the CCD temperature data. source : str, optional Source of the CCD temperature data. If 'maude', override the fetch_sci data source to 'maude allow_subset=False'. Else, use the cheta default. + pad : int, optional + Pad the query times by this amount in seconds. Default is 600. + This extra time is not included in the returned data. + smooth_window : int, optional + 3-sample Median filter and then smooth the data using a hanning window of this length. + Generally should not be more than pad * 2 / 32.8 . Returns ------- @@ -37,13 +43,33 @@ def get_tccd_data(start: CxoTimeLike, stop: CxoTimeLike, source="maude"): t_ccd_times : np.array Times corresponding to the CCD temperature values. """ + fetch_start = CxoTime(times[0]).secs + fetch_stop = CxoTime(times[-1]).secs + if pad > 0: + fetch_start -= pad + fetch_stop += pad + if source == "maude": # Override the data_source to be explicit about maude source. with fetch_sci.data_source("maude allow_subset=False"): - dat = fetch_sci.Msid("aacccdpt", start, stop) + dat = fetch_sci.Msid("aacccdpt", fetch_start, fetch_stop) else: - dat = fetch_sci.Msid("aacccdpt", start, stop) - return dat.vals, dat.times + dat = fetch_sci.Msid("aacccdpt", fetch_start, fetch_stop) + + xin = dat.times + yin = dat.vals.copy() + + if smooth_window > 0: + # Filter the data using a 3 point median filter from scipy + yin = scipy.signal.medfilt(yin, kernel_size=3) + # And then smooth it with hanning and window_len = smooth_window + yin = smooth(yin, window_len=smooth_window) + + xout = CxoTime(times).secs + # Interpolate the data to the requested times + vals = np.interp(xout, xin, yin) + + return vals, times @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) @@ -142,8 +168,17 @@ def _get_dcsub_aca_images( assert img_dark.shape == (1024, 1024) + # if the img_table has masked values in TIME, filter the table to remove + # the full rows before proceeding. + if hasattr(img_table["TIME"], "mask"): + # Filter the table in place so that we can unmask/compress TIME safely + img_table = img_table[~img_table["TIME"].mask] + times = img_table["TIME"].compressed() + else: + times = img_table["TIME"] + if t_ccd_vals is None: - t_ccd_vals, t_ccd_times = get_tccd_data(start, stop, source=source) + t_ccd_vals, t_ccd_times = get_tccd_data(times, source=source) imgs_dark = get_dark_current_imgs( img_table, img_dark, tccd_dark, t_ccd_vals, t_ccd_times @@ -153,7 +188,7 @@ def _get_dcsub_aca_images( img_table["IMGBGSUB"] = imgs_bgsub img_table["IMGDARK"] = imgs_dark - img_table["AACCCDPT"] = np.interp(img_table["TIME"], t_ccd_times, t_ccd_vals) + img_table["AACCCDPT"] = t_ccd_vals return img_table diff --git a/chandra_aca/tests/test_dark_subtract.py b/chandra_aca/tests/test_dark_subtract.py index 144096b1..962df866 100644 --- a/chandra_aca/tests/test_dark_subtract.py +++ b/chandra_aca/tests/test_dark_subtract.py @@ -7,6 +7,7 @@ from astropy.table import Table from mica.archive.aca_dark import get_dark_cal_props from numpy import ma +from cxotime import CxoTime from chandra_aca import maude_decom from chandra_aca.dark_subtract import ( @@ -165,15 +166,17 @@ def test_dcsub_aca_images(mock_dc, mock_img_table, dc_imgs_dn): def test_get_tccd_data(): start = "2021:001:00:00:00.000" stop = "2021:001:00:00:30.000" - t_ccd_maude, t_ccd_times_maude = get_tccd_data(start, stop, source="maude") + times = CxoTime.linspace(start, stop, 8) + t_ccd_maude, t_ccd_times_maude = get_tccd_data(times, source="maude") + assert np.all(t_ccd_times_maude == times) assert np.isclose(t_ccd_maude[0], -6.55137778) # Confirm the t_ccd values are the same for maude as default fetch data source # but only bother if cxc is in the sources. Technically this should be a skipif. if "cxc" in cheta.fetch.data_source.sources(): - t_ccd, t_ccd_times = get_tccd_data(start, stop) + t_ccd, t_ccd_times = get_tccd_data(times) assert np.allclose(t_ccd_maude, t_ccd) - assert np.allclose(t_ccd_times_maude, t_ccd_times) + assert np.allclose(t_ccd_times_maude.secs, t_ccd_times.secs) def test_get_aca_images(mock_dc, mock_img_table, dc_imgs_dn): From f341fb48c0bc632988029ea0fd79df0fe4f15375 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 6 Nov 2024 10:42:48 -0500 Subject: [PATCH 08/28] Reorg and progress --- chandra_aca/aca_image.py | 55 ++++++ chandra_aca/dark_subtract.py | 246 +++++------------------- chandra_aca/tests/test_dark_subtract.py | 82 ++++---- 3 files changed, 147 insertions(+), 236 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 44178f89..331301f0 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -5,8 +5,16 @@ from math import floor from pathlib import Path +import mica.archive.aca_dark +import mica.archive.aca_l0 import numba import numpy as np +import requests +from astropy.table import Table +from ska_helpers import retry + +import chandra_aca.dark_subtract +import chandra_aca.maude_decom __all__ = ["ACAImage", "centroid_fm", "AcaPsfLibrary", "EIGHT_LABELS"] @@ -828,3 +836,50 @@ def get_psf_image( out = ACAImage(psf, row0=row0, col0=col0) if aca_image else (psf, row0, col0) return out + + +@retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) +def get_aca_image_table(start, stop, bgsub=True, source="maude", **maude_kwargs): + if source == "maude": + imgs_table = chandra_aca.maude_decom.get_aca_images(start, stop, **maude_kwargs) + t_ccds = chandra_aca.dark_subtract.get_tccd_data( + imgs_table["TIME"], source=source, **maude_kwargs + ) + elif source == "mica": + imgs_table = mica.archive.aca_l0.get_aca_images(start, stop) + t_ccds = chandra_aca.dark_subtract.get_tccd_data( + imgs_table["TIME"], source=source + ) + else: + raise ValueError(f"source must be 'maude' or 'mica', not {source}") + + # Do some unmasking until the maude_decom and aca_l0 interfaces are updated with + # fewer masked columns + imgs_table_unmasked = Table() + for col in imgs_table.colnames: + if col in ["TIME", "INTEG", "IMGROW0_8X8", "IMGCOL0_8X8"]: + imgs_table_unmasked[col] = imgs_table[col].data.data + else: + imgs_table_unmasked[col] = imgs_table[col] + imgs_table = imgs_table_unmasked + + if bgsub: + dark_data = mica.archive.aca_dark.get_dark_cal_props( + imgs_table["TIME"].min(), + select="nearest", + include_image=True, + aca_image=False, + ) + img_dark = dark_data["image"] + tccd_dark = dark_data["ccd_temp"] + imgs_dark = chandra_aca.dark_subtract.get_dark_current_imgs( + imgs_table, img_dark, tccd_dark, t_ccds + ) + imgs_bgsub = imgs_table["IMG"] - imgs_dark + imgs_bgsub.clip(0, None) + + imgs_table["BGSUB"] = imgs_bgsub + imgs_table["DARK"] = imgs_dark + imgs_table["T_CCD"] = t_ccds + + return imgs_table diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index 262a1080..63581978 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -2,213 +2,97 @@ import numpy as np import requests.exceptions import scipy.signal -from astropy.table import Table from cheta import fetch_sci -from cxotime import CxoTime, CxoTimeLike -from mica.archive import aca_l0 -from mica.archive.aca_dark import get_dark_cal_props from ska_helpers import retry from ska_numpy import smooth -from chandra_aca import maude_decom from chandra_aca.dark_model import dark_temp_scale_img @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) -def get_tccd_data(times: CxoTimeLike, pad=600, smooth_window=30, source="maude"): +def get_tccd_data( + times, smooth_window=30, median_window=3, source="maude", maude_channel=None +): """ - Get the CCD temperature for a given time range. + Get the CCD temperature for given times and interpolate and smooth. This is a light wrapper around cheta.fetch_sci.Msid("aacccdpt", start, stop) to handle an option to use the maude data source in an explicit way. Parameters ---------- - times: + times: np.array float of Chandra seconds The times for which to get the CCD temperature data. source : str, optional Source of the CCD temperature data. If 'maude', override the fetch_sci data source to 'maude allow_subset=False'. Else, use the cheta default. - pad : int, optional - Pad the query times by this amount in seconds. Default is 600. - This extra time is not included in the returned data. + median_window : int, optional + 3-sample Median filter to to remove outliers. smooth_window : int, optional - 3-sample Median filter and then smooth the data using a hanning window of this length. - Generally should not be more than pad * 2 / 32.8 . + Smooth the data using a hanning window of this length in samples. + Returns ------- - t_ccd_vals : np.array + vals : np.array CCD temperature values. - t_ccd_times : np.array - Times corresponding to the CCD temperature values. """ - fetch_start = CxoTime(times[0]).secs - fetch_stop = CxoTime(times[-1]).secs - if pad > 0: - fetch_start -= pad - fetch_stop += pad + + pad = 600 # 600 seconds default pad + + # increase the pad if the smooth window is large + pad = np.max([smooth_window * 32.8 / 2, pad]) + + fetch_start = times[0] - pad + fetch_stop = times[-1] + pad if source == "maude": - # Override the data_source to be explicit about maude source. - with fetch_sci.data_source("maude allow_subset=False"): + # Override the cheta data_source to be explicit about maude source. + data_source = "maude allow_subset=False" + if maude_channel is not None: + data_source += f" channel={maude_channel}" + with fetch_sci.data_source(data_source): dat = fetch_sci.Msid("aacccdpt", fetch_start, fetch_stop) else: dat = fetch_sci.Msid("aacccdpt", fetch_start, fetch_stop) - xin = dat.times yin = dat.vals.copy() - if smooth_window > 0: + if median_window > 0: # Filter the data using a 3 point median filter from scipy - yin = scipy.signal.medfilt(yin, kernel_size=3) + yin = scipy.signal.medfilt(yin, kernel_size=median_window) + + if smooth_window > 0: # And then smooth it with hanning and window_len = smooth_window yin = smooth(yin, window_len=smooth_window) - xout = CxoTime(times).secs # Interpolate the data to the requested times - vals = np.interp(xout, xin, yin) - - return vals, times - - -@retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) -def _get_dcsub_aca_images( - start: CxoTimeLike = None, - stop: CxoTimeLike = None, - aca_image_table=None, - t_ccd_vals=None, - t_ccd_times=None, - img_dark=None, - tccd_dark=None, - source="maude", -): - """ - Get dark current subtracted ACA images for a given time range. - - This private method has expanded options for unit testing. - Users should use get_dcsub_aca_images instead. - - Parameters - ---------- - start : CxoTimeLike, optional - Start time of the time range. Default is None. - stop : CxoTimeLike, optional - Stop time of the time range. Default is None. - aca_image_table : astropy table of ACA images, optional - - t_ccd_vals : array-like, optional - Array of CCD temperatures. Default is None. - t_ccd_times : array-like, optional - Array of CCD temperature times. Default is None. - img_dark : array-like, optional - Dark current image. Default is None. - tccd_dark : float, optional - Dark current CCD temperature. Default is None. - source : str, optional - Source of the ACA images. Can be 'maude' or 'mica'. Default is 'maude'. - - Returns - ------- - astropy.table.Table - Table of aca image data including raw and dark current subtracted images. - Raw images are in column 'IMG', dark current subtracted images are in column 'IMGBGSUB', - dark current cut-out images are in column 'IMGDARK', and interpolated - CCD temperature values are in column 'AACCCDPT'. - - """ + vals = np.interp(times, dat.times, yin) - if start is None and aca_image_table is None: - raise ValueError("start must be defined if aca_image_table is not defined") - if stop is None and aca_image_table is None: - raise ValueError("stop must be defined if aca_image_table is not defined") - - if aca_image_table is not None: - # If supplied an image table, use that, but confirm it has the right pieces. - assert type(aca_image_table) is Table - - # Require aca_image_table have columns IMG, IMGROW0_8X8, IMGCOL0_8X8, TIME. - # If this is a table of mica data that doesn't have those columns, see - # mica.archive.aca_l0.get_aca_images for how to get the data or use that method. - assert "IMG" in aca_image_table.colnames - assert "IMGROW0_8X8" in aca_image_table.colnames - assert "IMGCOL0_8X8" in aca_image_table.colnames - assert "TIME" in aca_image_table.colnames - # require that IMGROW0_8X8 and IMGCOL0_8X8 are masked columns - assert hasattr(aca_image_table["IMGROW0_8X8"], "mask") - assert hasattr(aca_image_table["IMGCOL0_8X8"], "mask") - # Make a copy so we don't modify the input table - img_table = aca_image_table.copy() - - # And it doesn't make sense to also define start and stop so raise errors if so - if start is not None: - raise ValueError("Do not define start if aca_image_table is defined") - if stop is not None: - raise ValueError("Do not define stop if aca_image_table is defined") - elif source == "maude": - img_table = maude_decom.get_aca_images(start, stop) - elif source == "mica": - img_table = aca_l0.get_aca_images(start, stop) - else: - raise ValueError( - "aca_image_table must be defined or source must be maude or mica" - ) - - if start is None: - start = img_table["TIME"].min() - if stop is None: - stop = img_table["TIME"].max() - - if img_dark is None: - dark_data = get_dark_cal_props( - start, select="nearest", include_image=True, aca_image=False - ) - img_dark = dark_data["image"] - tccd_dark = dark_data["ccd_temp"] - - assert img_dark.shape == (1024, 1024) - - # if the img_table has masked values in TIME, filter the table to remove - # the full rows before proceeding. - if hasattr(img_table["TIME"], "mask"): - # Filter the table in place so that we can unmask/compress TIME safely - img_table = img_table[~img_table["TIME"].mask] - times = img_table["TIME"].compressed() - else: - times = img_table["TIME"] + return vals - if t_ccd_vals is None: - t_ccd_vals, t_ccd_times = get_tccd_data(times, source=source) - imgs_dark = get_dark_current_imgs( - img_table, img_dark, tccd_dark, t_ccd_vals, t_ccd_times - ) +def get_aca_images_bgd_sub(img_table, t_ccd_vals, img_dark, tccd_dark): + imgs_dark = get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccd_vals) imgs_bgsub = img_table["IMG"] - imgs_dark imgs_bgsub.clip(0, None) - img_table["IMGBGSUB"] = imgs_bgsub - img_table["IMGDARK"] = imgs_dark - img_table["AACCCDPT"] = t_ccd_vals - - return img_table + return imgs_bgsub, imgs_dark -def get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccd_vals, t_ccd_times): +def get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccds): """ Get the scaled dark current values for a table of ACA images. Parameters ---------- img_table : astropy.table.Table - A table containing the ACA images and expected metadata. img_dark : 1024x1024 array The dark calibration image. tccd_dark : float The reference temperature of the dark calibration image. - t_ccd_vals : array - The cheta temperature values covering the time range of the img_table. - t_ccd_times : array - The corresponding times for the temperature values. + t_ccds : array + The cheta temperature values at the times of img_table. Returns ------- @@ -217,21 +101,23 @@ def get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccd_vals, t_ccd_time in the img_table in e-/s. """ - assert img_dark.shape == (1024, 1024) - assert len(t_ccd_vals) == len(t_ccd_times) + if len(img_table) != len(t_ccds): + raise ValueError("img_table and t_ccds must have the same length") + + if img_dark.shape != (1024, 1024): + raise ValueError("Dark image shape is not 1024x1024") imgs_dark_unscaled = get_dark_backgrounds( img_dark, - img_table["IMGROW0_8X8"].data.filled(1025), - img_table["IMGCOL0_8X8"].data.filled(1025), + img_table["IMGROW0_8X8"], + img_table["IMGCOL0_8X8"], ) - dt_ccd = np.interp(img_table["TIME"], t_ccd_times, t_ccd_vals) imgs_dark = np.zeros((len(img_table), 8, 8), dtype=np.float64) # Scale the dark current at each dark cal 8x8 image to the ccd temperature and e-/s for i, (img_dark_unscaled, t_ccd_i, img) in enumerate( - zip(imgs_dark_unscaled, dt_ccd, img_table, strict=True) + zip(imgs_dark_unscaled, t_ccds, img_table, strict=True) ): img_scaled = dark_temp_scale_img(img_dark_unscaled, tccd_dark, t_ccd_i) # Default to using 1.696 integ time if not present in the image table @@ -276,49 +162,3 @@ def staggered_aca_slice(array_in, array_out, row, col): raw_dark_img.astype(float), imgs_dark, 512 + imgrow0, 512 + imgcol0 ) return imgs_dark - - -@retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) -def get_dcsub_aca_images( - start: CxoTimeLike = None, - stop: CxoTimeLike = None, - aca_image_table=None, - source="maude", -): - """ - Get aca images for a given time range. - - One must supply either aca_image_table or start and stop times. If start and stop - times are supplied, an aca image table will be fetched from the maude or mica data - and a table will be returned that includes dark current background subtracted images - in one of the astropy table columns. If an aca_image_table is supplied, the images - within that table will be used as the raw images, and a new table that includes - dark current background subtracted images will be returned. - - - Parameters - ---------- - start : CxoTimeLike, optional - Start time of the time range. Default is None. - stop : CxoTimeLike, optional - Stop time of the time range. Default is None. - aca_image_table : astropy.table.Table - Table including ACA images. - source : str, optional - Source of the ACA images. Can be 'maude' or 'mica'. Default is 'maude'. - - Returns - ------- - astropy.table.Table - Table of aca image data including raw and dark current subtracted images. - Raw images are in column 'IMG', dark current subtracted images are in column 'IMGBGSUB', - dark current cut-out images are in column 'IMGDARK', and interpolated - CCD temperature values are in column 'AACCCDPT'. - - """ - return _get_dcsub_aca_images( - start=start, - stop=stop, - aca_image_table=aca_image_table, - source=source, - ) diff --git a/chandra_aca/tests/test_dark_subtract.py b/chandra_aca/tests/test_dark_subtract.py index 962df866..86d92127 100644 --- a/chandra_aca/tests/test_dark_subtract.py +++ b/chandra_aca/tests/test_dark_subtract.py @@ -5,13 +5,13 @@ import numpy as np import pytest from astropy.table import Table -from mica.archive.aca_dark import get_dark_cal_props -from numpy import ma from cxotime import CxoTime +from mica.archive.aca_dark import get_dark_cal_props from chandra_aca import maude_decom +from chandra_aca.aca_image import get_aca_image_table from chandra_aca.dark_subtract import ( - _get_dcsub_aca_images, + get_aca_images_bgd_sub, get_dark_backgrounds, get_dark_current_imgs, get_tccd_data, @@ -44,8 +44,8 @@ def mock_img_table(): """ img_table = Table() img_table["TIME"] = np.arange(8) - img_table["IMGROW0_8X8"] = ma.arange(8) - 512 - img_table["IMGCOL0_8X8"] = ma.arange(8) - 512 + img_table["IMGROW0_8X8"] = np.arange(8) - 512 + img_table["IMGCOL0_8X8"] = np.arange(8) - 512 img_table["IMG"] = np.ones((8, 8, 8)) * 16 * 1.696 / 5 img_table["IMGNUM"] = 0 return img_table @@ -149,14 +149,12 @@ def test_dcsub_aca_images(mock_dc, mock_img_table, dc_imgs_dn): Confirm that the pattern of background subtraction is correct. """ - table_bgsub = _get_dcsub_aca_images( - aca_image_table=mock_img_table, + imgs_bgsub, _ = get_aca_images_bgd_sub( + img_table=mock_img_table, img_dark=mock_dc, tccd_dark=-10, t_ccd_vals=np.repeat(-10, 8), - t_ccd_times=mock_img_table["TIME"], ) - imgs_bgsub = table_bgsub["IMGBGSUB"] assert imgs_bgsub.shape == (8, 8, 8) # Note that the mock unsubtracted data is originally 16 * 1.696 / 5 exp = 16 - dc_imgs_dn @@ -167,35 +165,29 @@ def test_get_tccd_data(): start = "2021:001:00:00:00.000" stop = "2021:001:00:00:30.000" times = CxoTime.linspace(start, stop, 8) - t_ccd_maude, t_ccd_times_maude = get_tccd_data(times, source="maude") - assert np.all(t_ccd_times_maude == times) + t_ccd_maude = get_tccd_data(times.secs, source="maude") assert np.isclose(t_ccd_maude[0], -6.55137778) # Confirm the t_ccd values are the same for maude as default fetch data source # but only bother if cxc is in the sources. Technically this should be a skipif. if "cxc" in cheta.fetch.data_source.sources(): - t_ccd, t_ccd_times = get_tccd_data(times) + t_ccd = get_tccd_data(times.secs) assert np.allclose(t_ccd_maude, t_ccd) - assert np.allclose(t_ccd_times_maude.secs, t_ccd_times.secs) def test_get_aca_images(mock_dc, mock_img_table, dc_imgs_dn): """ Confirm the pattern of dark current images matches the reference. """ - dat = _get_dcsub_aca_images( - aca_image_table=mock_img_table, + imgs_bgsub, _ = get_aca_images_bgd_sub( + img_table=mock_img_table, img_dark=mock_dc, tccd_dark=-10, t_ccd_vals=np.repeat(-10, 8), - t_ccd_times=mock_img_table["TIME"], ) - for col in ["IMG", "IMGROW0_8X8", "IMGCOL0_8X8", "TIME"]: - assert col in dat.colnames - assert np.allclose(dat[col], mock_img_table[col], atol=1e-6, rtol=0) - assert dat["IMGBGSUB"].shape == (8, 8, 8) + assert imgs_bgsub.shape == (8, 8, 8) exp = 16 - dc_imgs_dn - assert np.allclose(dat["IMGBGSUB"] * 5 / 1.696, exp, atol=1e-6, rtol=0) + assert np.allclose(imgs_bgsub * 5 / 1.696, exp, atol=1e-6, rtol=0) def test_get_dark_images(mock_dc, mock_img_table, dc_imgs_dn): @@ -207,8 +199,7 @@ def test_get_dark_images(mock_dc, mock_img_table, dc_imgs_dn): img_table=mock_img_table, img_dark=mock_dc, tccd_dark=-10, - t_ccd_vals=np.repeat(-10, 8), - t_ccd_times=mock_img_table["TIME"], + t_ccds=np.repeat(-10, 8), ) assert dc_imgs.shape == (8, 8, 8) dc_imgs_raw = dc_imgs * 5 / 1.696 @@ -221,8 +212,8 @@ def test_get_dark_backgrounds(mock_dc, mock_img_table, dc_imgs_dn): """ dc_imgs_raw = get_dark_backgrounds( mock_dc, - mock_img_table["IMGROW0_8X8"].data.filled(1025), - mock_img_table["IMGCOL0_8X8"].data.filled(1025), + mock_img_table["IMGROW0_8X8"], + mock_img_table["IMGCOL0_8X8"], ) assert dc_imgs_raw.shape == (8, 8, 8) assert np.allclose(dc_imgs_raw, dc_imgs_dn, atol=1e-6, rtol=0) @@ -234,20 +225,23 @@ def test_dc_consistent(): edges of a 6x6. """ tstart = 746484519.429 - img_table = maude_decom.get_aca_images(tstart, tstart + 30) + imgs_table_masked = maude_decom.get_aca_images(tstart, tstart + 30) + imgs_table = Table() + imgs_table["IMG"] = imgs_table_masked["IMG"] + for col in ["TIME", "IMGROW0_8X8", "IMGCOL0_8X8", "IMGNUM", "INTEG"]: + imgs_table[col] = imgs_table_masked[col].data.data # This slot and this time show an overall decent correlation of dark current # and edge pixel values in the 6x6 image - but don't make a great test. slot = 4 - img_table_slot = img_table[img_table["IMGNUM"] == slot] - img = img_table_slot[0]["IMG"] + imgs_table_slot = imgs_table[imgs_table["IMGNUM"] == slot] + img = imgs_table_slot[0]["IMG"] dc = get_dark_cal_props(tstart, "nearest", include_image=True, aca_image=False) dc_imgs = get_dark_current_imgs( - img_table_slot, + imgs_table_slot, img_dark=dc["image"], tccd_dark=dc["t_ccd"], - t_ccd_vals=np.repeat(-9.9, len(img_table_slot)), - t_ccd_times=img_table_slot["TIME"], + t_ccds=np.repeat(-9.9, len(imgs_table_slot)), ) edge_mask_6x6 = np.zeros((8, 8), dtype=bool) edge_mask_6x6[1, 2:6] = True @@ -270,8 +264,28 @@ def test_dcsub_aca_images_maude(): # This one is just a regression test tstart = 471139130.93277466 # dwell tstart for obsid 15600 - random obsid tstop = tstart + 20 - table_bgsub = _get_dcsub_aca_images(tstart, tstop, source="maude") - imgs_bgsub = table_bgsub["IMGBGSUB"] + imgs_table_masked = maude_decom.get_aca_images(tstart, tstop) + imgs_table = Table() + imgs_table["IMG"] = imgs_table_masked["IMG"] + for col in ["TIME", "IMGROW0_8X8", "IMGCOL0_8X8", "IMGNUM", "INTEG"]: + imgs_table[col] = imgs_table_masked[col].data.data + + import mica.archive.aca_dark + + import chandra_aca.dark_subtract + + t_ccds = chandra_aca.dark_subtract.get_tccd_data(imgs_table["TIME"], source="maude") + dark_data = mica.archive.aca_dark.get_dark_cal_props( + imgs_table["TIME"].min(), select="nearest", include_image=True, aca_image=False + ) + img_dark = dark_data["image"] + tccd_dark = dark_data["ccd_temp"] + + imgs_bgsub, _ = get_aca_images_bgd_sub( + imgs_table, img_dark=img_dark, tccd_dark=tccd_dark, t_ccd_vals=t_ccds + ) + + imgs_bgsub_table = get_aca_image_table(tstart, tstop, source="maude") exp0 = np.array( [ @@ -286,6 +300,7 @@ def test_dcsub_aca_images_maude(): ] ) np.allclose(exp0, imgs_bgsub[0][0].filled(0), atol=1, rtol=0) + np.allclose(exp0, imgs_bgsub_table["BGSUB"][0][0].filled(0), atol=1, rtol=0) # slot 4 is 6x6 and masked exp4 = np.array( @@ -301,3 +316,4 @@ def test_dcsub_aca_images_maude(): ] ) np.allclose(exp4, imgs_bgsub[4][0].filled(0), atol=1, rtol=0) + np.allclose(exp4, imgs_bgsub_table["BGSUB"][4][0].filled(0), atol=1, rtol=0) From 4080ac5227fe5483d8a773b71604ce9ccb1a4e6e Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Tue, 19 Nov 2024 15:51:49 -0500 Subject: [PATCH 09/28] Move some imports to avoid circular issues --- chandra_aca/aca_image.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 331301f0..c51bba1c 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -5,17 +5,12 @@ from math import floor from pathlib import Path -import mica.archive.aca_dark -import mica.archive.aca_l0 import numba import numpy as np import requests from astropy.table import Table from ska_helpers import retry -import chandra_aca.dark_subtract -import chandra_aca.maude_decom - __all__ = ["ACAImage", "centroid_fm", "AcaPsfLibrary", "EIGHT_LABELS"] EIGHT_LABELS = np.array( @@ -840,6 +835,12 @@ def get_psf_image( @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) def get_aca_image_table(start, stop, bgsub=True, source="maude", **maude_kwargs): + import mica.archive.aca_dark + import mica.archive.aca_l0 + + import chandra_aca.dark_subtract + import chandra_aca.maude_decom + if source == "maude": imgs_table = chandra_aca.maude_decom.get_aca_images(start, stop, **maude_kwargs) t_ccds = chandra_aca.dark_subtract.get_tccd_data( From f0d8d90451716179a59d70282f801cb4b04c94ea Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Tue, 19 Nov 2024 19:48:41 -0500 Subject: [PATCH 10/28] Add more docstrings --- chandra_aca/aca_image.py | 29 ++++++++++++++++++++++++++++- chandra_aca/dark_subtract.py | 22 ++++++++++++++++++++++ 2 files changed, 50 insertions(+), 1 deletion(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index c51bba1c..6bce06fd 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -9,6 +9,7 @@ import numpy as np import requests from astropy.table import Table +from cxotime import CxoTimeLike from ska_helpers import retry __all__ = ["ACAImage", "centroid_fm", "AcaPsfLibrary", "EIGHT_LABELS"] @@ -834,13 +835,38 @@ def get_psf_image( @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) -def get_aca_image_table(start, stop, bgsub=True, source="maude", **maude_kwargs): +def get_aca_image_table( + start: CxoTimeLike, stop: CxoTimeLike, bgsub=True, source="maude", **maude_kwargs +): + """ + Get ACA images and ancillary data from either the mica or maude data sources. + + Parameters + ---------- + start : CxoTimeLike + start time + stop : CxoTimeLike + stop time (CXC sec) + bgsub : bool + include background subtracted images in output table + source : str + 'maude' or 'mica' + maude_kwargs + additional kwargs for maude data source + + Returns + ------- + Table of ACA images and ancillary data. If bgsub is True then the table + will include columns for the background subtracted image ('BGSUB'), the dark current + image ('DARK'), and the CCD temperature ('T_CCD'). + """ import mica.archive.aca_dark import mica.archive.aca_l0 import chandra_aca.dark_subtract import chandra_aca.maude_decom + # Get aca images over the time range if source == "maude": imgs_table = chandra_aca.maude_decom.get_aca_images(start, stop, **maude_kwargs) t_ccds = chandra_aca.dark_subtract.get_tccd_data( @@ -864,6 +890,7 @@ def get_aca_image_table(start, stop, bgsub=True, source="maude", **maude_kwargs) imgs_table_unmasked[col] = imgs_table[col] imgs_table = imgs_table_unmasked + # Get background subtracted values if bgsub is True if bgsub: dark_data = mica.archive.aca_dark.get_dark_cal_props( imgs_table["TIME"].min(), diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index 63581978..0142456c 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -73,6 +73,28 @@ def get_tccd_data( def get_aca_images_bgd_sub(img_table, t_ccd_vals, img_dark, tccd_dark): + """ + Get the background subtracted ACA images. + + Parameters + ---------- + img_table : astropy.table.Table + The table of ACA images. + t_ccd_vals : np.array + The CCD temperature values at the times of the ACA images. + img_dark : np.array + The dark calibration image. Must be 1024x1024. + tccd_dark : float + The reference temperature of the dark calibration image. + + Returns + ------- + tuple (imgs_bgsub, imgs_dark) + imgs_bgsub : np.array + The background subtracted ACA images. + imgs_dark : np.array + The dark current images. + """ imgs_dark = get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccd_vals) imgs_bgsub = img_table["IMG"] - imgs_dark imgs_bgsub.clip(0, None) From 84d95c31472632f9382bc88d48be8681c536eaf2 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 20 Nov 2024 14:36:11 -0500 Subject: [PATCH 11/28] Add more tests --- chandra_aca/tests/test_aca_image.py | 113 ++++++++++++++++++++++++ chandra_aca/tests/test_dark_subtract.py | 9 +- 2 files changed, 117 insertions(+), 5 deletions(-) diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index af91b2a1..3e9122b1 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -1,6 +1,13 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +from pathlib import Path + +import mica.common import numpy as np import pytest +from cxotime import CxoTime +from mica.archive.aca_dark import get_dark_cal_props + +import chandra_aca.aca_image from ..aca_image import ACAImage, centroid_fm @@ -482,3 +489,109 @@ def test_flicker_test_sequence(): ] ), ) + + +def images_check_range(start, stop, img_table): + tstart = CxoTime(start).secs + tstop = CxoTime(stop).secs + + # Check that start and stop of the range are within 4.1 secs of the requested times + assert np.all(np.abs(img_table["TIME"][0] - tstart) < 4.1) + assert np.all(np.abs(img_table["TIME"][-1] - tstop) < 4.1) + + # Check that all the times in the table are within the requested range + assert np.all((img_table["TIME"] >= tstart) & (img_table["TIME"] < tstop)) + + # Check that all the slots are in there + assert np.all(np.isin(np.arange(8), img_table["IMGNUM"])) + + # Check that if the table has BGSUB column that IMG - DARK = BGSUB + if "BGSUB" in img_table.colnames: + assert np.allclose(img_table["IMG"] - img_table["DARK"], img_table["BGSUB"]) + + # If there's a DARK column, then check it against the dark image from the dark cal + if "DARK" in img_table.colnames: + from mica.archive.aca_dark import get_dark_cal_props + + dc = get_dark_cal_props( + tstart, select="nearest", include_image=True, aca_image=True + ) + full_img_dark = dc["image"] + tccd_dark = dc["ccd_temp"] + for row in img_table: + dark_row = row["DARK"] + tccd_row = row["T_CCD"] + row8x8 = row["IMGROW0_8X8"] + col8x8 = row["IMGCOL0_8X8"] + dark_ref = full_img_dark.aca[row8x8 : row8x8 + 8, col8x8 : col8x8 + 8] + from chandra_aca.dark_model import dark_temp_scale_img + + dark_scale = dark_temp_scale_img(dark_ref, tccd_dark, tccd_row) + # Default to using 1.696 integ time if not present in the image table + integ = row["INTEG"] if "INTEG" in img_table.colnames else 1.696 + img_dark = dark_scale * integ / 5 + assert np.allclose(img_dark, dark_row) + + +def test_get_aca_image_table_maude(): + """Test get_aca_image_table in the maude mode + + This checks that the dark images are reasonable and the IMG data matches with and without the bgsub. + """ + tstart = "2021:001:00:00:00" + tstop = "2021:001:00:01:00" + img_table_maude = chandra_aca.aca_image.get_aca_image_table( + tstart, tstop, source="maude", bgsub=False + ) + images_check_range(tstart, tstop, img_table_maude) + img_table_maude_bgsub = chandra_aca.aca_image.get_aca_image_table( + tstart, tstop, source="maude", bgsub=True + ) + images_check_range(tstart, tstop, img_table_maude_bgsub) + assert np.allclose(img_table_maude["IMG"], img_table_maude_bgsub["IMG"]) + assert np.allclose(img_table_maude["TIME"], img_table_maude_bgsub["TIME"]) + + +HAS_ACA0_ARCHIVE = (Path(mica.common.MICA_ARCHIVE) / "aca0").exists() + + +@pytest.mark.skipif(not HAS_ACA0_ARCHIVE, reason="No ACA0 archive") +def test_get_aca_image_table_mica(): + """Test get_aca_image_table in the mica mode + + This checks that the dark images are reasonable and the answers match maude. + """ + tstart = "2012:180:00:00:00" + tstop = "2012:180:00:01:00" + + img_table_mica = chandra_aca.aca_image.get_aca_image_table( + tstart, tstop, source="mica", bgsub=False + ) + images_check_range(tstart, tstop, img_table_mica) + img_table_mica_bgsub = chandra_aca.aca_image.get_aca_image_table( + tstart, tstop, source="mica", bgsub=True + ) + images_check_range(tstart, tstop, img_table_mica_bgsub) + + # Get maude data too for comparison + img_table_maude = chandra_aca.aca_image.get_aca_image_table( + tstart, tstop, source="maude", bgsub=False + ) + img_table_maude_bgsub = chandra_aca.aca_image.get_aca_image_table( + tstart, tstop, source="maude", bgsub=True + ) + + # Check that the two mica tables are the same in the key columns + assert np.allclose(img_table_mica["IMG"], img_table_mica_bgsub["IMG"]) + assert np.allclose(img_table_mica["TIME"], img_table_mica_bgsub["TIME"]) + + # Check that the tables are the same + img_table_maude.sort(["TIME", "IMGNUM"]) + img_table_maude_bgsub.sort(["TIME", "IMGNUM"]) + img_table_mica.sort(["TIME", "IMGNUM"]) + img_table_mica_bgsub.sort(["TIME", "IMGNUM"]) + + assert np.allclose(img_table_maude["IMG"], img_table_mica["IMG"]) + assert np.allclose(img_table_maude["TIME"], img_table_mica["TIME"]) + assert np.allclose(img_table_maude_bgsub["IMG"], img_table_mica_bgsub["IMG"]) + assert np.allclose(img_table_maude_bgsub["TIME"], img_table_mica_bgsub["TIME"]) diff --git a/chandra_aca/tests/test_dark_subtract.py b/chandra_aca/tests/test_dark_subtract.py index 86d92127..3a703aae 100644 --- a/chandra_aca/tests/test_dark_subtract.py +++ b/chandra_aca/tests/test_dark_subtract.py @@ -1,7 +1,4 @@ -from pathlib import Path - import cheta.fetch -import mica.common import numpy as np import pytest from astropy.table import Table @@ -17,8 +14,6 @@ get_tccd_data, ) -HAS_ACA0_ARCHIVE = (Path(mica.common.MICA_ARCHIVE) / "aca0").exists() - # Make a test fixture for the mock dark current image @pytest.fixture @@ -206,6 +201,10 @@ def test_get_dark_images(mock_dc, mock_img_table, dc_imgs_dn): assert np.allclose(dc_imgs_raw, dc_imgs_dn, atol=1e-6, rtol=0) +def test_dark_map_at_edges(): + pass + + def test_get_dark_backgrounds(mock_dc, mock_img_table, dc_imgs_dn): """ Confirm the pattern of dark current matches the reference. From 149c9c2a302b91f455d6774c9c3f3bb9a3ae53f3 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 20 Nov 2024 14:40:40 -0500 Subject: [PATCH 12/28] Update pre-commit config --- .pre-commit-config.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index cafcca2c..ffda0758 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,12 +1,12 @@ repos: - repo: https://github.com/psf/black - rev: 23.11.0 + rev: 24.10.0 hooks: - id: black language_version: python3.11 - repo: https://github.com/pycqa/isort - rev: 5.12.0 + rev: 5.13.2 hooks: - id: isort name: isort (python) From de58e09f8904dc2cff04251e50778828f3cfe058 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 20 Nov 2024 14:41:58 -0500 Subject: [PATCH 13/28] Remove old import --- chandra_aca/tests/test_aca_image.py | 1 - 1 file changed, 1 deletion(-) diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index 3e9122b1..4d435dec 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -511,7 +511,6 @@ def images_check_range(start, stop, img_table): # If there's a DARK column, then check it against the dark image from the dark cal if "DARK" in img_table.colnames: - from mica.archive.aca_dark import get_dark_cal_props dc = get_dark_cal_props( tstart, select="nearest", include_image=True, aca_image=True From 3edd838b3f9d5aa1570684da51713f5dbd76ae7b Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 20 Nov 2024 14:45:55 -0500 Subject: [PATCH 14/28] Fix long line --- chandra_aca/tests/test_aca_image.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index 4d435dec..4c8fa526 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -535,7 +535,8 @@ def images_check_range(start, stop, img_table): def test_get_aca_image_table_maude(): """Test get_aca_image_table in the maude mode - This checks that the dark images are reasonable and the IMG data matches with and without the bgsub. + This checks that the dark images are reasonable and the IMG data matches with + and without the bgsub. """ tstart = "2021:001:00:00:00" tstop = "2021:001:00:01:00" From f263a89c0b5c6075ce08987da9692398e1a771be Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 25 Nov 2024 11:07:52 -0500 Subject: [PATCH 15/28] Minor renaming and doc improvements from feedback --- chandra_aca/aca_image.py | 40 ++++++++++++++++++------- chandra_aca/dark_subtract.py | 36 ++++++++++++++-------- chandra_aca/tests/test_aca_image.py | 20 ++++++------- chandra_aca/tests/test_dark_subtract.py | 4 +-- docs/dark_subtract.rst | 5 ++++ docs/index.rst | 9 ++++++ 6 files changed, 79 insertions(+), 35 deletions(-) create mode 100644 docs/dark_subtract.rst diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 6bce06fd..8a4da1d4 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -12,7 +12,7 @@ from cxotime import CxoTimeLike from ska_helpers import retry -__all__ = ["ACAImage", "centroid_fm", "AcaPsfLibrary", "EIGHT_LABELS"] +__all__ = ["ACAImage", "centroid_fm", "AcaPsfLibrary", "EIGHT_LABELS", "get_aca_images"] EIGHT_LABELS = np.array( [ @@ -835,11 +835,29 @@ def get_psf_image( @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) -def get_aca_image_table( +def get_aca_images( start: CxoTimeLike, stop: CxoTimeLike, bgsub=True, source="maude", **maude_kwargs -): +) -> Table: """ - Get ACA images and ancillary data from either the mica or maude data sources. + Get ACA images and ancillary data from either the maude or cxc data sources. + + The returned table of ACA images and ancillary data will include the default + columns returned by chandra_aca.maude_decom.get_aca_images or + mica.archive.aca_l0.get_aca_images. If bgsub is True (the default) then the + table will also include columns:: + + name dtype unit + --------------------- ------- ----------- + IMG_BGSUB float64 DN + IMG_DARK float64 DN + T_CCD_SMOOTH float32 degC + + where: + + - 'IMG_BGSUB': background subtracted image + - 'IMG_DARK': dark current image + - 'T_CCD_SMOOTH': smoothed CCD temperature + Parameters ---------- @@ -850,15 +868,15 @@ def get_aca_image_table( bgsub : bool include background subtracted images in output table source : str - 'maude' or 'mica' + Data source for image and temperature telemetry ('maude' or 'cxc'). For 'cxc', + the image telemetry is from mica and temperature telemetry is from CXC L0 via cheta. maude_kwargs additional kwargs for maude data source Returns ------- - Table of ACA images and ancillary data. If bgsub is True then the table - will include columns for the background subtracted image ('BGSUB'), the dark current - image ('DARK'), and the CCD temperature ('T_CCD'). + imgs_table : astropy.table.Table + Table of ACA images and ancillary data. """ import mica.archive.aca_dark import mica.archive.aca_l0 @@ -906,8 +924,8 @@ def get_aca_image_table( imgs_bgsub = imgs_table["IMG"] - imgs_dark imgs_bgsub.clip(0, None) - imgs_table["BGSUB"] = imgs_bgsub - imgs_table["DARK"] = imgs_dark - imgs_table["T_CCD"] = t_ccds + imgs_table["IMG_BGSUB"] = imgs_bgsub + imgs_table["IMG_DARK"] = imgs_dark + imgs_table["T_CCD_SMOOTH"] = t_ccds return imgs_table diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index 0142456c..bfe84413 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -8,6 +8,13 @@ from chandra_aca.dark_model import dark_temp_scale_img +__all__ = [ + "get_tccd_data", + "get_aca_images_bgd_sub", + "get_dark_current_imgs", + "get_dark_backgrounds", +] + @retry.retry(exceptions=requests.exceptions.RequestException, delay=5, tries=3) def get_tccd_data( @@ -25,11 +32,13 @@ def get_tccd_data( The times for which to get the CCD temperature data. source : str, optional Source of the CCD temperature data. If 'maude', override the fetch_sci data source - to 'maude allow_subset=False'. Else, use the cheta default. + to 'maude allow_subset=False'. If 'cxc' use the local cxc archive. Default is 'maude'. median_window : int, optional - 3-sample Median filter to to remove outliers. + Median filter window to remove outliers. Default is 3. smooth_window : int, optional - Smooth the data using a hanning window of this length in samples. + Smooth the data using a hanning window of this length in samples. Default is 30. + maude_channel : str, optional + The maude channel to use. Default is None. Returns @@ -53,13 +62,16 @@ def get_tccd_data( data_source += f" channel={maude_channel}" with fetch_sci.data_source(data_source): dat = fetch_sci.Msid("aacccdpt", fetch_start, fetch_stop) + elif source == "cxc": + with fetch_sci.data_source("cxc"): + dat = fetch_sci.Msid("aacccdpt", fetch_start, fetch_stop) else: - dat = fetch_sci.Msid("aacccdpt", fetch_start, fetch_stop) + raise ValueError(f"Unknown source: {source}") yin = dat.vals.copy() if median_window > 0: - # Filter the data using a 3 point median filter from scipy + # Filter the data using a median filter from scipy yin = scipy.signal.medfilt(yin, kernel_size=median_window) if smooth_window > 0: @@ -81,19 +93,19 @@ def get_aca_images_bgd_sub(img_table, t_ccd_vals, img_dark, tccd_dark): img_table : astropy.table.Table The table of ACA images. t_ccd_vals : np.array - The CCD temperature values at the times of the ACA images. + The CCD temperature values at the times of the ACA images in deg C. img_dark : np.array - The dark calibration image. Must be 1024x1024. + The dark calibration image. Must be 1024x1024. In e-/s. tccd_dark : float - The reference temperature of the dark calibration image. + The reference temperature of the dark calibration image in deg C. Returns ------- tuple (imgs_bgsub, imgs_dark) imgs_bgsub : np.array - The background subtracted ACA images. + The background subtracted ACA images in DN. imgs_dark : np.array - The dark current images. + The dark current images in DN. """ imgs_dark = get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccd_vals) imgs_bgsub = img_table["IMG"] - imgs_dark @@ -120,7 +132,7 @@ def get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccds): ------- imgs_dark : np.array An array containing the temperature scaled dark current for each ACA image - in the img_table in e-/s. + in the img_table in DN. """ if len(img_table) != len(t_ccds): @@ -167,7 +179,7 @@ def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): Returns ------- imgs_dark : np.array - The dark backgrounds for the ACA images sampled from raw_dark_img. + The dark backgrounds in e-/s for the ACA images sampled from raw_dark_img. This will have the same length as imgrow0 and imgcol0, with shape (len(imgrow0), size, size). These pixels have not been scaled. """ diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index 4c8fa526..ecc6626c 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -532,19 +532,19 @@ def images_check_range(start, stop, img_table): assert np.allclose(img_dark, dark_row) -def test_get_aca_image_table_maude(): - """Test get_aca_image_table in the maude mode +def test_get_aca_images_maude(): + """Test get_aca_images in the maude mode This checks that the dark images are reasonable and the IMG data matches with and without the bgsub. """ tstart = "2021:001:00:00:00" tstop = "2021:001:00:01:00" - img_table_maude = chandra_aca.aca_image.get_aca_image_table( + img_table_maude = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="maude", bgsub=False ) images_check_range(tstart, tstop, img_table_maude) - img_table_maude_bgsub = chandra_aca.aca_image.get_aca_image_table( + img_table_maude_bgsub = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="maude", bgsub=True ) images_check_range(tstart, tstop, img_table_maude_bgsub) @@ -556,28 +556,28 @@ def test_get_aca_image_table_maude(): @pytest.mark.skipif(not HAS_ACA0_ARCHIVE, reason="No ACA0 archive") -def test_get_aca_image_table_mica(): - """Test get_aca_image_table in the mica mode +def test_get_aca_images_mica(): + """Test get_aca_images in the mica mode This checks that the dark images are reasonable and the answers match maude. """ tstart = "2012:180:00:00:00" tstop = "2012:180:00:01:00" - img_table_mica = chandra_aca.aca_image.get_aca_image_table( + img_table_mica = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="mica", bgsub=False ) images_check_range(tstart, tstop, img_table_mica) - img_table_mica_bgsub = chandra_aca.aca_image.get_aca_image_table( + img_table_mica_bgsub = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="mica", bgsub=True ) images_check_range(tstart, tstop, img_table_mica_bgsub) # Get maude data too for comparison - img_table_maude = chandra_aca.aca_image.get_aca_image_table( + img_table_maude = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="maude", bgsub=False ) - img_table_maude_bgsub = chandra_aca.aca_image.get_aca_image_table( + img_table_maude_bgsub = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="maude", bgsub=True ) diff --git a/chandra_aca/tests/test_dark_subtract.py b/chandra_aca/tests/test_dark_subtract.py index 3a703aae..414d9f02 100644 --- a/chandra_aca/tests/test_dark_subtract.py +++ b/chandra_aca/tests/test_dark_subtract.py @@ -6,7 +6,7 @@ from mica.archive.aca_dark import get_dark_cal_props from chandra_aca import maude_decom -from chandra_aca.aca_image import get_aca_image_table +from chandra_aca.aca_image import get_aca_images from chandra_aca.dark_subtract import ( get_aca_images_bgd_sub, get_dark_backgrounds, @@ -284,7 +284,7 @@ def test_dcsub_aca_images_maude(): imgs_table, img_dark=img_dark, tccd_dark=tccd_dark, t_ccd_vals=t_ccds ) - imgs_bgsub_table = get_aca_image_table(tstart, tstop, source="maude") + imgs_bgsub_table = get_aca_images(tstart, tstop, source="maude") exp0 = np.array( [ diff --git a/docs/dark_subtract.rst b/docs/dark_subtract.rst new file mode 100644 index 00000000..c5ade0b6 --- /dev/null +++ b/docs/dark_subtract.rst @@ -0,0 +1,5 @@ +chandra_aca.dark_subtract +========================= + +.. automodule:: chandra_aca.dark_subtract + :members: diff --git a/docs/index.rst b/docs/index.rst index e1d653ae..db10b901 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -45,6 +45,15 @@ Dark model dark_model.rst + +Dark background subtraction +--------------------------- + +.. toctree:: + :maxdepth: 2 + + dark_subtract.rst + MAUDE Decom ----------- From 71eb86a2ab5df346055f82becc8a7d5881bd2b41 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 25 Nov 2024 13:49:10 -0500 Subject: [PATCH 16/28] Fix some mica -> cxc pieces --- chandra_aca/aca_image.py | 4 ++-- chandra_aca/tests/test_aca_image.py | 32 ++++++++++++------------- chandra_aca/tests/test_dark_subtract.py | 4 ++-- 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 8a4da1d4..5701399b 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -890,13 +890,13 @@ def get_aca_images( t_ccds = chandra_aca.dark_subtract.get_tccd_data( imgs_table["TIME"], source=source, **maude_kwargs ) - elif source == "mica": + elif source == "cxc": imgs_table = mica.archive.aca_l0.get_aca_images(start, stop) t_ccds = chandra_aca.dark_subtract.get_tccd_data( imgs_table["TIME"], source=source ) else: - raise ValueError(f"source must be 'maude' or 'mica', not {source}") + raise ValueError(f"source must be 'maude' or 'cxc', not {source}") # Do some unmasking until the maude_decom and aca_l0 interfaces are updated with # fewer masked columns diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index ecc6626c..a3fca986 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -556,22 +556,22 @@ def test_get_aca_images_maude(): @pytest.mark.skipif(not HAS_ACA0_ARCHIVE, reason="No ACA0 archive") -def test_get_aca_images_mica(): - """Test get_aca_images in the mica mode +def test_get_aca_images_cxc(): + """Test get_aca_images in the cxc mode This checks that the dark images are reasonable and the answers match maude. """ tstart = "2012:180:00:00:00" tstop = "2012:180:00:01:00" - img_table_mica = chandra_aca.aca_image.get_aca_images( - tstart, tstop, source="mica", bgsub=False + img_table_cxc = chandra_aca.aca_image.get_aca_images( + tstart, tstop, source="cxc", bgsub=False ) - images_check_range(tstart, tstop, img_table_mica) - img_table_mica_bgsub = chandra_aca.aca_image.get_aca_images( - tstart, tstop, source="mica", bgsub=True + images_check_range(tstart, tstop, img_table_cxc) + img_table_cxc_bgsub = chandra_aca.aca_image.get_aca_images( + tstart, tstop, source="cxc", bgsub=True ) - images_check_range(tstart, tstop, img_table_mica_bgsub) + images_check_range(tstart, tstop, img_table_cxc_bgsub) # Get maude data too for comparison img_table_maude = chandra_aca.aca_image.get_aca_images( @@ -582,16 +582,16 @@ def test_get_aca_images_mica(): ) # Check that the two mica tables are the same in the key columns - assert np.allclose(img_table_mica["IMG"], img_table_mica_bgsub["IMG"]) - assert np.allclose(img_table_mica["TIME"], img_table_mica_bgsub["TIME"]) + assert np.allclose(img_table_cxc["IMG"], img_table_cxc_bgsub["IMG"]) + assert np.allclose(img_table_cxc["TIME"], img_table_cxc_bgsub["TIME"]) # Check that the tables are the same img_table_maude.sort(["TIME", "IMGNUM"]) img_table_maude_bgsub.sort(["TIME", "IMGNUM"]) - img_table_mica.sort(["TIME", "IMGNUM"]) - img_table_mica_bgsub.sort(["TIME", "IMGNUM"]) + img_table_cxc.sort(["TIME", "IMGNUM"]) + img_table_cxc_bgsub.sort(["TIME", "IMGNUM"]) - assert np.allclose(img_table_maude["IMG"], img_table_mica["IMG"]) - assert np.allclose(img_table_maude["TIME"], img_table_mica["TIME"]) - assert np.allclose(img_table_maude_bgsub["IMG"], img_table_mica_bgsub["IMG"]) - assert np.allclose(img_table_maude_bgsub["TIME"], img_table_mica_bgsub["TIME"]) + assert np.allclose(img_table_maude["IMG"], img_table_cxc["IMG"]) + assert np.allclose(img_table_maude["TIME"], img_table_cxc["TIME"]) + assert np.allclose(img_table_maude_bgsub["IMG"], img_table_cxc_bgsub["IMG"]) + assert np.allclose(img_table_maude_bgsub["TIME"], img_table_cxc_bgsub["TIME"]) diff --git a/chandra_aca/tests/test_dark_subtract.py b/chandra_aca/tests/test_dark_subtract.py index 414d9f02..7901e45a 100644 --- a/chandra_aca/tests/test_dark_subtract.py +++ b/chandra_aca/tests/test_dark_subtract.py @@ -299,7 +299,7 @@ def test_dcsub_aca_images_maude(): ] ) np.allclose(exp0, imgs_bgsub[0][0].filled(0), atol=1, rtol=0) - np.allclose(exp0, imgs_bgsub_table["BGSUB"][0][0].filled(0), atol=1, rtol=0) + np.allclose(exp0, imgs_bgsub_table["IMG_BGSUB"][0][0].filled(0), atol=1, rtol=0) # slot 4 is 6x6 and masked exp4 = np.array( @@ -315,4 +315,4 @@ def test_dcsub_aca_images_maude(): ] ) np.allclose(exp4, imgs_bgsub[4][0].filled(0), atol=1, rtol=0) - np.allclose(exp4, imgs_bgsub_table["BGSUB"][4][0].filled(0), atol=1, rtol=0) + np.allclose(exp4, imgs_bgsub_table["IMG_BGSUB"][4][0].filled(0), atol=1, rtol=0) From b3f58b17d2e4acae22c164aa76c27681c728acc8 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 25 Nov 2024 14:15:32 -0500 Subject: [PATCH 17/28] Be more strict about bounds --- chandra_aca/aca_image.py | 3 +++ chandra_aca/dark_subtract.py | 10 ++++++++-- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 5701399b..b0068ec5 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -858,6 +858,9 @@ def get_aca_images( - 'IMG_DARK': dark current image - 'T_CCD_SMOOTH': smoothed CCD temperature + The IMG_DARK individual values are only calculated if within the 1024x1024 + dark current map, otherwise they are set to 0. In practice this is not an + issue in that IMG and IMG_BGSUB must be within the CCD to be tracked. Parameters ---------- diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index bfe84413..b9f326c6 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -181,14 +181,20 @@ def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): imgs_dark : np.array The dark backgrounds in e-/s for the ACA images sampled from raw_dark_img. This will have the same length as imgrow0 and imgcol0, with shape - (len(imgrow0), size, size). These pixels have not been scaled. + (len(imgrow0), size, size). These pixels have not been scaled. Pixels + outside the raw_dark_img are set to 0. """ # Borrowed from the agasc code @numba.jit(nopython=True) def staggered_aca_slice(array_in, array_out, row, col): for i in np.arange(len(row)): - if row[i] + size < 1024 and col[i] + size < 1024: + if ( + row[i] >= 0 + and row[i] + size < 1024 + and col[i] >= 0 + and col[i] + size < 1024 + ): array_out[i] = array_in[row[i] : row[i] + size, col[i] : col[i] + size] imgs_dark = np.zeros([len(imgrow0), size, size], dtype=np.float64) From 81bc8f9e10be10ce2c5365cd4ad68be7cf078f8f Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 27 Nov 2024 10:57:18 -0500 Subject: [PATCH 18/28] Remove unmasking as no longer needed --- chandra_aca/aca_image.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index b0068ec5..fb58554e 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -901,16 +901,6 @@ def get_aca_images( else: raise ValueError(f"source must be 'maude' or 'cxc', not {source}") - # Do some unmasking until the maude_decom and aca_l0 interfaces are updated with - # fewer masked columns - imgs_table_unmasked = Table() - for col in imgs_table.colnames: - if col in ["TIME", "INTEG", "IMGROW0_8X8", "IMGCOL0_8X8"]: - imgs_table_unmasked[col] = imgs_table[col].data.data - else: - imgs_table_unmasked[col] = imgs_table[col] - imgs_table = imgs_table_unmasked - # Get background subtracted values if bgsub is True if bgsub: dark_data = mica.archive.aca_dark.get_dark_cal_props( From 5cba36ba775bdc9f39aac8fc91f378b76d7a4c37 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 27 Nov 2024 12:24:15 -0500 Subject: [PATCH 19/28] Wrong type in docstring --- chandra_aca/aca_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index fb58554e..4e388a8f 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -850,7 +850,7 @@ def get_aca_images( --------------------- ------- ----------- IMG_BGSUB float64 DN IMG_DARK float64 DN - T_CCD_SMOOTH float32 degC + T_CCD_SMOOTH float64 degC where: From 6e20ccf80970ac9dad70a01e4008a5381d108d0c Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 27 Nov 2024 12:24:59 -0500 Subject: [PATCH 20/28] Add some handling for zero length img_table --- chandra_aca/aca_image.py | 12 ++++++++++-- chandra_aca/dark_subtract.py | 4 ++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 4e388a8f..4e16ecbc 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -901,8 +901,10 @@ def get_aca_images( else: raise ValueError(f"source must be 'maude' or 'cxc', not {source}") - # Get background subtracted values if bgsub is True - if bgsub: + # Get background subtracted values if bgsub is True. + # There's nothing to do if there are no images (len(imgs_table) == 0), + # so special case that. + if bgsub and len(imgs_table) > 0: dark_data = mica.archive.aca_dark.get_dark_cal_props( imgs_table["TIME"].min(), select="nearest", @@ -921,4 +923,10 @@ def get_aca_images( imgs_table["IMG_DARK"] = imgs_dark imgs_table["T_CCD_SMOOTH"] = t_ccds + if bgsub and len(imgs_table) == 0: + # Add the columns to the table even if there are no rows + imgs_table["IMG_BGSUB"] = np.zeros(0) + imgs_table["IMG_DARK"] = np.zeros(0) + imgs_table["T_CCD_SMOOTH"] = np.zeros(0) + return imgs_table diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index b9f326c6..81f269fe 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -47,6 +47,10 @@ def get_tccd_data( CCD temperature values. """ + # If no times are given, return an empty array. + if len(times) == 0: + return np.array([]) + pad = 600 # 600 seconds default pad # increase the pad if the smooth window is large From 9b79677a6568836b1217d5def1b60ea8f2ce675e Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Mon, 2 Dec 2024 06:59:35 -0500 Subject: [PATCH 21/28] Change shapes for zero-length output --- chandra_aca/aca_image.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index 4e16ecbc..aff4fd78 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -925,8 +925,8 @@ def get_aca_images( if bgsub and len(imgs_table) == 0: # Add the columns to the table even if there are no rows - imgs_table["IMG_BGSUB"] = np.zeros(0) - imgs_table["IMG_DARK"] = np.zeros(0) - imgs_table["T_CCD_SMOOTH"] = np.zeros(0) + imgs_table["IMG_BGSUB"] = np.zeros(shape=(0, 8, 8)) + imgs_table["IMG_DARK"] = np.zeros(shape=(0, 8, 8)) + imgs_table["T_CCD_SMOOTH"] = np.zeros(shape=0) return imgs_table From 28625f40a73dcbe3fdaefca9d9d6b5b596106c63 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Mon, 2 Dec 2024 07:01:07 -0500 Subject: [PATCH 22/28] Update tests (now failing) --- chandra_aca/tests/test_aca_image.py | 52 +++++++++++------------------ 1 file changed, 20 insertions(+), 32 deletions(-) diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index a3fca986..abfb323e 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -491,7 +491,7 @@ def test_flicker_test_sequence(): ) -def images_check_range(start, stop, img_table): +def images_check_range(start, stop, img_table, *, bgsub): tstart = CxoTime(start).secs tstop = CxoTime(stop).secs @@ -506,11 +506,11 @@ def images_check_range(start, stop, img_table): assert np.all(np.isin(np.arange(8), img_table["IMGNUM"])) # Check that if the table has BGSUB column that IMG - DARK = BGSUB - if "BGSUB" in img_table.colnames: - assert np.allclose(img_table["IMG"] - img_table["DARK"], img_table["BGSUB"]) + if bgsub: + assert np.allclose(img_table["IMG"] - img_table["IMG_DARK"], img_table["IMG_BGSUB"]) # If there's a DARK column, then check it against the dark image from the dark cal - if "DARK" in img_table.colnames: + if bgsub: dc = get_dark_cal_props( tstart, select="nearest", include_image=True, aca_image=True @@ -518,8 +518,8 @@ def images_check_range(start, stop, img_table): full_img_dark = dc["image"] tccd_dark = dc["ccd_temp"] for row in img_table: - dark_row = row["DARK"] - tccd_row = row["T_CCD"] + dark_row = row["IMG_DARK"] + tccd_row = row["T_CCD_SMOOTH"] row8x8 = row["IMGROW0_8X8"] col8x8 = row["IMGCOL0_8X8"] dark_ref = full_img_dark.aca[row8x8 : row8x8 + 8, col8x8 : col8x8 + 8] @@ -532,54 +532,42 @@ def images_check_range(start, stop, img_table): assert np.allclose(img_dark, dark_row) -def test_get_aca_images_maude(): - """Test get_aca_images in the maude mode - - This checks that the dark images are reasonable and the IMG data matches with - and without the bgsub. - """ - tstart = "2021:001:00:00:00" - tstop = "2021:001:00:01:00" - img_table_maude = chandra_aca.aca_image.get_aca_images( - tstart, tstop, source="maude", bgsub=False - ) - images_check_range(tstart, tstop, img_table_maude) - img_table_maude_bgsub = chandra_aca.aca_image.get_aca_images( - tstart, tstop, source="maude", bgsub=True - ) - images_check_range(tstart, tstop, img_table_maude_bgsub) - assert np.allclose(img_table_maude["IMG"], img_table_maude_bgsub["IMG"]) - assert np.allclose(img_table_maude["TIME"], img_table_maude_bgsub["TIME"]) - - HAS_ACA0_ARCHIVE = (Path(mica.common.MICA_ARCHIVE) / "aca0").exists() @pytest.mark.skipif(not HAS_ACA0_ARCHIVE, reason="No ACA0 archive") -def test_get_aca_images_cxc(): +def test_get_aca_images_cxc_and_maude(): """Test get_aca_images in the cxc mode This checks that the dark images are reasonable and the answers match maude. """ - tstart = "2012:180:00:00:00" - tstop = "2012:180:00:01:00" + tstart = "2012:270:02:44:00" + tstop = "2012:270:02:47:00" + # Get CXC data and check that it looks reasonable img_table_cxc = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="cxc", bgsub=False ) - images_check_range(tstart, tstop, img_table_cxc) + images_check_range(tstart, tstop, img_table_cxc, bgsub=False) + img_table_cxc_bgsub = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="cxc", bgsub=True ) - images_check_range(tstart, tstop, img_table_cxc_bgsub) + images_check_range(tstart, tstop, img_table_cxc_bgsub, bgsub=True) - # Get maude data too for comparison + # Get MAUDE data and check that it looks reasonable img_table_maude = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="maude", bgsub=False ) + images_check_range(tstart, tstop, img_table_maude, bgsub=False) + img_table_maude_bgsub = chandra_aca.aca_image.get_aca_images( tstart, tstop, source="maude", bgsub=True ) + images_check_range(tstart, tstop, img_table_maude_bgsub, bgsub=True) + + assert np.allclose(img_table_maude["IMG"], img_table_maude_bgsub["IMG"]) + assert np.allclose(img_table_maude["TIME"], img_table_maude_bgsub["TIME"]) # Check that the two mica tables are the same in the key columns assert np.allclose(img_table_cxc["IMG"], img_table_cxc_bgsub["IMG"]) From 07d1c35e3ee63eb37d0dbf657a727d1bded8cb1d Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Mon, 2 Dec 2024 07:01:30 -0500 Subject: [PATCH 23/28] Docstring updates --- chandra_aca/aca_image.py | 15 +++++---- chandra_aca/dark_subtract.py | 65 ++++++++++++++++++------------------ 2 files changed, 40 insertions(+), 40 deletions(-) diff --git a/chandra_aca/aca_image.py b/chandra_aca/aca_image.py index aff4fd78..c236fd6f 100644 --- a/chandra_aca/aca_image.py +++ b/chandra_aca/aca_image.py @@ -839,7 +839,7 @@ def get_aca_images( start: CxoTimeLike, stop: CxoTimeLike, bgsub=True, source="maude", **maude_kwargs ) -> Table: """ - Get ACA images and ancillary data from either the maude or cxc data sources. + Get ACA images and ancillary data from either the MAUDE or CXC data sources. The returned table of ACA images and ancillary data will include the default columns returned by chandra_aca.maude_decom.get_aca_images or @@ -865,16 +865,17 @@ def get_aca_images( Parameters ---------- start : CxoTimeLike - start time + Start time. stop : CxoTimeLike - stop time (CXC sec) + Stop time (CXC sec). bgsub : bool - include background subtracted images in output table + Include background subtracted images in output table. source : str Data source for image and temperature telemetry ('maude' or 'cxc'). For 'cxc', - the image telemetry is from mica and temperature telemetry is from CXC L0 via cheta. - maude_kwargs - additional kwargs for maude data source + the image telemetry is from mica and temperature telemetry is from CXC L0 via + cheta. + **maude_kwargs + Additional kwargs for maude data source. Returns ------- diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index 81f269fe..b4eda771 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -28,18 +28,16 @@ def get_tccd_data( Parameters ---------- - times: np.array float of Chandra seconds - The times for which to get the CCD temperature data. + times: np.array (float) + Sampling times for CCD temperature data (CXC seconds). source : str, optional - Source of the CCD temperature data. If 'maude', override the fetch_sci data source - to 'maude allow_subset=False'. If 'cxc' use the local cxc archive. Default is 'maude'. + Source of CCD temperature data ('maude' (default) or 'cxc' for cheta archive). median_window : int, optional - Median filter window to remove outliers. Default is 3. + Median filter window to remove outliers (default=3). smooth_window : int, optional - Smooth the data using a hanning window of this length in samples. Default is 30. + Smooth data using a hanning window of this length in samples (default=30). maude_channel : str, optional - The maude channel to use. Default is None. - + Maude channel to use (default is flight). Returns ------- @@ -95,21 +93,21 @@ def get_aca_images_bgd_sub(img_table, t_ccd_vals, img_dark, tccd_dark): Parameters ---------- img_table : astropy.table.Table - The table of ACA images. + Table of ACA images with columns 'IMG', 'IMGROW0_8X8', 'IMGCOL0_8X8', 'INTEG'. t_ccd_vals : np.array - The CCD temperature values at the times of the ACA images in deg C. + CCD temperature values at the times of the ACA images (deg C). img_dark : np.array - The dark calibration image. Must be 1024x1024. In e-/s. + Dark calibration image. Must be 1024x1024 (e-/s). tccd_dark : float - The reference temperature of the dark calibration image in deg C. + Reference temperature of the dark calibration image (deg C). Returns ------- tuple (imgs_bgsub, imgs_dark) imgs_bgsub : np.array - The background subtracted ACA images in DN. + Background subtracted ACA images (DN). imgs_dark : np.array - The dark current images in DN. + Dark current images (DN). """ imgs_dark = get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccd_vals) imgs_bgsub = img_table["IMG"] - imgs_dark @@ -122,28 +120,31 @@ def get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccds): """ Get the scaled dark current values for a table of ACA images. + This scales the dark current to the appropriate temperature and integration time, + returning the dark current in DN matching the ACA images in img_table. + Parameters ---------- img_table : astropy.table.Table + Table of ACA images with columns 'IMG', 'IMGROW0_8X8', 'IMGCOL0_8X8', 'INTEG'. img_dark : 1024x1024 array - The dark calibration image. + Dark calibration image. tccd_dark : float - The reference temperature of the dark calibration image. + Reference temperature of the dark calibration image. t_ccds : array - The cheta temperature values at the times of img_table. + Cheta temperature values at the times of img_table. Returns ------- - imgs_dark : np.array - An array containing the temperature scaled dark current for each ACA image - in the img_table in DN. + imgs_dark : np.array (len(img_table), 8, 8) + Temperature scaled dark current for each ACA image in img_table (DN). """ if len(img_table) != len(t_ccds): raise ValueError("img_table and t_ccds must have the same length") if img_dark.shape != (1024, 1024): - raise ValueError("Dark image shape is not 1024x1024") + raise ValueError("dark image shape is not 1024x1024") imgs_dark_unscaled = get_dark_backgrounds( img_dark, @@ -167,26 +168,24 @@ def get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccds): def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): """ - Get the dark background for a stack/table of ACA image. + Get dark background cutouts at a set of ACA image positions. Parameters ---------- raw_dark_img : np.array - The dark calibration image. - imgrow0 : np.array - The row of the ACA image. - imgcol0 : np.array - The column of the ACA image. + Dark calibration image. + imgrow0 : np.array (int) + Row of ACA image. + imgcol0 : np.array (int) + Column of ACA image. size : int, optional - The size of the ACA image. Default is 8. + Size of ACA image (default=8). Returns ------- - imgs_dark : np.array - The dark backgrounds in e-/s for the ACA images sampled from raw_dark_img. - This will have the same length as imgrow0 and imgcol0, with shape - (len(imgrow0), size, size). These pixels have not been scaled. Pixels - outside the raw_dark_img are set to 0. + imgs_dark : np.array (len(imgrow0), size, size) + Dark backgrounds for image locations sampled from raw_dark_img (e-/s). + Pixels outside dark_img are set to 0.0. """ # Borrowed from the agasc code From 4a8a73ca644cec7de9e7a1dfd068211eafc131a6 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 2 Dec 2024 15:42:13 -0500 Subject: [PATCH 24/28] Update dark map extraction to be per-pixel to work near edges --- chandra_aca/dark_subtract.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index b4eda771..92bf5870 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -185,20 +185,20 @@ def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): ------- imgs_dark : np.array (len(imgrow0), size, size) Dark backgrounds for image locations sampled from raw_dark_img (e-/s). - Pixels outside dark_img are set to 0.0. + Pixels outside raw_dark_img are set to 0.0. """ + CCD_MAX = 1024 + CCD_MIN = 0 - # Borrowed from the agasc code @numba.jit(nopython=True) def staggered_aca_slice(array_in, array_out, row, col): - for i in np.arange(len(row)): - if ( - row[i] >= 0 - and row[i] + size < 1024 - and col[i] >= 0 - and col[i] + size < 1024 - ): - array_out[i] = array_in[row[i] : row[i] + size, col[i] : col[i] + size] + for idx in np.arange(len(row)): + row_out = np.zeros((size, size), dtype=np.float64) + for i, r in enumerate(np.arange(row[idx], row[idx] + size)): + for j, c in enumerate(np.arange(col[idx], col[idx] + size)): + if r >= CCD_MIN and r < CCD_MAX and c >= CCD_MIN and c < CCD_MAX: + row_out[i, j] = array_in[r, c] + array_out[idx] = row_out imgs_dark = np.zeros([len(imgrow0), size, size], dtype=np.float64) staggered_aca_slice( From 3ea97077b1c5991d7e9dc60be6eb80feea31b030 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 2 Dec 2024 15:43:17 -0500 Subject: [PATCH 25/28] Update test to work at positive ccd edges --- chandra_aca/tests/test_aca_image.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index abfb323e..a822b5fd 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -522,14 +522,26 @@ def images_check_range(start, stop, img_table, *, bgsub): tccd_row = row["T_CCD_SMOOTH"] row8x8 = row["IMGROW0_8X8"] col8x8 = row["IMGCOL0_8X8"] + + # If row8x8 or col8x8 is < 0 just skip the test for now + if row8x8 < 0 or col8x8 < 0: + continue + dark_ref = full_img_dark.aca[row8x8 : row8x8 + 8, col8x8 : col8x8 + 8] from chandra_aca.dark_model import dark_temp_scale_img dark_scale = dark_temp_scale_img(dark_ref, tccd_dark, tccd_row) # Default to using 1.696 integ time if not present in the image table integ = row["INTEG"] if "INTEG" in img_table.colnames else 1.696 - img_dark = dark_scale * integ / 5 - assert np.allclose(img_dark, dark_row) + img_dark = ACAImage(dark_scale * integ / 5, row0=row8x8, col0=col8x8) + + # Make an 8x8 expected image (even if the dark current goes off the CCD) + expected_img = ( + ACAImage(np.zeros((8, 8)), row0=row8x8, col0=col8x8).aca + img_dark + ) + + # Compare + assert np.allclose(expected_img, dark_row) HAS_ACA0_ARCHIVE = (Path(mica.common.MICA_ARCHIVE) / "aca0").exists() From 45e01087f7efb2603d73d96b758babe764358d19 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 2 Dec 2024 15:43:38 -0500 Subject: [PATCH 26/28] Pre-ruff formatting fix --- chandra_aca/tests/test_aca_image.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index a822b5fd..e824e372 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -507,11 +507,12 @@ def images_check_range(start, stop, img_table, *, bgsub): # Check that if the table has BGSUB column that IMG - DARK = BGSUB if bgsub: - assert np.allclose(img_table["IMG"] - img_table["IMG_DARK"], img_table["IMG_BGSUB"]) + assert np.allclose( + img_table["IMG"] - img_table["IMG_DARK"], img_table["IMG_BGSUB"] + ) # If there's a DARK column, then check it against the dark image from the dark cal if bgsub: - dc = get_dark_cal_props( tstart, select="nearest", include_image=True, aca_image=True ) From 4c876652caab8949bc5916da03c6a9eeddda6ae4 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Mon, 2 Dec 2024 16:03:56 -0500 Subject: [PATCH 27/28] Oops, edge at -512 --- chandra_aca/tests/test_aca_image.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/chandra_aca/tests/test_aca_image.py b/chandra_aca/tests/test_aca_image.py index e824e372..fc42218f 100644 --- a/chandra_aca/tests/test_aca_image.py +++ b/chandra_aca/tests/test_aca_image.py @@ -524,8 +524,8 @@ def images_check_range(start, stop, img_table, *, bgsub): row8x8 = row["IMGROW0_8X8"] col8x8 = row["IMGCOL0_8X8"] - # If row8x8 or col8x8 is < 0 just skip the test for now - if row8x8 < 0 or col8x8 < 0: + # If row8x8 or col8x8 is < -512 just skip the test for now + if row8x8 < -512 or col8x8 < -512: continue dark_ref = full_img_dark.aca[row8x8 : row8x8 + 8, col8x8 : col8x8 + 8] From 37560e590d0990e2546b749f14c63fee830bfbd5 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Tue, 3 Dec 2024 07:24:29 -0500 Subject: [PATCH 28/28] Move the numba function to module scope This improves performance substantially. --- chandra_aca/dark_subtract.py | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/chandra_aca/dark_subtract.py b/chandra_aca/dark_subtract.py index 92bf5870..734f79d9 100644 --- a/chandra_aca/dark_subtract.py +++ b/chandra_aca/dark_subtract.py @@ -166,6 +166,24 @@ def get_dark_current_imgs(img_table, img_dark, tccd_dark, t_ccds): return imgs_dark +CCD_MAX = 1024 +CCD_MIN = 0 + + +@numba.jit(nopython=True) +def staggered_aca_slice(array_in, array_out, row, col): + """Make cutouts of array_in at positions row, col and put them in array_out. + + array_out must be 3D with shape (N, sz_r, sz_c) and be initialized to 0. + row and col must be 1D arrays of length N. + """ + for idx in np.arange(array_out.shape[0]): + for i, r in enumerate(np.arange(row[idx], row[idx] + array_out.shape[1])): + for j, c in enumerate(np.arange(col[idx], col[idx] + array_out.shape[2])): + if r >= CCD_MIN and r < CCD_MAX and c >= CCD_MIN and c < CCD_MAX: + array_out[idx, i, j] = array_in[r, c] + + def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): """ Get dark background cutouts at a set of ACA image positions. @@ -187,19 +205,6 @@ def get_dark_backgrounds(raw_dark_img, imgrow0, imgcol0, size=8): Dark backgrounds for image locations sampled from raw_dark_img (e-/s). Pixels outside raw_dark_img are set to 0.0. """ - CCD_MAX = 1024 - CCD_MIN = 0 - - @numba.jit(nopython=True) - def staggered_aca_slice(array_in, array_out, row, col): - for idx in np.arange(len(row)): - row_out = np.zeros((size, size), dtype=np.float64) - for i, r in enumerate(np.arange(row[idx], row[idx] + size)): - for j, c in enumerate(np.arange(col[idx], col[idx] + size)): - if r >= CCD_MIN and r < CCD_MAX and c >= CCD_MIN and c < CCD_MAX: - row_out[i, j] = array_in[r, c] - array_out[idx] = row_out - imgs_dark = np.zeros([len(imgrow0), size, size], dtype=np.float64) staggered_aca_slice( raw_dark_img.astype(float), imgs_dark, 512 + imgrow0, 512 + imgcol0