Skip to content

Commit

Permalink
JP-2693: Flag groups do_not_use after detected ramp jumps (#110)
Browse files Browse the repository at this point in the history 
* initial flag n groups after jump

* adding log messages

* updating to 2 thresholds and 2 # of groups

* parameter passing works and adding DN thresholds

* after_jump flagging is working

* fully working!

* bug in the non-working code fixed

* moving after_flags parameters to keywords and adding tests

* adjusting for jumps in electrons not DNs

* fix too long line

* adding to change log

* changes in right section now

* fixing keywords

* fixing tests for revised keywords

* responding to PR comments

* Address in PR comments
  • Loading branch information
@karllark
karllark committed Aug 17, 2022
1 parent 4bbffcc commit 73a88c5
Show file tree
Hide file tree
Showing 4 changed files with 208 additions and 8 deletions.
6 changes: 5 additions & 1 deletion CHANGES.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,11 @@ ramp_fitting
Changes to API
--------------

jump
~~~~

- Added flagging after detected ramp jumps based on two DN thresholds and
two number of groups to flag [#110]

1.0.0 (2022-06-24)
==================
Expand Down Expand Up @@ -144,7 +148,7 @@ saturation
- Pin astropy min version to 5.0.4. [#82]

- Fix for jumps in first good group after dropping groups [#84]


0.6.2 (22-03-29)
================
Expand Down
39 changes: 35 additions & 4 deletions src/stcal/jump/jump.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,11 @@
def detect_jumps(frames_per_group, data, gdq, pdq, err,
gain_2d, readnoise_2d, rejection_thresh,
three_grp_thresh, four_grp_thresh, max_cores, max_jump_to_flag_neighbors,
min_jump_to_flag_neighbors, flag_4_neighbors, dqflags):
min_jump_to_flag_neighbors, flag_4_neighbors, dqflags,
after_jump_flag_dn1=0.0,
after_jump_flag_n1=0,
after_jump_flag_dn2=0.0,
after_jump_flag_n2=0):
"""
This is the high-level controlling routine for the jump detection process.
It loads and sets the various input data and parameters needed by each of
Expand Down Expand Up @@ -88,7 +92,21 @@ def detect_jumps(frames_per_group, data, gdq, pdq, err,
A dictionary with at least the following keywords:
DO_NOT_USE, SATURATED, JUMP_DET, NO_GAIN_VALUE, GOOD
after_jump_flag_dn1 : float
Jumps with amplitudes above the specified DN value will have subsequent
groups flagged with the number determined by the after_jump_flag_n1
after_jump_flag_n1 : int
Gives the number of groups to flag after jumps with DN values above that
given by after_jump_flag_dn1
after_jump_flag_dn2 : float
Jumps with amplitudes above the specified DN value will have subsequent
groups flagged with the number determined by the after_jump_flag_n2
after_jump_flag_n2 : int
Gives the number of groups to flag after jumps with DN values above that
given by after_jump_flag_dn2
Returns
-------
Expand Down Expand Up @@ -116,6 +134,9 @@ def detect_jumps(frames_per_group, data, gdq, pdq, err,
data *= gain_2d
err *= gain_2d
readnoise_2d *= gain_2d
# also apply to the after_jump thresholds
after_jump_flag_e1 = after_jump_flag_dn1 * gain_2d
after_jump_flag_e2 = after_jump_flag_dn2 * gain_2d

# Apply the 2-point difference method as a first pass
log.info('Executing two-point difference method')
Expand Down Expand Up @@ -149,7 +170,11 @@ def detect_jumps(frames_per_group, data, gdq, pdq, err,
twopt.find_crs(data, gdq, readnoise_2d, rejection_thresh,
three_grp_thresh, four_grp_thresh, frames_per_group,
flag_4_neighbors, max_jump_to_flag_neighbors,
min_jump_to_flag_neighbors, dqflags)
min_jump_to_flag_neighbors, dqflags,
after_jump_flag_e1=after_jump_flag_e1,
after_jump_flag_n1=after_jump_flag_n1,
after_jump_flag_e2=after_jump_flag_e2,
after_jump_flag_n2=after_jump_flag_n2)

elapsed = time.time() - start
else:
Expand All @@ -174,7 +199,10 @@ def detect_jumps(frames_per_group, data, gdq, pdq, err,
rejection_thresh, three_grp_thresh, four_grp_thresh,
frames_per_group, flag_4_neighbors,
max_jump_to_flag_neighbors,
min_jump_to_flag_neighbors, dqflags, copy_arrs))
min_jump_to_flag_neighbors, dqflags,
after_jump_flag_e1, after_jump_flag_n1,
after_jump_flag_e2, after_jump_flag_n2,
copy_arrs))

# last slice get the rest
slices.insert(n_slices - 1, (data[:, :, (n_slices - 1) * yinc:n_rows, :],
Expand All @@ -183,7 +211,10 @@ def detect_jumps(frames_per_group, data, gdq, pdq, err,
rejection_thresh, three_grp_thresh,
four_grp_thresh, frames_per_group,
flag_4_neighbors, max_jump_to_flag_neighbors,
min_jump_to_flag_neighbors, dqflags, copy_arrs))
min_jump_to_flag_neighbors, dqflags,
after_jump_flag_e1, after_jump_flag_n1,
after_jump_flag_e2, after_jump_flag_n2,
copy_arrs))
log.info("Creating %d processes for jump detection " % n_slices)
pool = multiprocessing.Pool(processes=n_slices)
# Starts each slice in its own process. Starmap allows more than one
Expand Down
55 changes: 52 additions & 3 deletions src/stcal/jump/twopoint_difference.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,12 @@
def find_crs(dataa, group_dq, read_noise, rejection_thresh,
two_diff_rej_thresh, three_diff_rej_thresh, nframes,
flag_4_neighbors, max_jump_to_flag_neighbors,
min_jump_to_flag_neighbors, dqflags, copy_arrs=True):
min_jump_to_flag_neighbors, dqflags,
after_jump_flag_e1=0.0,
after_jump_flag_n1=0,
after_jump_flag_e2=0.0,
after_jump_flag_n2=0,
copy_arrs=True):

"""
Find CRs/Jumps in each integration within the input data array. The input
Expand Down Expand Up @@ -54,6 +59,26 @@ def find_crs(dataa, group_dq, read_noise, rejection_thresh,
neighbors (marginal detections). Any primary jump below this value will
not have its neighbors flagged.
dqflags: dict
A dictionary with at least the following keywords:
DO_NOT_USE, SATURATED, JUMP_DET, NO_GAIN_VALUE, GOOD
after_jump_flag_e1 : float
Jumps with amplitudes above the specified e value will have subsequent
groups flagged with the number determined by the after_jump_flag_n1
after_jump_flag_n1 : int
Gives the number of groups to flag after jumps with DN values above that
given by after_jump_flag_dn1
after_jump_flag_e2 : float
Jumps with amplitudes above the specified e value will have subsequent
groups flagged with the number determined by the after_jump_flag_n2
after_jump_flag_n2 : int
Gives the number of groups to flag after jumps with DN values above that
given by after_jump_flag_dn2
copy_arrs : bool
Flag for making internal copies of the arrays so the input isn't modified,
defaults to True.
Expand Down Expand Up @@ -122,7 +147,8 @@ def find_crs(dataa, group_dq, read_noise, rejection_thresh,
# compute 'ratio' for each group. this is the value that will be
# compared to 'threshold' to classify jumps. subtract the median of
# first_diffs from first_diffs, take the abs. value and divide by sigma.
ratio = np.abs(first_diffs - median_diffs[np.newaxis, :, :]) / sigma[np.newaxis, :, :]
e_jump = first_diffs - median_diffs[np.newaxis, :, :]
ratio = np.abs(e_jump) / sigma[np.newaxis, :, :]

# create a 2d array containing the value of the largest 'ratio' for each group
max_ratio = np.nanmax(ratio, axis=0)
Expand All @@ -138,7 +164,7 @@ def find_crs(dataa, group_dq, read_noise, rejection_thresh,
max_ratio > two_diff_rej_thresh))

log_str = 'From highest outlier, two-point found {} pixels with at least one CR from {} groups.'
log.info(log_str.format(len(row4cr), 'five'))
log.info(log_str.format(len(row4cr), 'five or more'))
log.info(log_str.format(len(row3cr), 'four'))
log.info(log_str.format(len(row2cr), 'three'))

Expand Down Expand Up @@ -255,6 +281,29 @@ def find_crs(dataa, group_dq, read_noise, rejection_thresh,
if (gdq[integ, group, row, col + 1] & dnu_flag) == 0:
gdq[integ, group, row, col + 1] =\
np.bitwise_or(gdq[integ, group, row, col + 1], jump_flag)

# flag n groups after jumps above the specified thresholds to account for
# the transient seen after ramp jumps
flag_e_threshold = [after_jump_flag_e1, after_jump_flag_e2]
flag_groups = [after_jump_flag_n1, after_jump_flag_n2]

cr_group, cr_row, cr_col = np.where(np.bitwise_and(gdq[integ], jump_flag))

for cthres, cgroup in zip(flag_e_threshold, flag_groups):
if cgroup > 0:
log.info(f"Flagging {cgroup} groups after detected jumps with e >= {np.mean(cthres)}.")

for j in range(len(cr_group)):
group = cr_group[j]
row = cr_row[j]
col = cr_col[j]
if e_jump[group - 1, row, col] >= cthres[row, col]:
for kk in range(group, min(group + cgroup + 1, ngroups)):
if (gdq[integ, kk, row, col] & sat_flag) == 0:
if (gdq[integ, kk, row, col] & dnu_flag) == 0:
gdq[integ, kk, row, col] =\
np.bitwise_or(gdq[integ, kk, row, col], jump_flag)

return gdq, row_below_gdq, row_above_gdq


Expand Down
116 changes: 116 additions & 0 deletions tests/test_twopoint_difference.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -835,6 +835,122 @@ def test_first_last_3group(setup_cube):
assert outgdq[0, 1, 0, 0] == 0


def test_10grps_1cr_afterjump(setup_cube):
ngroups = 10
data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, readnoise=10)
nframes = 1
data[0, 0, 100, 100] = 0
data[0, 1, 100, 100] = 10
data[0, 2, 100, 100] = 21
data[0, 3, 100, 100] = 33
data[0, 4, 100, 100] = 46
data[0, 5, 100, 100] = 60
data[0, 6, 100, 100] = 1160
data[0, 7, 100, 100] = 1175
data[0, 8, 100, 100] = 1190
data[0, 9, 100, 100] = 1209

after_jump_flag_e1 = np.full(data.shape[2:4], 1.0) * 0.0
out_gdq, row_below_gdq, row_above_gdq = find_crs(data, gdq, read_noise, rej_threshold,
rej_threshold, rej_threshold, nframes,
False, 200, 10, DQFLAGS,
after_jump_flag_e1=after_jump_flag_e1,
after_jump_flag_n1=10)
# all groups after CR should be flagged
for k in range(6, 10):
assert 4 == out_gdq[0, k, 100, 100], f"after jump flagging failed in group {k}"


def test_10grps_1cr_afterjump_2group(setup_cube):
ngroups = 10
data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, readnoise=10)
nframes = 1
data[0, 0, 100, 100] = 0
data[0, 1, 100, 100] = 10
data[0, 2, 100, 100] = 21
data[0, 3, 100, 100] = 33
data[0, 4, 100, 100] = 46
data[0, 5, 100, 100] = 60
data[0, 6, 100, 100] = 1160
data[0, 7, 100, 100] = 1175
data[0, 8, 100, 100] = 1190
data[0, 9, 100, 100] = 1209

after_jump_flag_e1 = np.full(data.shape[2:4], 1.0) * 0.0
out_gdq, row_below_gdq, row_above_gdq = find_crs(data, gdq, read_noise, rej_threshold,
rej_threshold, rej_threshold, nframes,
False, 200, 10, DQFLAGS,
after_jump_flag_e1=after_jump_flag_e1,
after_jump_flag_n1=2)

# 2 groups after CR should be flagged
for k in range(6, 9):
assert 4 == out_gdq[0, k, 100, 100], f"after jump flagging failed in group {k}"

# rest not flagged
for k in range(9, 10):
assert 0 == out_gdq[0, k, 100, 100], f"after jump flagging incorrect in group {k}"


def test_10grps_1cr_afterjump_toosmall(setup_cube):
ngroups = 10
data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, readnoise=10)
nframes = 1
data[0, 0, 100, 100] = 0
data[0, 1, 100, 100] = 10
data[0, 2, 100, 100] = 21
data[0, 3, 100, 100] = 33
data[0, 4, 100, 100] = 46
data[0, 5, 100, 100] = 60
data[0, 6, 100, 100] = 1160
data[0, 7, 100, 100] = 1175
data[0, 8, 100, 100] = 1190
data[0, 9, 100, 100] = 1209

after_jump_flag_e1 = np.full(data.shape[2:4], 1.0) * 10000.0
out_gdq, row_below_gdq, row_above_gdq = find_crs(data, gdq, read_noise, rej_threshold,
rej_threshold, rej_threshold, nframes,
False, 200, 10, DQFLAGS,
after_jump_flag_e1=after_jump_flag_e1,
after_jump_flag_n1=10)
# all groups after CR should be flagged
for k in range(7, 10):
assert 0 == out_gdq[0, k, 100, 100], f"after jump flagging incorrect in group {k}"


def test_10grps_1cr_afterjump_twothresholds(setup_cube):
ngroups = 10
data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, readnoise=10)
nframes = 1
data[0, 0, 100, 100] = 0
data[0, 1, 100, 100] = 10
data[0, 2, 100, 100] = 121
data[0, 3, 100, 100] = 133
data[0, 4, 100, 100] = 146
data[0, 5, 100, 100] = 160
data[0, 6, 100, 100] = 1160
data[0, 7, 100, 100] = 1175
data[0, 8, 100, 100] = 1190
data[0, 9, 100, 100] = 1209

after_jump_flag_e1 = np.full(data.shape[2:4], 1.0) * 500.
after_jump_flag_e2 = np.full(data.shape[2:4], 1.0) * 10.
out_gdq, row_below_gdq, row_above_gdq = find_crs(data, gdq, read_noise, rej_threshold,
rej_threshold, rej_threshold, nframes,
False, 200, 10, DQFLAGS,
after_jump_flag_e1=after_jump_flag_e1,
after_jump_flag_n1=10,
after_jump_flag_e2=after_jump_flag_e2,
after_jump_flag_n2=2)
# 2 groups after CR should be flagged
for k in range(2, 5):
assert 4 == out_gdq[0, k, 100, 100], f"after jump flagging incorrect in group {k}"

# all groups after CR should be flagged
for k in range(6, 10):
assert 4 == out_gdq[0, k, 100, 100], f"after jump flagging incorrect in group {k}"


def test_median_func():

"""
Expand Down

0 comments on commit 73a88c5

Please sign in to comment.