Skip to content

Commit

Permalink
JP-3560: Updating CHARGELOSS Flagging (#8336)
Browse files Browse the repository at this point in the history
  • Loading branch information
@kmacdonald-stsci
kmacdonald-stsci authored Mar 11, 2024
1 parent cf166a0 commit bf3b06f
Show file tree
Hide file tree
Showing 3 changed files with 164 additions and 154 deletions.
9 changes: 9 additions & 0 deletions CHANGES.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,15 @@ background
- Updated to allow multi-integration (rateints) background exposures to have
a different value of NINTS than the science exposure. [#8326]

charge_migration
----------------

- Updated the CHARGELOSS flagging. In an integration ramp, the first group in
the SCI data is found that is above the CHARGELOSS threshold and not flagged
as DO_NOT_USE. This group, and all subsequent groups, are then flagged as
CHARGELOSS and DO_NOT_USE. The four nearest pixel neighbor are then flagged
in the same group. [#8336]

cube_build
----------

Expand Down
92 changes: 25 additions & 67 deletions jwst/charge_migration/charge_migration.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@

from stdatamodels.jwst.datamodels import dqflags


log = logging.getLogger(__name__)
log.setLevel(logging.DEBUG)

Expand Down Expand Up @@ -72,72 +73,29 @@ def flag_pixels(data, gdq, signal_threshold):
updated group dq array
"""
n_ints, n_grps, n_rows, n_cols = gdq.shape

new_gdq = gdq.copy() # Updated gdq

# Flag all exceedances with CHARGELOSS and NO_NOT_USE
chargeloss_pix = (data > signal_threshold) & (gdq != DNU)
new_gdq[chargeloss_pix] = np.bitwise_or(new_gdq[chargeloss_pix], CHLO | DNU)

# Reset groups previously flagged as DNU
gdq_orig = gdq.copy() # For resetting to previously flagged DNU
wh_gdq_DNU = np.bitwise_and(gdq_orig, DNU)

# Get indices for exceedances
arg_where = np.argwhere(new_gdq == CHLO_DNU)

a_int = arg_where[:, 0] # array of integrations
a_grp = arg_where[:, 1] # array of groups
a_row = arg_where[:, 2] # array of rows
a_col = arg_where[:, 3] # array of columns

# Process the 4 nearest neighbors of each exceedance
# Pixel to the east
xx_max_p1 = a_col[a_col < (n_cols-1)] + 1
i_int = a_int[a_col < (n_cols-1)]
i_grp = a_grp[a_col < (n_cols-1)]
i_row = a_row[a_col < (n_cols-1)]

if len(xx_max_p1) > 0:
new_gdq[i_int, i_grp, i_row, xx_max_p1] = \
np.bitwise_or(new_gdq[i_int, i_grp, i_row, xx_max_p1], CHLO | DNU)

new_gdq[wh_gdq_DNU == 1] = gdq_orig[wh_gdq_DNU == 1] # reset for earlier DNUs

# Pixel to the west
xx_m1 = a_col[a_col > 0] - 1
i_int = a_int[a_col > 0]
i_grp = a_grp[a_col > 0]
i_row = a_row[a_col > 0]

if len(xx_m1) > 0:
new_gdq[i_int, i_grp, i_row, xx_m1] = \
np.bitwise_or(new_gdq[i_int, i_grp, i_row, xx_m1], CHLO | DNU)

new_gdq[wh_gdq_DNU == 1] = gdq_orig[wh_gdq_DNU == 1] # reset for earlier DNUs

# Pixel to the north
yy_m1 = a_row[a_row > 0] - 1
i_int = a_int[a_row > 0]
i_grp = a_grp[a_row > 0]
i_col = a_col[a_row > 0]

if len(yy_m1) > 0:
new_gdq[i_int, i_grp, yy_m1, i_col] = \
np.bitwise_or(new_gdq[i_int, i_grp, yy_m1, i_col], CHLO | DNU)

new_gdq[wh_gdq_DNU == 1] = gdq_orig[wh_gdq_DNU == 1] # reset for earlier DNUs

# Pixel to the south
yy_max_p1 = a_row[a_row < (n_rows-1)] + 1
i_int = a_int[a_row < (n_rows-1)]
i_grp = a_grp[a_row < (n_rows-1)]
i_col = a_col[a_row < (n_rows-1)]

if len(yy_max_p1) > 0:
new_gdq[i_int, i_grp, yy_max_p1, i_col] = \
np.bitwise_or(new_gdq[i_int, i_grp, yy_max_p1, i_col], CHLO | DNU)

new_gdq[wh_gdq_DNU == 1] = gdq_orig[wh_gdq_DNU == 1] # reset for earlier DNUs
chargeloss_pix = np.where((data > signal_threshold) & (gdq != DNU))

new_gdq = gdq.copy()

for k in range(len(chargeloss_pix[0])):
integ, group = chargeloss_pix[0][k], chargeloss_pix[1][k]
row, col = chargeloss_pix[2][k], chargeloss_pix[3][k]
new_gdq[integ, group:, row, col] |= CHLO_DNU

# North
if row > 0:
new_gdq[integ, group:, row-1, col] |= CHLO_DNU

# South
if row < (n_rows-1):
new_gdq[integ, group:, row+1, col] |= CHLO_DNU

# East
if col < (n_cols-1):
new_gdq[integ, group:, row, col+1] |= CHLO_DNU

# West
if col > 0:
new_gdq[integ, group:, row, col-1] |= CHLO_DNU

return new_gdq
217 changes: 130 additions & 87 deletions jwst/charge_migration/tests/test_charge_migration.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@

import numpy.testing as npt


test_dq_flags = dqflags.pixel
GOOD = test_dq_flags["GOOD"]
DNU = test_dq_flags["DO_NOT_USE"]
Expand Down Expand Up @@ -78,120 +79,162 @@ def test_pix_1():
true_out_gdq[0, 4:, 0, 0] = CHLO_DNU

out_model = charge_migration(ramp_model, signal_threshold)

out_data = out_model.data
out_gdq = out_model.groupdq

npt.assert_array_equal(true_out_data, out_data)
npt.assert_array_equal(out_gdq, true_out_gdq)


def test_too_few_groups():
def test_pix_2():
"""
Test that processing for datasets having too few (<3) groups per integration
are skipped.
Test a later group being below the threshold.
"""
ngroups, nints, nrows, ncols = 2, 1, 1, 1
ngroups, nints, nrows, ncols = 10, 1, 1, 1
ramp_model, pixdq, groupdq, err = create_mod_arrays(
ngroups, nints, nrows, ncols)

ramp_model.data[0, :, 0, 0] = 20000.
sig_thresh = 100.
signal_threshold = 4000.

result = ChargeMigrationStep.call(ramp_model, skip=False,
signal_threshold=sig_thresh)
status = result.meta.cal_step.charge_migration
arr = [1000., 2000., 4005., 4500., 5000., 5500., 3500., 6000., 6500., 3700.]
ramp_model.data[0, :, 0, 0] = np.array(arr, dtype=np.float32)
arr = [0, DNU, 0, 0, 0, 0, 0, 0, 0, 0]
ramp_model.groupdq[0, :, 0, 0] = np.array(arr, dtype=np.uint8)

npt.assert_string_equal(status, "SKIPPED")
out_model = charge_migration(ramp_model, signal_threshold)

truth_arr = [0, DNU, CHLO_DNU, CHLO_DNU, CHLO_DNU, CHLO_DNU, CHLO_DNU, CHLO_DNU, CHLO_DNU, CHLO_DNU]
truth_gdq = np.array(truth_arr, dtype=np.uint8)

def test_flag_neighbors():
npt.assert_array_equal(truth_gdq, out_model.groupdq[0, :, 0, 0])



def nearest_neighbor_base(chg_thresh, pixel):
"""
Test flagging of 4 nearest neighbors of exceedances. Tests pixels on
array edges, Tests exclusion of groups previously flagged as DO_NOT_USE.
Set up ramp array that is 5, 5 with 10 groups.
The flagging starts in group 3 (zero based) in the pixel tested.
"""
ngroups, nints, nrows, ncols = 6, 1, 4, 3
nints, ngroups, nrows, ncols = 1, 10, 5, 5
ramp_model, pixdq, groupdq, err = create_mod_arrays(
ngroups, nints, nrows, ncols)

signal_threshold = 4400.
# Set up dummy data
base = chg_thresh * 0.05
base_arr = [float(k+1) * base for k in range(ngroups)]
for row in range(nrows):
for col in range(ncols):
ramp_model.data[0, :, row, col] = np.array(base_arr, dtype=np.float32)

# Populate pixel-specific SCI and GROUPDQ arrays
ramp_model.data[0, :, :, :] = \
np.array([[
[1900., 2666., 2100.],
[3865., 2300., 3177.],
[3832., 3044., 3588.],
[3799., 3233., 3000.]],

[[2100., 2866., 2300.],
[4065., 2500., 3377.],
[4032., 3244., 3788.],
[3999., 3433., 3200.]],

[[2300., 3066., 2500.],
[4265., 2700., 3577.],
[4232., 3444., 3988.],
[4199., 3633., 3400.]],

[[2500., 3266., 2700.],
[4465., 2900., 3777.],
[4432., 3644., 4188.],
[4399., 3833., 3600.]],

[[2700., 3466., 2900.],
[4665., 3100., 3977.],
[4632., 3844., 4388.],
[4599., 4033., 3800.]],

[[2900., 3666., 3100.],
[4865., 3300., 4177.],
[4832., 4044., 4588.],
[4799., 4233., 4000.]]], dtype=np.float32)

# These group DQ values should propagate unchanged to the output
ramp_model.groupdq[:, 4, 2, 0] = [DNU]
ramp_model.groupdq[:, 1, 2, 2] = [DNU]
ramp_model.groupdq[:, 2, 1, 1] = [DROU + DNU]
# Make CHARGELOSS threshold starting at group 3
in_row, in_col = pixel
ramp_model.data[0, 3:, in_row, in_col] += chg_thresh

out_model = charge_migration(ramp_model, signal_threshold)
return ramp_model, pixdq, groupdq, err

out_gdq = out_model.groupdq

true_out_gdq = ramp_model.groupdq.copy()
true_out_gdq[0, :, :, :] = \
np.array([[
[0, 0, 0],
[0, 0, 0],
[0, 0, 0],
[0, 0, 0]],

[[0, 0, 0],
[0, 0, 0],
[0, 0, DNU],
[0, 0, 0]],

[[0, 0, 0],
[0, 9, 0],
[0, 0, 0],
[0, 0, 0]],

[[CHLO_DNU, 0, 0],
[CHLO_DNU, CHLO_DNU, 0],
[CHLO_DNU, CHLO_DNU, 0],
[CHLO_DNU, 0, 0]],

[[CHLO_DNU, 0, 0],
[CHLO_DNU, CHLO_DNU, 0],
[DNU, 0, 0],
[CHLO_DNU, CHLO_DNU, 0]],

[[CHLO_DNU, 0, 0],
[CHLO_DNU, CHLO_DNU, CHLO_DNU],
[CHLO_DNU, CHLO_DNU, CHLO_DNU],
[CHLO_DNU, CHLO_DNU, CHLO_DNU]]], dtype=np.uint8)
def test_nearest_neighbor_1():
"""
CHARGELOSS center
The flagging starts in group 3 (zero based) in the pixel tested.
"""
chg_thresh = 4000.
pixel = (2, 2)
ramp_model, pixdq, groupdq, err = nearest_neighbor_base(chg_thresh, pixel)
gdq_check = ramp_model.groupdq.copy()
ngroups = gdq_check.shape[1]

out_model = charge_migration(ramp_model, chg_thresh)

check_pattern = [
[GOOD, GOOD, GOOD, GOOD, GOOD],
[GOOD, GOOD, CHLO_DNU, GOOD, GOOD],
[GOOD, CHLO_DNU, CHLO_DNU, CHLO_DNU, GOOD],
[GOOD, GOOD, CHLO_DNU, GOOD, GOOD],
[GOOD, GOOD, GOOD, GOOD, GOOD],
]
check = np.array(check_pattern, dtype=gdq_check.dtype)
for group in range(3, ngroups):
gdq_check[0, group, :, :] = check

npt.assert_array_equal(out_model.data, ramp_model.data)
npt.assert_array_equal(out_gdq, true_out_gdq)
npt.assert_array_equal(out_model.groupdq, gdq_check)


def test_nearest_neighbor_2():
"""
CHARGELOSS corner
The flagging starts in group 3 (zero based) in the pixel tested.
"""
chg_thresh = 4000.
pixel = (0, 0)
ramp_model, pixdq, groupdq, err = nearest_neighbor_base(chg_thresh, pixel)
gdq_check = ramp_model.groupdq.copy()
ngroups = gdq_check.shape[1]

out_model = charge_migration(ramp_model, chg_thresh)

check_pattern = [
[CHLO_DNU, CHLO_DNU, GOOD, GOOD, GOOD],
[CHLO_DNU, GOOD, GOOD, GOOD, GOOD],
[GOOD, GOOD, GOOD, GOOD, GOOD],
[GOOD, GOOD, GOOD, GOOD, GOOD],
[GOOD, GOOD, GOOD, GOOD, GOOD],
]
check = np.array(check_pattern, dtype=gdq_check.dtype)
for group in range(3, ngroups):
gdq_check[0, group, :, :] = check

npt.assert_array_equal(out_model.data, ramp_model.data)
npt.assert_array_equal(out_model.groupdq, gdq_check)


def test_nearest_neighbor_3():
"""
CHARGELOSS Edge
The flagging starts in group 3 (zero based) in the pixel tested.
"""
chg_thresh = 4000.
pixel = (2, 4)
ramp_model, pixdq, groupdq, err = nearest_neighbor_base(chg_thresh, pixel)
gdq_check = ramp_model.groupdq.copy()
ngroups = gdq_check.shape[1]

out_model = charge_migration(ramp_model, chg_thresh)

check_pattern = [
[GOOD, GOOD, GOOD, GOOD, GOOD],
[GOOD, GOOD, GOOD, GOOD, CHLO_DNU],
[GOOD, GOOD, GOOD, CHLO_DNU, CHLO_DNU],
[GOOD, GOOD, GOOD, GOOD, CHLO_DNU],
[GOOD, GOOD, GOOD, GOOD, GOOD],
]
check = np.array(check_pattern, dtype=gdq_check.dtype)
for group in range(3, ngroups):
gdq_check[0, group, :, :] = check

npt.assert_array_equal(out_model.data, ramp_model.data)
npt.assert_array_equal(out_model.groupdq, gdq_check)


def test_too_few_groups():
"""
Test that processing for datasets having too few (<3) groups per integration
are skipped.
"""
ngroups, nints, nrows, ncols = 2, 1, 1, 1
ramp_model, pixdq, groupdq, err = create_mod_arrays(
ngroups, nints, nrows, ncols)

ramp_model.data[0, :, 0, 0] = 20000.
sig_thresh = 100.

result = ChargeMigrationStep.call(ramp_model, skip=False,
signal_threshold=sig_thresh)
status = result.meta.cal_step.charge_migration

npt.assert_string_equal(status, "SKIPPED")


def create_mod_arrays(ngroups, nints, nrows, ncols):
Expand Down

0 comments on commit bf3b06f

Please sign in to comment.