Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

JP-3560: Updating CHARGELOSS Flagging #8336

Merged
merged 3 commits into from
Mar 11, 2024
Merged

JP-3560: Updating CHARGELOSS Flagging #8336

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1820072
Updating the CHARGELOSS flagging to flag all later groups in a ramp. …
kmacdonald-stsci Mar 8, 2024
de5d2b6
Updating the change log.
kmacdonald-stsci Mar 8, 2024
e66abdc
Updating the change log with a more complete description of the PR.
kmacdonald-stsci Mar 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions CHANGES.rst
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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