From d3225c67af42c5e951f135101503bd2563adeea7 Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Wed, 15 May 2024 14:05:51 +0100 Subject: [PATCH 01/11] Initial changes, splitting up flagging and fixing of bad pixels --- spaceKLIP/imagetools.py | 459 ++++++++++++++++++++++++++++++---------- 1 file changed, 347 insertions(+), 112 deletions(-) diff --git a/spaceKLIP/imagetools.py b/spaceKLIP/imagetools.py index b762b8fd..afd75dc8 100644 --- a/spaceKLIP/imagetools.py +++ b/spaceKLIP/imagetools.py @@ -26,6 +26,7 @@ from scipy.ndimage import gaussian_filter, median_filter, fourier_shift from scipy.ndimage import shift as spline_shift from scipy.optimize import leastsq, minimize +from scipy.interpolate import griddata from skimage.registration import phase_cross_correlation from spaceKLIP import utils as ut from spaceKLIP.psf import JWST_PSF @@ -766,40 +767,39 @@ def subtract_background(self, pass - - def fix_bad_pixels(self, - method='timemed+dqmed+medfilt', - bpclean_kwargs={}, - custom_kwargs={}, - timemed_kwargs={}, - dqmed_kwargs={}, - medfilt_kwargs={}, - types=['SCI', 'SCI_TA', 'SCI_BG', 'REF', 'REF_TA', 'REF_BG'], - subdir='bpcleaned', - restrict_to=None): + def find_bad_pixels(self, + method='dqarr', + dqarr_kwargs={}, + sigclip_kwargs={}, + custom_kwargs={}, + types=['SCI', 'SCI_TA', 'SCI_BG', 'REF', 'REF_TA', 'REF_BG'], + subdir='bpfound', + restrict_to=None): """ - Identify and fix bad pixels. + Identify bad pixels for cleaning Parameters ---------- method : str, optional - Sequence of bad pixel identification and cleaning methods to be run - on the data. Different methods must be joined by a '+' sign without + Sequence of bad pixel cleaning methods to be run on the data. + Different methods must be joined by a '+' sign without whitespace. Available methods are: + + - dqarr: uses DQ array to identify bad pixels + + - sigclip: use sigma clipping to identify additional bad pixels. - - bpclean: use sigma clipping to identify additional bad pixels. - - custom: use a custom bad pixel map. - - timemed: replace pixels which are only bad in some frames with - their median value from the good frames. - - dqmed: replace bad pixels with the median value of their - surrounding good pixels. - - medfilt: replace bad pixels with an image plane median filter. + - custom: use a custom bad pixel map - The default is 'timemed+dqmed+medfilt'. - bpclean_kwargs : dict, optional - Keyword arguments for the 'bpclean' method. Available keywords are: + The default is 'dqarr'. + dqarr_kwargs : dict, optional + Keyword arguments for the 'dqarr' identification method. Available keywords are: + + The default is {}. + sigclip_kwargs : dict, optional + Keyword arguments for the 'sigclip' identification methods. Available keywords are: - - sigclip : float, optional + - sigma: float, optional Sigma clipping threshold. The default is 5. - shift_x : list of int, optional Pixels in x-direction to which each pixel shall be compared to. @@ -814,22 +814,112 @@ def fix_bad_pixels(self, match the keys of the observations database and the dictionary content must be binary bad pixel maps (1 = bad, 0 = good) with the same shape as the corresponding data. The default is {}. + + The default is {}. + types : list of str, optional + List of data types for which bad pixels shall be identified. + The default is ['SCI', 'SCI_TA', 'SCI_BG', 'REF', 'REF_TA', 'REF_BG']. + subdir : str, optional + Name of the directory where the data products shall be saved. The + default is 'bpfound'. + + Returns + ------- + None + + """ + # Set output directory. + output_dir = os.path.join(self.database.output_dir, subdir) + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + # Loop through concatenations. + for i, key in enumerate(self.database.obs.keys()): + # if we limit to only processing some concatenations, + # check whether this concatenation matches the pattern + if (restrict_to is not None) and (restrict_to not in key): + continue + + log.info('--> Concatenation ' + key) + + # Loop through FITS files. + nfitsfiles = len(self.database.obs[key]) + for j in range(nfitsfiles): + + # Read FITS file and PSF mask. + fitsfile = self.database.obs[key]['FITSFILE'][j] + data, erro, pxdq, head_pri, head_sci, is2d, imshifts, maskoffs = ut.read_obs(fitsfile) + maskfile = self.database.obs[key]['MASKFILE'][j] + mask = ut.read_msk(maskfile) + + # Make copy of DQ array filled with zeros, i.e. all good pixels + pxdq_temp = np.zeros_like(pxdq) + + # Skip file types that are not in the list of types. + if self.database.obs[key]['TYPE'][j] in types: + # Call bad pixel identification routines. + method_split = method.split('+') + for k in range(len(method_split)): + head, tail = os.path.split(fitsfile) + if method_split[k] == 'dqarr': + log.info(' --> Method ' + method_split[k] + ': ' + tail) + # Flag any pixels marked as DO_NOT_USE that aren't NONSCIENCE + pxdq_temp = (np.isnan(data) | (pxdq_temp & 1 == 1)) \ + & np.logical_not(pxdq_temp & 512 == 512) + elif method_split[k] == 'sigclip': + log.info(' --> Method ' + method_split[k] + ': ' + tail) + self.find_bad_pixels_sigclip(data, erro, pxdq_temp, pxdq & 512 == 512, sigclip_kwargs) + elif method_split[k] == 'custom': + log.info(' --> Method ' + method_split[k] + ': ' + tail) + if self.database.obs[key]['TYPE'][j] not in ['SCI_TA', 'REF_TA']: + self.find_bad_pixels_custom(data, erro, pxdq_temp, key, custom_kwargs) + else: + log.info(' --> Method ' + method_split[k] + ': skipped because TA file') + else: + log.info(' --> Unknown method ' + method_split[k] + ': skipped') + + # The new DQ will just be the pxdq_temp we've been modifying + new_dq = pxdq_temp.astype(np.uint32) + + # Write FITS file and PSF mask. + fitsfile = ut.write_obs(fitsfile, output_dir, data, erro, new_dq, head_pri, head_sci, is2d, imshifts, maskoffs) + maskfile = ut.write_msk(maskfile, mask, fitsfile) + + # Update spaceKLIP database. + self.database.update_obs(key, j, fitsfile, maskfile) + + pass + + def fix_bad_pixels(self, + method='timemed+medfilt', + timemed_kwargs={}, + localmed_kwargs={}, + medfilt_kwargs={}, + interp2d_kwargs={}, + types=['SCI', 'SCI_TA', 'SCI_BG', 'REF', 'REF_TA', 'REF_BG'], + subdir='bpcleaned', + restrict_to=None): + """ + Clean bad pixels. + + Parameters + ---------- + method : str, optional + Sequence of bad pixel cleaning methods to be run on the data. + Different methods must be joined by a '+' sign without + whitespace. Available methods are: + + - timemed: replace pixels which are only bad in some frames with + their median value from the good frames. + + - medfilt: replace bad pixels with an image plane median filter. + + The default is 'timemed+medfilt'. timemed_kwargs : dict, optional Keyword arguments for the 'timemed' method. Available keywords are: - n/a - The default is {}. - dqmed_kwargs : dict, optional - Keyword arguments for the 'dqmed' method. Available keywords are: - - - shift_x : list of int, optional - Pixels in x-direction from which the median shall be computed. - The default is [-1, 0, 1]. - - shift_y : list of int, optional - Pixels in y-direction from which the median shall be computed. - The default is [-1, 0, 1]. - The default is {}. medfilt_kwargs : dict, optional Keyword arguments for the 'medfilt' method. Available keywords are: @@ -874,55 +964,65 @@ def fix_bad_pixels(self, data, erro, pxdq, head_pri, head_sci, is2d, imshifts, maskoffs = ut.read_obs(fitsfile) maskfile = self.database.obs[key]['MASKFILE'][j] mask = ut.read_msk(maskfile) + + fig = plt.figure() + ax = plt.gca() + ax.hist(data.flatten(), + bins=int(np.sqrt(len(data.flatten()))), + histtype='step', + label='Pre Cleaning') + + # Make copy of DQ array + pxdq_temp = pxdq.copy() # Skip file types that are not in the list of types. if self.database.obs[key]['TYPE'][j] in types: - - # Call bad pixel cleaning routines. - pxdq_temp = pxdq.copy() - # if self.database.obs[key]['TELESCOP'][j] == 'JWST' and self.database.obs[key]['INSTRUME'][j] == 'NIRCAM': - # pxdq_temp = (pxdq_temp != 0) & np.logical_not(pxdq_temp & 512 == 512) - # else: - pxdq_temp = (np.isnan(data) | (pxdq_temp & 1 == 1)) & np.logical_not(pxdq_temp & 512 == 512) method_split = method.split('+') + + spatial = ['localmed', 'medfilt', 'interp2d'] + # If localmed and medfilt in cleaning, can't run both + if len(set(method_split) & set(spatial)) > 1: + log.info(' --> The localmed/medfilt/interp2d methods are redundant,') + log.info(' only the first method you listed will be affect the data.') + + # Loop over methods for k in range(len(method_split)): head, tail = os.path.split(fitsfile) - if method_split[k] == 'bpclean': - log.info(' --> Method ' + method_split[k] + ': ' + tail) - self.find_bad_pixels_bpclean(data, erro, pxdq_temp, pxdq & 512 == 512, bpclean_kwargs) - elif method_split[k] == 'custom': - log.info(' --> Method ' + method_split[k] + ': ' + tail) - if self.database.obs[key]['TYPE'][j] not in ['SCI_TA', 'REF_TA']: - self.find_bad_pixels_custom(data, erro, pxdq_temp, key, custom_kwargs) - else: - log.info(' --> Method ' + method_split[k] + ': skipped because TA file') - elif method_split[k] == 'timemed': - log.info(' --> Method ' + method_split[k] + ': ' + tail) + log.info(' --> Method ' + method_split[k] + ': ' + tail) + if method_split[k] == 'timemed': self.fix_bad_pixels_timemed(data, erro, pxdq_temp, timemed_kwargs) - elif method_split[k] == 'dqmed': - log.info(' --> Method ' + method_split[k] + ': ' + tail) - self.fix_bad_pixels_dqmed(data, erro, pxdq_temp, dqmed_kwargs) + elif method_split[k] == 'localmed': + self.fix_bad_pixels_localmed(data, erro, pxdq_temp, localmed_kwargs) elif method_split[k] == 'medfilt': - log.info(' --> Method ' + method_split[k] + ': ' + tail) self.fix_bad_pixels_medfilt(data, erro, pxdq_temp, medfilt_kwargs) + elif method_split[k] == 'interp2d': + self.fix_bad_pixels_interp2d(data, erro, pxdq_temp, interp2d_kwargs) else: log.info(' --> Unknown method ' + method_split[k] + ': skipped') - # if self.database.obs[key]['TELESCOP'][j] == 'JWST' and self.database.obs[key]['INSTRUME'][j] == 'NIRCAM': - # pxdq[(pxdq != 0) & np.logical_not(pxdq & 512 == 512) & (pxdq_temp == 0)] = 0 - # else: - # pxdq[(pxdq & 1 == 1) & np.logical_not(pxdq & 512 == 512) & (pxdq_temp == 0)] = 0 # update the pixel DQ bit flags for the output files. # The pxdq variable here is effectively just the DO_NOT_USE flag, discarding other bits. # We want to make a new dq which retains the other bits as much as possible. # first, retain all the other bits (bits greater than 1), then add in the new/cleaned DO_NOT_USE bit - import jwst.datamodels do_not_use = jwst.datamodels.dqflags.pixel['DO_NOT_USE'] new_dq = np.bitwise_and(pxdq.copy(), np.invert(do_not_use)) # retain all other bits except the do_not_use bit new_dq = np.bitwise_or(new_dq, pxdq_temp) # add in the do_not_use bit from the cleaned version new_dq = new_dq.astype(np.uint32) # ensure correct output type for saving # (the bitwise steps otherwise return np.int64 which isn't FITS compatible) + # Finish figure for this file + ax.hist(data.flatten(), + bins=int(np.sqrt(len(data.flatten()))), + histtype='step', + label='Post Cleaning') + ax.legend() + #ax.set_xscale('log') + ax.set_yscale('log') + ax.tick_params(which='both', direction='in', top=True, right=True, labelsize=12) + + output_file = os.path.join(output_dir, tail.replace('.fits','_hist.png')) + plt.savefig(output_file) + # Write FITS file and PSF mask. fitsfile = ut.write_obs(fitsfile, output_dir, data, erro, new_dq, head_pri, head_sci, is2d, imshifts, maskoffs) maskfile = ut.write_msk(maskfile, mask, fitsfile) @@ -932,12 +1032,12 @@ def fix_bad_pixels(self, pass - def find_bad_pixels_bpclean(self, + def find_bad_pixels_sigclip(self, data, erro, pxdq, NON_SCIENCE, - bpclean_kwargs={}): + sigclip_kwargs={}): """ Use an iterative sigma clipping algorithm to identify additional bad pixels in the data. @@ -954,11 +1054,13 @@ def find_bad_pixels_bpclean(self, NON_SCIENCE : 3D-array Input binary non-science pixel maps (1 = bad, 0 = good). Will not be modified by the routine. - bpclean_kwargs : dict, optional - Keyword arguments for the 'bpclean' method. Available keywords are: + sigclip_kwargs : dict, optional + Keyword arguments for the 'sigclip' method. Available keywords are: - - sigclip : float, optional + - sigma : float, optional Sigma clipping threshold. The default is 5. + - neg_sigma : float, optional + Sigma clipping threshold for negative outliers. The default is 1. - shift_x : list of int, optional Pixels in x-direction to which each pixel shall be compared to. The default is [-1, 0, 1]. @@ -975,26 +1077,28 @@ def find_bad_pixels_bpclean(self, """ # Check input. - if 'sigclip' not in bpclean_kwargs.keys(): - bpclean_kwargs['sigclip'] = 5. - if 'shift_x' not in bpclean_kwargs.keys(): - bpclean_kwargs['shift_x'] = [-1, 0, 1] - if 'shift_y' not in bpclean_kwargs.keys(): - bpclean_kwargs['shift_y'] = [-1, 0, 1] - if 0 not in bpclean_kwargs['shift_x']: - bpclean_kwargs['shift_x'] += [0] - if 0 not in bpclean_kwargs['shift_y']: - bpclean_kwargs['shift_y'] += [0] + if 'sigma' not in sigclip_kwargs.keys(): + sigclip_kwargs['sigma'] = 5. + if 'neg_sigma' not in sigclip_kwargs.keys(): + sigclip_kwargs['neg_sigma'] = 1. + if 'shift_x' not in sigclip_kwargs.keys(): + sigclip_kwargs['shift_x'] = [-1, 0, 1] + if 'shift_y' not in sigclip_kwargs.keys(): + sigclip_kwargs['shift_y'] = [-1, 0, 1] + if 0 not in sigclip_kwargs['shift_x']: + sigclip_kwargs['shift_x'] += [0] + if 0 not in sigclip_kwargs['shift_y']: + sigclip_kwargs['shift_y'] += [0] # Pad data. - pad_left = np.abs(np.min(bpclean_kwargs['shift_x'])) - pad_right = np.abs(np.max(bpclean_kwargs['shift_x'])) + pad_left = np.abs(np.min(sigclip_kwargs['shift_x'])) + pad_right = np.abs(np.max(sigclip_kwargs['shift_x'])) if pad_right == 0: right = None else: right = -pad_right - pad_bottom = np.abs(np.min(bpclean_kwargs['shift_y'])) - pad_top = np.abs(np.max(bpclean_kwargs['shift_y'])) + pad_bottom = np.abs(np.min(sigclip_kwargs['shift_y'])) + pad_top = np.abs(np.max(sigclip_kwargs['shift_y'])) if pad_top == 0: top = None else: @@ -1016,32 +1120,56 @@ def find_bad_pixels_bpclean(self, bg_ind = data[i] < (bg_med + 10. * bg_std) # clip bright PSFs for final calculation bg_med = np.nanmedian(data_temp[i][bg_ind]) bg_std = robust.medabsdev(data_temp[i][bg_ind]) - + # Create initial mask of large negative values. - ww[i] = ww[i] | (data[i] < bg_med - bpclean_kwargs['sigclip'] * bg_std) + ww[i] = ww[i] | (data[i] < bg_med - sigclip_kwargs['neg_sigma'] * bg_std) + ww[i][NON_SCIENCE[i]] = 0 # Loop through max 10 iterations. for it in range(10): data_temp[i][ww[i]] = np.nan erro_temp[i][ww[i]] = np.nan - # Shift data. + # Shift data and calculate median and standard deviation of neighbours pad_data = np.pad(data_temp[i], pad_vals, mode='edge') pad_erro = np.pad(erro_temp[i], pad_vals, mode='edge') data_arr = [] erro_arr = [] - for ix in bpclean_kwargs['shift_x']: - for iy in bpclean_kwargs['shift_y']: + for ix in sigclip_kwargs['shift_x']: + for iy in sigclip_kwargs['shift_y']: if ix != 0 or iy != 0: data_arr += [np.roll(pad_data, (iy, ix), axis=(0, 1))] erro_arr += [np.roll(pad_erro, (iy, ix), axis=(0, 1))] data_arr = np.array(data_arr) - data_arr = data_arr[:, pad_bottom:top, pad_left:right] - data_med = np.nanmedian(data_arr, axis=0) + data_arr_trim = data_arr[:, pad_bottom:top, pad_left:right] + data_med = np.nanmedian(data_arr_trim, axis=0) diff = data[i] - data_med - data_std = np.nanstd(data_arr, axis=0) - # data_std = robust.medabsdev(data_arr, axis=0) - mask_new = diff > bpclean_kwargs['sigclip'] * data_std + data_std = np.nanstd(data_arr_trim, axis=0) + + # fig, ax = plt.subplots(1, 2) + # ax[0].imshow(diff) + # ax[1].imshow(data_std) + # plt.show() + + # # Do the same for the diff array we just made + # pad_diff = np.pad(diff, pad_vals, mode='edge') + # diff_arr = [] + # for ix in sigclip_kwargs['shift_x']: + # for iy in sigclip_kwargs['shift_y']: + # if ix != 0 or iy != 0: + # diff_arr += [np.roll(pad_diff, (iy, ix), axis=(0, 1))] + # diff_arr = np.array(diff_arr) + # diff_arr = diff_arr[:, pad_bottom:top, pad_left:right] + # diff_med = np.nanmedian(diff_arr, axis=0) + # doublediff = data[i] - data_med - diff_med + # diff_std = np.nanstd(diff_arr, axis=0) + + # Find values N standard deviations above the mean of neighbors + threshold = sigclip_kwargs['sigma'] * data_std + mask_new = diff > threshold + + # data_temp[i][mask_new] = np.nan + nmask_new = np.sum(mask_new & np.logical_not(ww[i])) # print('Iteration %.0f: %.0f bad pixels identified, %.0f are new' % (it + 1, np.sum(mask_new), nmask_new)) sys.stdout.write('\rFrame %.0f/%.0f, iteration %.0f' % (i + 1, ww.shape[0], it + 1)) @@ -1052,7 +1180,7 @@ def find_bad_pixels_bpclean(self, ww[i][NON_SCIENCE[i]] = 0 pxdq[i][ww[i]] = 1 print('') - log.info(' --> Method bpclean: identified %.0f additional bad pixel(s) -- %.2f%%' % (np.sum(pxdq) - np.sum(pxdq_orig), 100. * (np.sum(pxdq) - np.sum(pxdq_orig)) / np.prod(pxdq.shape))) + log.info(' --> Method sigclip: identified %.0f additional bad pixel(s) -- %.2f%%' % (np.sum(pxdq) - np.sum(pxdq_orig), 100. * (np.sum(pxdq) - np.sum(pxdq_orig)) / np.prod(pxdq.shape))) pass @@ -1140,11 +1268,11 @@ def fix_bad_pixels_timemed(self, pass - def fix_bad_pixels_dqmed(self, + def fix_bad_pixels_localmed(self, data, erro, pxdq, - dqmed_kwargs={}): + localmed_kwargs={}): """ Replace bad pixels with the median value of their surrounding good pixels. @@ -1158,8 +1286,8 @@ def fix_bad_pixels_dqmed(self, pxdq : 3D-array Input binary bad pixel maps (1 = bad, 0 = good). Will be updated by the routine to exclude the fixed bad pixels. - dqmed_kwargs : dict, optional - Keyword arguments for the 'dqmed' method. Available keywords are: + localmed_kwargs : dict, optional + Keyword arguments for the 'localmed' method. Available keywords are: - shift_x : list of int, optional Pixels in x-direction from which the median shall be computed. @@ -1177,24 +1305,24 @@ def fix_bad_pixels_dqmed(self, """ # Check input. - if 'shift_x' not in dqmed_kwargs.keys(): - dqmed_kwargs['shift_x'] = [-1, 0, 1] - if 'shift_y' not in dqmed_kwargs.keys(): - dqmed_kwargs['shift_y'] = [-1, 0, 1] - if 0 not in dqmed_kwargs['shift_x']: - dqmed_kwargs['shift_x'] += [0] - if 0 not in dqmed_kwargs['shift_y']: - dqmed_kwargs['shift_y'] += [0] + if 'shift_x' not in localmed_kwargs.keys(): + localmed_kwargs['shift_x'] = [-1, 0, 1] + if 'shift_y' not in localmed_kwargs.keys(): + localmed_kwargs['shift_y'] = [-1, 0, 1] + if 0 not in localmed_kwargs['shift_x']: + localmed_kwargs['shift_x'] += [0] + if 0 not in localmed_kwargs['shift_y']: + localmed_kwargs['shift_y'] += [0] # Pad data. - pad_left = np.abs(np.min(dqmed_kwargs['shift_x'])) - pad_right = np.abs(np.max(dqmed_kwargs['shift_x'])) + pad_left = np.abs(np.min(localmed_kwargs['shift_x'])) + pad_right = np.abs(np.max(localmed_kwargs['shift_x'])) if pad_right == 0: right = None else: right = -pad_right - pad_bottom = np.abs(np.min(dqmed_kwargs['shift_y'])) - pad_top = np.abs(np.max(dqmed_kwargs['shift_y'])) + pad_bottom = np.abs(np.min(localmed_kwargs['shift_y'])) + pad_top = np.abs(np.max(localmed_kwargs['shift_y'])) if pad_top == 0: top = None else: @@ -1212,8 +1340,8 @@ def fix_bad_pixels_dqmed(self, for i in range(ww.shape[0]): data_arr = [] erro_arr = [] - for ix in dqmed_kwargs['shift_x']: - for iy in dqmed_kwargs['shift_y']: + for ix in localmed_kwargs['shift_x']: + for iy in localmed_kwargs['shift_y']: if ix != 0 or iy != 0: data_arr += [np.roll(pad_data[i], (iy, ix), axis=(0, 1))] erro_arr += [np.roll(pad_erro[i], (iy, ix), axis=(0, 1))] @@ -1227,7 +1355,7 @@ def fix_bad_pixels_dqmed(self, erro_med = np.nanmedian(erro_arr, axis=0) erro[i][ww[i]] = erro_med[ww[i]] pxdq[i][ww[i]] = 0 - log.info(' --> Method dqmed: fixing %.0f bad pixel(s) -- %.2f%%' % (np.sum(ww), 100. * np.sum(ww) / np.prod(ww.shape))) + log.info(' --> Method localmed: fixing %.0f bad pixel(s) -- %.2f%%' % (np.sum(ww), 100. * np.sum(ww) / np.prod(ww.shape))) pass @@ -1279,6 +1407,113 @@ def fix_bad_pixels_medfilt(self, pxdq[i][ww[i]] = 0 pass + + def fix_bad_pixels_interp2d(self, + data, + erro, + pxdq, + interp2d_kwargs={}): + """ + Replace bad pixels with an interpolation of neighbouring pixels. + + Parameters + ---------- + data : 3D-array + Input images. + erro : 3D-array + Input image uncertainties. + pxdq : 3D-array + Input binary bad pixel maps (1 = bad, 0 = good). Will be updated by + the routine to exclude the fixed bad pixels. + interp2d_kwargs : dict, optional + Keyword arguments for the 'interp2d' method. Available keywords are: + + - size : int, optional + Kernel size of the median filter to be used. The default is 4. + + The default is {}. + + Returns + ------- + None. + + """ + + # Check input. + if 'size' not in interp2d_kwargs.keys(): + interp2d_kwargs['size'] = 5 + + # Fix bad pixels using interpolation of neighbors. + ww = (pxdq != 0) & np.logical_not(pxdq & 512 == 512) + log.info(' --> Method interp2d: fixing %.0f bad pixel(s) -- %.2f%%' % (np.sum(ww), 100. * np.sum(ww) / np.prod(ww.shape))) + + # NaN pixels to be replaced with interpolation + data_temp = data.copy() + data_temp[np.where(np.isnan(data_temp))] = 0 + data_temp[ww] = np.nan + erro_temp = erro.copy() + erro_temp[np.where(np.isnan(erro_temp))] = 0 + erro_temp[ww] = np.nan + + rows, cols = data_temp[0].shape + half_box = interp2d_kwargs['size'] // 2 + for i in range(ww.shape[0]): + for ri in range(rows): + for ci in range(cols): + if np.isnan(data_temp[i][ri, ci]): + # Calculate the indices of the NxN box centered around the NaN pixel + x_min = max(0, ci - half_box) + x_max = min(cols, ci + half_box + 1) + y_min = max(0, ri - half_box) + y_max = min(rows, ri + half_box + 1) + + # Extract a NxN box within the valid range + box = data_temp[i][y_min:y_max, x_min:x_max] + ebox = erro_temp[i][y_min:y_max, x_min:x_max] + + # Extract coordinates and values from the box + box_coords = np.array(np.where(~np.isnan(box))).T \ + + np.array([[x_min, y_min]]) + box_values = box[~np.isnan(box)] + + ebox_coords = np.array(np.where(~np.isnan(ebox))).T \ + + np.array([[x_min, y_min]]) + ebox_values = ebox[~np.isnan(ebox)] + + # Perform interpolation if there are valid values in the box + if len(box_values) > interp2d_kwargs['size'] \ + and len(ebox_values) > interp2d_kwargs['size']: + # Extract x and y coordinates of valid values, same coords for + # data and err + x_coords = box_coords[:, 0] + y_coords = box_coords[:, 1] + ex_coords = ebox_coords[:, 0] + ey_coords = ebox_coords[:, 1] + + # Perform interpolation of data + data_interp = griddata((x_coords, y_coords), + box_values, + (ci, ri), + method='linear', + fill_value=np.nanmedian(box_values)) + + # Replace data pixel with interpolated value + data[i][ri, ci] = data_interp + + # Perform interpolation of error + err_interp = griddata((ex_coords, ey_coords), + ebox_values, + (ci, ri), + method='linear', + fill_value=np.nanmedian(ebox_values)) + + # Replace error pixel + erro[i][ri, ci] = err_interp + + + pxdq[i][ww[i]] = 0 + + pass def replace_nans(self, cval=0., From 0579e4b7edb1199c8281db7a136efb53ba85b08c Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Thu, 16 May 2024 10:33:21 +0100 Subject: [PATCH 02/11] Playing with pixel cleaning --- spaceKLIP/imagetools.py | 94 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 85 insertions(+), 9 deletions(-) diff --git a/spaceKLIP/imagetools.py b/spaceKLIP/imagetools.py index 6ab173d3..e1f31b6d 100644 --- a/spaceKLIP/imagetools.py +++ b/spaceKLIP/imagetools.py @@ -1007,6 +1007,8 @@ def find_bad_pixels(self, self.find_bad_pixels_custom(data, erro, pxdq_temp, key, custom_kwargs) else: log.info(' --> Method ' + method_split[k] + ': skipped because TA file') + elif method_split[k] == 'timeints': + self.find_bad_pixels_timeints(data, erro, pxdq_temp, key, custom_kwargs) else: log.info(' --> Unknown method ' + method_split[k] + ': skipped') @@ -1276,12 +1278,8 @@ def find_bad_pixels_sigclip(self, data_arr_trim = data_arr[:, pad_bottom:top, pad_left:right] data_med = np.nanmedian(data_arr_trim, axis=0) diff = data[i] - data_med - data_std = np.nanstd(data_arr_trim, axis=0) - # fig, ax = plt.subplots(1, 2) - # ax[0].imshow(diff) - # ax[1].imshow(data_std) - # plt.show() + data_std = np.nanstd(data_arr_trim, axis=0) # # Do the same for the diff array we just made # pad_diff = np.pad(diff, pad_vals, mode='edge') @@ -1300,7 +1298,12 @@ def find_bad_pixels_sigclip(self, threshold = sigclip_kwargs['sigma'] * data_std mask_new = diff > threshold - # data_temp[i][mask_new] = np.nan + data_temp[i][mask_new] = np.nan + + # fig, ax = plt.subplots(1, 2) + # ax[0].imshow(data_temp[i]) + # ax[1].imshow(data_std) + # plt.show() nmask_new = np.sum(mask_new & np.logical_not(ww[i])) # print('Iteration %.0f: %.0f bad pixels identified, %.0f are new' % (it + 1, np.sum(mask_new), nmask_new)) @@ -1315,7 +1318,80 @@ def find_bad_pixels_sigclip(self, log.info(' --> Method sigclip: identified %.0f additional bad pixel(s) -- %.2f%%' % (np.sum(pxdq) - np.sum(pxdq_orig), 100. * (np.sum(pxdq) - np.sum(pxdq_orig)) / np.prod(pxdq.shape))) pass - + + def find_bad_pixels_timeints(self, + data, + erro, + pxdq, + NON_SCIENCE, + timeints_kwargs={}): + """ + Identify bad pixels from temporal variations across integrations. + + Parameters + ---------- + data : 3D-array + Input images. + erro : 3D-array + Input image uncertainties. + pxdq : 3D-array + Input binary bad pixel maps (1 = bad, 0 = good). Will be updated by + the routine to include the newly identified bad pixels. + NON_SCIENCE : 3D-array + Input binary non-science pixel maps (1 = bad, 0 = good). Will not + be modified by the routine. + timeints_kwargs : dict, optional + Keyword arguments for the 'timeints' method. Available keywords are: + + - sigma : float, optional + Sigma clipping threshold. The default is 5. + The default is {}. + + Returns + ------- + None. + + """ + + # Check input. + if 'sigma' not in timeints_kwargs.keys(): + timeints_kwargs['sigma'] = 10. + + pxdq_orig = pxdq.copy() + ww = pxdq != 0 + data_temp = data.copy() + data_temp[ww] = np.nan + + # Find bad pixels across the cube + + med_ints = np.nanmedian(data_temp, axis=0) + std_ints = np.nanstd(data_temp, axis=0) + + std2_ints = robust.medabsdev(data_temp, axis=0) + + diff = np.abs((data_temp - med_ints)) / std2_ints + + mask_new = diff > timeints_kwargs['sigma'] + + + data_temp[mask_new] = 9999 + plt.imshow(data_temp[1]) + plt.show() + # plt.hist(diff.flatten(), + # bins=int(np.sqrt(len(diff.flatten()))), + # histtype='step', + # label='Pre Cleaning') + # plt.yscale('log') + # plt.show() + + + ww = ww | mask_new + pxdq[ww] = 1 + print('') + log.info(' --> Method timeints: identified %.0f additional bad pixel(s) -- %.2f%%' % (np.sum(pxdq) - np.sum(pxdq_orig), 100. * (np.sum(pxdq) - np.sum(pxdq_orig)) / np.prod(pxdq.shape))) + + pass + def find_bad_pixels_custom(self, data, erro, @@ -1627,7 +1703,7 @@ def fix_bad_pixels_interp2d(self, box_values, (ci, ri), method='linear', - fill_value=np.nanmedian(box_values)) + fill_value=np.nan) # Replace data pixel with interpolated value data[i][ri, ci] = data_interp @@ -1637,7 +1713,7 @@ def fix_bad_pixels_interp2d(self, ebox_values, (ci, ri), method='linear', - fill_value=np.nanmedian(ebox_values)) + fill_value=np.nan) # Replace error pixel erro[i][ri, ci] = err_interp From dbe82b69a331bc6d2fe985aba6a741fad1c01861 Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Fri, 17 May 2024 10:57:14 +0100 Subject: [PATCH 03/11] End of Exeter visit: Seba gradient method incorporated --- spaceKLIP/coron1pipeline.py | 268 +++++++++++++++++++++++++++++++++--- spaceKLIP/imagetools.py | 96 ++++++++++++- spaceKLIP/pyklippipeline.py | 4 +- 3 files changed, 345 insertions(+), 23 deletions(-) diff --git a/spaceKLIP/coron1pipeline.py b/spaceKLIP/coron1pipeline.py index 45a3782a..76713f90 100644 --- a/spaceKLIP/coron1pipeline.py +++ b/spaceKLIP/coron1pipeline.py @@ -24,6 +24,7 @@ from .expjumpramp import ExperimentalJumpRampStep from webbpsf_ext import robust +from scipy.interpolate import interp1d from skimage.metrics import structural_similarity import logging @@ -894,16 +895,39 @@ def run_obs(database, # Need to do some preparation steps for group masking before running pipeline if 'mask_groups' in steps: - if ('groups_to_mask' not in steps['mask_groups']) and (groupmaskflag == 0): + if ('mask_array' not in steps['mask_groups']) and (groupmaskflag == 0): + # set a flag that we are running the group optimisation groupmaskflag = 1 - steps = prepare_group_masking(steps, database.obs[key], quiet) - + + # set method to basic if it wasn't set + if 'mask_method' not in steps['mask_groups']: + steps['mask_groups']['mask_method'] = 'basic' + # Also ensure certain files are skipped if necessary - if 'types' in steps['mask_groups']: - if database.obs[key]['TYPE'][j] not in steps['mask_groups']['types']: - steps['mask_groups']['skip'] = True - else: - steps['mask_groups']['skip'] = False + if not 'types' in steps['mask_groups']: + steps['mask_groups']['types'] = ['REF', 'REF_BG'] + + if database.obs[key]['TYPE'][j] not in steps['mask_groups']['types']: + steps['mask_groups']['skip'] = True + elif database.obs[key]['TYPE'][j] not in ['REF', 'REF_BG']: + log.info(' --> Group masking only works for reference images! Skipping...') + steps['mask_groups']['skip'] = True + else: + steps['mask_groups']['skip'] = False + + # Don't run function prep function if we don't need to + if not steps['mask_groups']['skip']: + if steps['mask_groups']['mask_method'] == 'basic': + steps = prepare_group_masking_basic(steps, + database.obs[key], + quiet) + elif steps['mask_groups']['mask_method'] == 'advanced': + fitstype = database.obs[key]['TYPE'][j] + steps = prepare_group_masking_advanced(steps, + database.obs[key], + fitspath, + fitstype, + quiet) # Get expected output file name outfile_name = tail.replace('uncal.fits', 'rateints.fits') @@ -920,18 +944,23 @@ def run_obs(database, if (j == jtervals[-1]) and (groupmaskflag == 1): '''This is the last file for this concatenation, and the groupmaskflag has been - set. This means we need to reset the groups_to_mask back to original state, - which was that it didn't exist. ''' + set. This means we need to reset the mask_array back to original state, + which was that it didn't exist, so that the routine is rerun. ''' groupmaskflag = 0 - del steps['mask_groups']['groups_to_mask'] + + if steps['mask_groups']['mask_method'] == 'basic': + del steps['mask_groups']['mask_array'] + elif steps['mask_groups']['mask_method'] == 'advanced': + del steps['mask_groups']['maxgrps_faint'] + del steps['mask_groups']['maxgrps_bright'] # Update spaceKLIP database. database.update_obs(key, j, fitsout_path) -def prepare_group_masking(steps, observations, quiet=False): +def prepare_group_masking_basic(steps, observations, quiet=False): - if 'groups_to_mask' not in steps['mask_groups']: + if 'mask_array' not in steps['mask_groups']: '''First time in the file loop, or groups_to_mask has been preset, run the optimisation and set groups to mask. ''' if not quiet: @@ -975,6 +1004,7 @@ def prepare_group_masking(steps, observations, quiet=False): elif observations['TYPE'][j] == 'REF': with fits.open(os.path.abspath(observations['FITSFILE'][j])) as hdul: ref_cube = hdul['SCI'].data.astype(float) + ref_shape = ref_cube.shape # Subtract a median so we focus on brightest pixels ref_cube -= np.nanmedian(ref_cube, axis=(2,3), keepdims=True) @@ -1019,12 +1049,214 @@ def prepare_group_masking(steps, observations, quiet=False): # Assemble array of groups to mask, starting one above the max group final_max_grp_to_use = int(np.nanmedian(max_grp_to_use)) groups_to_mask = np.arange(final_max_grp_to_use+1, ref_cubes[0].shape[1]) - + + # Make the mask array + mask_array = np.zeros(ref_shape, dtype=bool) + mask_array[:,groups_to_mask,:,:] = 1 + # Assign to steps so this stage doesn't get repeated. - steps['mask_groups']['groups_to_mask'] = groups_to_mask + steps['mask_groups']['mask_array'] = mask_array return steps +def prepare_group_masking_advanced(steps, observations, refpath, reftype, quiet=False): + + ''' + Advanced group masking method which computes the group mask on a pixel by pixel + and reference cube by reference cube basis + ''' + + if 'cropwidth' not in steps['mask_groups']: + steps['mask_groups']['cropwidth'] = 30 + if 'edgewidth' not in steps['mask_groups']: + steps['mask_groups']['edgewidth'] = 20 + if 'threshold' not in steps['mask_groups']: + steps['mask_groups']['threshold'] = 85 + + # Get crop width, part of image we care about + crop = steps['mask_groups']['cropwidth'] + edge = steps['mask_groups']['edgewidth'] + threshold = steps['mask_groups']['threshold'] + + + if ('maxgrps_faint' not in steps['mask_groups'] or + 'maxgrps_bright' not in steps['mask_groups'] or + 'cropmask' not in steps['mask_groups'] or + 'refbg_maxcounts' not in steps['mask_groups']): + '''First time in the file loop, or groups_to_mask has been preset, + run the optimisation and set groups to mask. ''' + + sci_frames = [] + sci_crpixs = [] + ref_cubes = [] + ref_crpixs = [] + + nfitsfiles = len(observations) + for j in range(nfitsfiles): + if observations['TYPE'][j] == 'SCI': + with fits.open(os.path.abspath(observations['FITSFILE'][j])) as hdul: + sci_frame = hdul['SCI'].data[:, -1, :, :].astype(float) + sci_crpix_x, sci_crpix_y = hdul["SCI"].header["CRPIX1"], hdul["SCI"].header["CRPIX2"] + sci_frames.append(sci_frame) + sci_crpixs.append([sci_crpix_x, sci_crpix_y]) + elif observations['TYPE'][j] == 'REF': + with fits.open(os.path.abspath(observations['FITSFILE'][j])) as hdul: + ref_cube = hdul['SCI'].data.astype(float) + ref_crpix_x, ref_crpix_y = hdul["SCI"].header["CRPIX1"], hdul["SCI"].header["CRPIX2"] + ref_cubes.append(ref_cube) + ref_crpixs.append([ref_crpix_x, ref_crpix_y]) + + # Crop and median science + sci_frames_cropped = [] + for i, scif in enumerate(sci_frames): + crpix_x, crpix_y = sci_crpixs[i] + xlo = int(crpix_x) - crop + xhi = int(crpix_x) + crop + ylo = int(crpix_y) - crop + yhi = int(crpix_y) + crop + sci_frames_cropped.append(scif[:, ylo:yhi, xlo:xhi]) + + sci_frames_modified = np.nanmedian(sci_frames_cropped, axis=1) #Median over integrations + + # Crop and median reference + ref_cubes_cropped = [] + for i, refc in enumerate(ref_cubes): + crpix_x, crpix_y = ref_crpixs[i] + xlo = int(crpix_x) - crop + xhi = int(crpix_x) + crop + ylo = int(crpix_y) - crop + yhi = int(crpix_y) + crop + ref_cubes_cropped.append(refc[:, :, ylo:yhi, xlo:xhi]) + + ref_cubes_modified = np.nanmedian(ref_cubes_cropped, axis=1) #Median over integrations + + # Median subtract the images + sci_frames_medsub = sci_frames_modified - np.nanmedian(sci_frames_modified, axis=(1,2), keepdims=True) + ref_cubes_medsub = ref_cubes_modified - np.nanmedian(ref_cubes_modified, axis=(2,3), keepdims=True) + + # Flatten array and find indices above percentile threshold value + sci_frames_flat = np.reshape(sci_frames_medsub, (sci_frames_medsub.shape[0], -1)) + per = np.percentile(sci_frames_flat, [threshold]) + above_threshold_indices = np.where(sci_frames_medsub > per) + + # Create an empty mask array + mask = np.zeros_like(sci_frames_medsub, dtype=bool) + + # Define function to expand indices and update mask + def expand_and_update_mask(indices, mask, xpad=2, ypad=2): + for z, x, y in zip(*indices): + mask[z, max(0, x-xpad):min(mask.shape[1], x+xpad), max(0, y-ypad):min(mask.shape[2], y+ypad)] = True + + # Expand indices and update mask for each 2D slice + for z_slice in range(sci_frames_medsub.shape[0]): + indices_slice = np.where(above_threshold_indices[0] == z_slice) + expand_and_update_mask((above_threshold_indices[0][indices_slice], above_threshold_indices[1][indices_slice], above_threshold_indices[2][indices_slice]), mask) + + # Okay now make some hollowed out cropped frames, to focus on fainter, but still bright, pixels + sci_frames_hollow = sci_frames_medsub.copy() + sci_frames_hollow[:, edge:-edge, edge:-edge] = np.nan + sci_frames_hollow[~mask] = np.nan + + ref_cubes_hollow = ref_cubes_medsub.copy() + ref_cubes_hollow[:, :, edge:-edge, edge:-edge] = np.nan + + # Create a 4D mask with zeros for reference + ref_cubes_shape = ref_cubes_hollow.shape + mask_4d = np.zeros(ref_cubes_shape, dtype=bool) + mask_ref = np.tile(mask[0:1], (ref_cubes_shape[0],1,1)) + + # FIX THIS TOMORROW! Average rolls? + + for i in range(mask_4d.shape[1]): + temp = mask_4d[:,i,:,:] + temp[mask_ref] = True + mask_4d[:,i,:,:] = temp + ref_cubes_hollow[~mask_4d] = np.nan + + # Now run the routine to figure out which groups to mask + best_faint_maxgrps = [] + best_bright_maxgrps = [] + ref_peak_pixels = [] + for i, scif in enumerate(sci_frames_hollow): + for j, refc in enumerate(ref_cubes_hollow): + # Need to save the peak pixel from each reference, as we'll use this for + # the REF_BG frame interpolations + if i == 0: + ref_peak_pixels.append(np.nanmax(ref_cubes_medsub[j][-1])) + this_faint_diffs= [] + this_bright_diffs = [] + for refg in refc: + faint_diff = np.abs(np.nanmedian(refg)-np.nanmedian(scif)) + this_faint_diffs.append(faint_diff) + for refg in ref_cubes_medsub[j]: + bright_diff = np.abs(np.nanmax(refg)-np.nanmax(sci_frames_medsub[i])) + this_bright_diffs.append(bright_diff) + + best_faint_maxgrp = np.argmin(this_faint_diffs) + best_faint_maxgrps.append(best_faint_maxgrp) + + best_bright_maxgrp = np.argmin(this_bright_diffs) + best_bright_maxgrps.append(best_bright_maxgrp) + + maxgrps_faint = int(np.nanmedian(best_faint_maxgrps)) + maxgrps_bright = int(np.nanmedian(best_bright_maxgrps)) + + steps['mask_groups']['maxgrps_faint'] = maxgrps_faint + steps['mask_groups']['maxgrps_bright'] = maxgrps_bright + steps['mask_groups']['cropmask'] = mask + steps['mask_groups']['refbg_maxcounts'] = np.nanmedian(ref_peak_pixels) + + # Now read in the specific reference file we're looking at. + with fits.open(refpath) as hdul: + refshape = hdul['SCI'].data.shape + ref_ints_slice = hdul['SCI'].data[:,-1,:,:].astype(float) + ref_slice = np.nanmedian(ref_ints_slice, axis=0) + ref_slice -= np.nanmedian(ref_slice) + ref_crpix_x, ref_crpix_y = hdul["SCI"].header["CRPIX1"], hdul["SCI"].header["CRPIX2"] + + # Get the peak pixel count in the last group, only look at the PSF core + xlo = int(ref_crpix_x) - crop + xhi = int(ref_crpix_x) + crop + ylo = int(ref_crpix_y) - crop + yhi = int(ref_crpix_y) + crop + ref_slice_cropped = ref_slice[ylo:yhi, xlo:xhi] + + if 'BG' in reftype: + maxcounts = steps['mask_groups']['refbg_maxcounts'] + else: + maxcounts = np.nanmax(ref_slice_cropped) + + # Get the median pixel count in our mask area from earlier + ref_slice_hollow = ref_slice_cropped.copy() + ref_slice_hollow[edge:-edge, edge:-edge] = np.nan + all_mincounts = [] + for saved_mask in steps['mask_groups']['cropmask']: + ref_slice_masked = ref_slice_hollow.copy() + ref_slice_masked[~saved_mask] = np.nan + all_mincounts.append(np.nanmedian(ref_slice_hollow)) + mincounts = np.nanmedian(all_mincounts) + + # Now make an interpolation connecting counts to the number of groups to be masked + maxgrps_interp = interp1d([maxcounts, mincounts], + [steps['mask_groups']['maxgrps_bright'],steps['mask_groups']['maxgrps_faint']], + kind='linear', + bounds_error=False, + fill_value=(steps['mask_groups']['maxgrps_bright'],steps['mask_groups']['maxgrps_faint'])) + + # Now use the interpolation to set the mask array, zero values will be included in the ramp fit + mask_array = np.zeros(refshape, dtype=bool) + for ri in range(mask_array.shape[2]): + for ci in range(mask_array.shape[3]): + if ref_slice[ri, ci] >= mincounts: + # Determine number of groups + this_grps = int(maxgrps_interp(ref_slice[ri, ci])) + groups_to_mask = np.arange(this_grps+1, refshape[1]) + mask_array[:,groups_to_mask,ri,ci] = 1 + + # Assign the mask array to the steps dictionary + steps['mask_groups']['mask_array'] = mask_array + + return steps class MaskGroupsStep(Step): """ @@ -1033,15 +1265,17 @@ class MaskGroupsStep(Step): class_alias = "maskgroups" spec = """ + mask_sigma_med = float(default=3) #Only mask pixels Nsigma above the median + mask_window = integer(default=2) #Also mask pixels within N pixels of a masked pixel """ def process(self, input): """Mask particular groups prior to ramp fitting""" with datamodels.open(input) as input_model: datamodel = input_model.copy() - + # Set particular groups to DO_NOT_USE - datamodel.groupdq[:,self.groups_to_mask,:,:] = 1 + datamodel.groupdq[self.mask_array] = 1 return datamodel diff --git a/spaceKLIP/imagetools.py b/spaceKLIP/imagetools.py index e1f31b6d..d4ea3414 100644 --- a/spaceKLIP/imagetools.py +++ b/spaceKLIP/imagetools.py @@ -904,6 +904,8 @@ def find_bad_pixels(self, dqarr_kwargs={}, sigclip_kwargs={}, custom_kwargs={}, + timeints_kwargs={}, + gradient_kwargs={}, types=['SCI', 'SCI_TA', 'SCI_BG', 'REF', 'REF_TA', 'REF_BG'], subdir='bpfound', restrict_to=None): @@ -1008,7 +1010,9 @@ def find_bad_pixels(self, else: log.info(' --> Method ' + method_split[k] + ': skipped because TA file') elif method_split[k] == 'timeints': - self.find_bad_pixels_timeints(data, erro, pxdq_temp, key, custom_kwargs) + self.find_bad_pixels_timeints(data, erro, pxdq_temp, key, timeints_kwargs) + elif method_split[k] == 'gradient': + self.find_bad_pixels_gradient(data, erro, pxdq_temp, key, gradient_kwargs) else: log.info(' --> Unknown method ' + method_split[k] + ': skipped') @@ -1374,9 +1378,9 @@ def find_bad_pixels_timeints(self, mask_new = diff > timeints_kwargs['sigma'] - data_temp[mask_new] = 9999 - plt.imshow(data_temp[1]) - plt.show() + # data_temp[mask_new] = 9999 + # plt.imshow(data_temp[1]) + # plt.show() # plt.hist(diff.flatten(), # bins=int(np.sqrt(len(diff.flatten()))), # histtype='step', @@ -1392,6 +1396,90 @@ def find_bad_pixels_timeints(self, pass + def find_bad_pixels_gradient(self, + data, + erro, + pxdq, + key, + gradient_kwargs={}): + + # Check input. + if 'sigma' not in gradient_kwargs.keys(): + gradient_kwargs['sigma'] = 0.5 + if 'threshold' not in gradient_kwargs.keys(): + gradient_kwargs['threshold'] = 0.05 + if 'negative' not in gradient_kwargs.keys(): + gradient_kwargs['negative'] = True + + sig = gradient_kwargs['sigma'] + threshold = gradient_kwargs['threshold'] + negative = gradient_kwargs['negative'] + + pxdq_orig = pxdq.copy() + ww = pxdq != 0 + data_temp = data.copy() + data_temp[ww] = np.nan + + # Loop over the images + for i in range(ww.shape[0]): + image = data_temp[i] + + ### remove nans + x=np.arange(0, image.shape[1]) + y=np.arange(0, image.shape[0]) + + xx, yy = np.meshgrid(x, y) + + # mask nans + image = np.ma.masked_invalid(image) + + xvalid=xx[~image.mask] + yvalid=yy[~image.mask] + + newimage=image[~image.mask] + + image_no_nans = griddata((xvalid, yvalid), + newimage.ravel(), + (xx, yy), + method='linear') + + ### get smooth image + smimage=gaussian_filter(image_no_nans, sigma=sig) + + ### get sharp image + shimage=image_no_nans-smimage + + ### get gradients + image_to_gradient=shimage/smimage + + gr=np.gradient((image_to_gradient)) + gr_dx=gr[1] + gr_dy=gr[0] + + ### pad gradient adding 1 extra pixel at beginning and end + gr_dxp=np.pad(gr_dx,(1,1)) + gr_dyp=np.pad(gr_dy,(1,1)) + + ### identify bad pixels + # positive + bad_pixels = (gr_dxp[1:-1,2:]<-threshold) & (gr_dxp[1:-1,:-2]>threshold) & (gr_dyp[2:,1:-1]<-threshold) & (gr_dyp[:-2,1:-1]>threshold) + # negative + if negative: + bad_pixels_n = (gr_dxp[1:-1,2:]>threshold) & (gr_dxp[1:-1,:-2]<-threshold) & (gr_dyp[2:,1:-1]>threshold) & (gr_dyp[:-2,1:-1] Method gradient: identified %.0f additional bad pixel(s) -- %.2f%%' % (np.sum(pxdq) - np.sum(pxdq_orig), 100. * (np.sum(pxdq) - np.sum(pxdq_orig)) / np.prod(pxdq.shape))) + pass + def find_bad_pixels_custom(self, data, erro, diff --git a/spaceKLIP/pyklippipeline.py b/spaceKLIP/pyklippipeline.py index ddf8ff11..bac3de6b 100644 --- a/spaceKLIP/pyklippipeline.py +++ b/spaceKLIP/pyklippipeline.py @@ -207,9 +207,9 @@ def run_obs(database, ww = [k for k in range(len(dataset._filenames)) if fitsfile in dataset._filenames[k]] hdul = fits.open(datapath) if dataset.allints.shape[1] == 1: - hdul[0].data = np.median(dataset.allints[:, :, ww, :, :], axis=(1, 2)) + hdul[0].data = np.nanmedian(dataset.allints[:, :, ww, :, :], axis=(1, 2)) else: - hdul[0].data = np.median(dataset.allints[:, :, ww, :, :], axis=2) + hdul[0].data = np.nanmedian(dataset.allints[:, :, ww, :, :], axis=2) hdul[0].header['NINTS'] = database.obs[key]['NINTS'][j] hdul[0].header['WCSAXES'] = head_sci['WCSAXES'] hdul[0].header['CRVAL1'] = head_sci['CRVAL1'] From b06cab0f6f9bd90042ec9066b662fd2a2bf1cabc Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Fri, 17 May 2024 11:02:11 +0100 Subject: [PATCH 04/11] Fixing file issue: --- spaceKLIP/coron1pipeline.py | 270 +++--------------------------------- 1 file changed, 17 insertions(+), 253 deletions(-) diff --git a/spaceKLIP/coron1pipeline.py b/spaceKLIP/coron1pipeline.py index 76713f90..eac4f5d6 100644 --- a/spaceKLIP/coron1pipeline.py +++ b/spaceKLIP/coron1pipeline.py @@ -1,5 +1,3 @@ -from __future__ import division - import matplotlib # ============================================================================= @@ -24,7 +22,6 @@ from .expjumpramp import ExperimentalJumpRampStep from webbpsf_ext import robust -from scipy.interpolate import interp1d from skimage.metrics import structural_similarity import logging @@ -895,39 +892,16 @@ def run_obs(database, # Need to do some preparation steps for group masking before running pipeline if 'mask_groups' in steps: - if ('mask_array' not in steps['mask_groups']) and (groupmaskflag == 0): - # set a flag that we are running the group optimisation + if ('groups_to_mask' not in steps['mask_groups']) and (groupmaskflag == 0): groupmaskflag = 1 - - # set method to basic if it wasn't set - if 'mask_method' not in steps['mask_groups']: - steps['mask_groups']['mask_method'] = 'basic' - + steps = prepare_group_masking(steps, database.obs[key], quiet) + # Also ensure certain files are skipped if necessary - if not 'types' in steps['mask_groups']: - steps['mask_groups']['types'] = ['REF', 'REF_BG'] - - if database.obs[key]['TYPE'][j] not in steps['mask_groups']['types']: - steps['mask_groups']['skip'] = True - elif database.obs[key]['TYPE'][j] not in ['REF', 'REF_BG']: - log.info(' --> Group masking only works for reference images! Skipping...') - steps['mask_groups']['skip'] = True - else: - steps['mask_groups']['skip'] = False - - # Don't run function prep function if we don't need to - if not steps['mask_groups']['skip']: - if steps['mask_groups']['mask_method'] == 'basic': - steps = prepare_group_masking_basic(steps, - database.obs[key], - quiet) - elif steps['mask_groups']['mask_method'] == 'advanced': - fitstype = database.obs[key]['TYPE'][j] - steps = prepare_group_masking_advanced(steps, - database.obs[key], - fitspath, - fitstype, - quiet) + if 'types' in steps['mask_groups']: + if database.obs[key]['TYPE'][j] not in steps['mask_groups']['types']: + steps['mask_groups']['skip'] = True + else: + steps['mask_groups']['skip'] = False # Get expected output file name outfile_name = tail.replace('uncal.fits', 'rateints.fits') @@ -944,23 +918,18 @@ def run_obs(database, if (j == jtervals[-1]) and (groupmaskflag == 1): '''This is the last file for this concatenation, and the groupmaskflag has been - set. This means we need to reset the mask_array back to original state, - which was that it didn't exist, so that the routine is rerun. ''' + set. This means we need to reset the groups_to_mask back to original state, + which was that it didn't exist. ''' groupmaskflag = 0 - - if steps['mask_groups']['mask_method'] == 'basic': - del steps['mask_groups']['mask_array'] - elif steps['mask_groups']['mask_method'] == 'advanced': - del steps['mask_groups']['maxgrps_faint'] - del steps['mask_groups']['maxgrps_bright'] + del steps['mask_groups']['groups_to_mask'] # Update spaceKLIP database. database.update_obs(key, j, fitsout_path) -def prepare_group_masking_basic(steps, observations, quiet=False): +def prepare_group_masking(steps, observations, quiet=False): - if 'mask_array' not in steps['mask_groups']: + if 'groups_to_mask' not in steps['mask_groups']: '''First time in the file loop, or groups_to_mask has been preset, run the optimisation and set groups to mask. ''' if not quiet: @@ -1004,7 +973,6 @@ def prepare_group_masking_basic(steps, observations, quiet=False): elif observations['TYPE'][j] == 'REF': with fits.open(os.path.abspath(observations['FITSFILE'][j])) as hdul: ref_cube = hdul['SCI'].data.astype(float) - ref_shape = ref_cube.shape # Subtract a median so we focus on brightest pixels ref_cube -= np.nanmedian(ref_cube, axis=(2,3), keepdims=True) @@ -1049,214 +1017,12 @@ def prepare_group_masking_basic(steps, observations, quiet=False): # Assemble array of groups to mask, starting one above the max group final_max_grp_to_use = int(np.nanmedian(max_grp_to_use)) groups_to_mask = np.arange(final_max_grp_to_use+1, ref_cubes[0].shape[1]) - - # Make the mask array - mask_array = np.zeros(ref_shape, dtype=bool) - mask_array[:,groups_to_mask,:,:] = 1 - + # Assign to steps so this stage doesn't get repeated. - steps['mask_groups']['mask_array'] = mask_array + steps['mask_groups']['groups_to_mask'] = groups_to_mask return steps -def prepare_group_masking_advanced(steps, observations, refpath, reftype, quiet=False): - - ''' - Advanced group masking method which computes the group mask on a pixel by pixel - and reference cube by reference cube basis - ''' - - if 'cropwidth' not in steps['mask_groups']: - steps['mask_groups']['cropwidth'] = 30 - if 'edgewidth' not in steps['mask_groups']: - steps['mask_groups']['edgewidth'] = 20 - if 'threshold' not in steps['mask_groups']: - steps['mask_groups']['threshold'] = 85 - - # Get crop width, part of image we care about - crop = steps['mask_groups']['cropwidth'] - edge = steps['mask_groups']['edgewidth'] - threshold = steps['mask_groups']['threshold'] - - - if ('maxgrps_faint' not in steps['mask_groups'] or - 'maxgrps_bright' not in steps['mask_groups'] or - 'cropmask' not in steps['mask_groups'] or - 'refbg_maxcounts' not in steps['mask_groups']): - '''First time in the file loop, or groups_to_mask has been preset, - run the optimisation and set groups to mask. ''' - - sci_frames = [] - sci_crpixs = [] - ref_cubes = [] - ref_crpixs = [] - - nfitsfiles = len(observations) - for j in range(nfitsfiles): - if observations['TYPE'][j] == 'SCI': - with fits.open(os.path.abspath(observations['FITSFILE'][j])) as hdul: - sci_frame = hdul['SCI'].data[:, -1, :, :].astype(float) - sci_crpix_x, sci_crpix_y = hdul["SCI"].header["CRPIX1"], hdul["SCI"].header["CRPIX2"] - sci_frames.append(sci_frame) - sci_crpixs.append([sci_crpix_x, sci_crpix_y]) - elif observations['TYPE'][j] == 'REF': - with fits.open(os.path.abspath(observations['FITSFILE'][j])) as hdul: - ref_cube = hdul['SCI'].data.astype(float) - ref_crpix_x, ref_crpix_y = hdul["SCI"].header["CRPIX1"], hdul["SCI"].header["CRPIX2"] - ref_cubes.append(ref_cube) - ref_crpixs.append([ref_crpix_x, ref_crpix_y]) - - # Crop and median science - sci_frames_cropped = [] - for i, scif in enumerate(sci_frames): - crpix_x, crpix_y = sci_crpixs[i] - xlo = int(crpix_x) - crop - xhi = int(crpix_x) + crop - ylo = int(crpix_y) - crop - yhi = int(crpix_y) + crop - sci_frames_cropped.append(scif[:, ylo:yhi, xlo:xhi]) - - sci_frames_modified = np.nanmedian(sci_frames_cropped, axis=1) #Median over integrations - - # Crop and median reference - ref_cubes_cropped = [] - for i, refc in enumerate(ref_cubes): - crpix_x, crpix_y = ref_crpixs[i] - xlo = int(crpix_x) - crop - xhi = int(crpix_x) + crop - ylo = int(crpix_y) - crop - yhi = int(crpix_y) + crop - ref_cubes_cropped.append(refc[:, :, ylo:yhi, xlo:xhi]) - - ref_cubes_modified = np.nanmedian(ref_cubes_cropped, axis=1) #Median over integrations - - # Median subtract the images - sci_frames_medsub = sci_frames_modified - np.nanmedian(sci_frames_modified, axis=(1,2), keepdims=True) - ref_cubes_medsub = ref_cubes_modified - np.nanmedian(ref_cubes_modified, axis=(2,3), keepdims=True) - - # Flatten array and find indices above percentile threshold value - sci_frames_flat = np.reshape(sci_frames_medsub, (sci_frames_medsub.shape[0], -1)) - per = np.percentile(sci_frames_flat, [threshold]) - above_threshold_indices = np.where(sci_frames_medsub > per) - - # Create an empty mask array - mask = np.zeros_like(sci_frames_medsub, dtype=bool) - - # Define function to expand indices and update mask - def expand_and_update_mask(indices, mask, xpad=2, ypad=2): - for z, x, y in zip(*indices): - mask[z, max(0, x-xpad):min(mask.shape[1], x+xpad), max(0, y-ypad):min(mask.shape[2], y+ypad)] = True - - # Expand indices and update mask for each 2D slice - for z_slice in range(sci_frames_medsub.shape[0]): - indices_slice = np.where(above_threshold_indices[0] == z_slice) - expand_and_update_mask((above_threshold_indices[0][indices_slice], above_threshold_indices[1][indices_slice], above_threshold_indices[2][indices_slice]), mask) - - # Okay now make some hollowed out cropped frames, to focus on fainter, but still bright, pixels - sci_frames_hollow = sci_frames_medsub.copy() - sci_frames_hollow[:, edge:-edge, edge:-edge] = np.nan - sci_frames_hollow[~mask] = np.nan - - ref_cubes_hollow = ref_cubes_medsub.copy() - ref_cubes_hollow[:, :, edge:-edge, edge:-edge] = np.nan - - # Create a 4D mask with zeros for reference - ref_cubes_shape = ref_cubes_hollow.shape - mask_4d = np.zeros(ref_cubes_shape, dtype=bool) - mask_ref = np.tile(mask[0:1], (ref_cubes_shape[0],1,1)) - - # FIX THIS TOMORROW! Average rolls? - - for i in range(mask_4d.shape[1]): - temp = mask_4d[:,i,:,:] - temp[mask_ref] = True - mask_4d[:,i,:,:] = temp - ref_cubes_hollow[~mask_4d] = np.nan - - # Now run the routine to figure out which groups to mask - best_faint_maxgrps = [] - best_bright_maxgrps = [] - ref_peak_pixels = [] - for i, scif in enumerate(sci_frames_hollow): - for j, refc in enumerate(ref_cubes_hollow): - # Need to save the peak pixel from each reference, as we'll use this for - # the REF_BG frame interpolations - if i == 0: - ref_peak_pixels.append(np.nanmax(ref_cubes_medsub[j][-1])) - this_faint_diffs= [] - this_bright_diffs = [] - for refg in refc: - faint_diff = np.abs(np.nanmedian(refg)-np.nanmedian(scif)) - this_faint_diffs.append(faint_diff) - for refg in ref_cubes_medsub[j]: - bright_diff = np.abs(np.nanmax(refg)-np.nanmax(sci_frames_medsub[i])) - this_bright_diffs.append(bright_diff) - - best_faint_maxgrp = np.argmin(this_faint_diffs) - best_faint_maxgrps.append(best_faint_maxgrp) - - best_bright_maxgrp = np.argmin(this_bright_diffs) - best_bright_maxgrps.append(best_bright_maxgrp) - - maxgrps_faint = int(np.nanmedian(best_faint_maxgrps)) - maxgrps_bright = int(np.nanmedian(best_bright_maxgrps)) - - steps['mask_groups']['maxgrps_faint'] = maxgrps_faint - steps['mask_groups']['maxgrps_bright'] = maxgrps_bright - steps['mask_groups']['cropmask'] = mask - steps['mask_groups']['refbg_maxcounts'] = np.nanmedian(ref_peak_pixels) - - # Now read in the specific reference file we're looking at. - with fits.open(refpath) as hdul: - refshape = hdul['SCI'].data.shape - ref_ints_slice = hdul['SCI'].data[:,-1,:,:].astype(float) - ref_slice = np.nanmedian(ref_ints_slice, axis=0) - ref_slice -= np.nanmedian(ref_slice) - ref_crpix_x, ref_crpix_y = hdul["SCI"].header["CRPIX1"], hdul["SCI"].header["CRPIX2"] - - # Get the peak pixel count in the last group, only look at the PSF core - xlo = int(ref_crpix_x) - crop - xhi = int(ref_crpix_x) + crop - ylo = int(ref_crpix_y) - crop - yhi = int(ref_crpix_y) + crop - ref_slice_cropped = ref_slice[ylo:yhi, xlo:xhi] - - if 'BG' in reftype: - maxcounts = steps['mask_groups']['refbg_maxcounts'] - else: - maxcounts = np.nanmax(ref_slice_cropped) - - # Get the median pixel count in our mask area from earlier - ref_slice_hollow = ref_slice_cropped.copy() - ref_slice_hollow[edge:-edge, edge:-edge] = np.nan - all_mincounts = [] - for saved_mask in steps['mask_groups']['cropmask']: - ref_slice_masked = ref_slice_hollow.copy() - ref_slice_masked[~saved_mask] = np.nan - all_mincounts.append(np.nanmedian(ref_slice_hollow)) - mincounts = np.nanmedian(all_mincounts) - - # Now make an interpolation connecting counts to the number of groups to be masked - maxgrps_interp = interp1d([maxcounts, mincounts], - [steps['mask_groups']['maxgrps_bright'],steps['mask_groups']['maxgrps_faint']], - kind='linear', - bounds_error=False, - fill_value=(steps['mask_groups']['maxgrps_bright'],steps['mask_groups']['maxgrps_faint'])) - - # Now use the interpolation to set the mask array, zero values will be included in the ramp fit - mask_array = np.zeros(refshape, dtype=bool) - for ri in range(mask_array.shape[2]): - for ci in range(mask_array.shape[3]): - if ref_slice[ri, ci] >= mincounts: - # Determine number of groups - this_grps = int(maxgrps_interp(ref_slice[ri, ci])) - groups_to_mask = np.arange(this_grps+1, refshape[1]) - mask_array[:,groups_to_mask,ri,ci] = 1 - - # Assign the mask array to the steps dictionary - steps['mask_groups']['mask_array'] = mask_array - - return steps class MaskGroupsStep(Step): """ @@ -1265,17 +1031,15 @@ class MaskGroupsStep(Step): class_alias = "maskgroups" spec = """ - mask_sigma_med = float(default=3) #Only mask pixels Nsigma above the median - mask_window = integer(default=2) #Also mask pixels within N pixels of a masked pixel """ def process(self, input): """Mask particular groups prior to ramp fitting""" with datamodels.open(input) as input_model: datamodel = input_model.copy() - + # Set particular groups to DO_NOT_USE - datamodel.groupdq[self.mask_array] = 1 + datamodel.groupdq[:,self.groups_to_mask,:,:] = 1 return datamodel From 4a00de3c6ac2a9c2495d49efae84e966009afd25 Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Thu, 27 Jun 2024 11:14:57 -0400 Subject: [PATCH 05/11] Revert "Merge branch 'miri_improvements' into badpixels" This reverts commit 9aba908f085594148d476ec259989c79f8cd3d24, reversing changes made to d3225c67af42c5e951f135101503bd2563adeea7. --- spaceKLIP/coron1pipeline.py | 190 +-- spaceKLIP/expjumpramp.py | 1301 ----------------- spaceKLIP/imagetools.py | 133 -- spaceKLIP/pyklippipeline.py | 1 - .../miri_bg_masks/godoy_mask_f1065c.fits | Bin 184320 -> 0 bytes .../miri_bg_masks/godoy_mask_f1140c.fits | Bin 184320 -> 0 bytes .../miri_bg_masks/godoy_mask_f1550c.fits | Bin 184320 -> 0 bytes spaceKLIP/utils.py | 61 +- 8 files changed, 15 insertions(+), 1671 deletions(-) delete mode 100644 spaceKLIP/expjumpramp.py delete mode 100644 spaceKLIP/resources/miri_bg_masks/godoy_mask_f1065c.fits delete mode 100644 spaceKLIP/resources/miri_bg_masks/godoy_mask_f1140c.fits delete mode 100644 spaceKLIP/resources/miri_bg_masks/godoy_mask_f1550c.fits diff --git a/spaceKLIP/coron1pipeline.py b/spaceKLIP/coron1pipeline.py index eac4f5d6..f13e171b 100644 --- a/spaceKLIP/coron1pipeline.py +++ b/spaceKLIP/coron1pipeline.py @@ -8,21 +8,17 @@ import pdb import sys -import astropy.io.fits as fits +import astropy.io.fits as pyfits import numpy as np from tqdm import trange from jwst.lib import reffile_utils -from jwst.stpipe import Step -from jwst import datamodels from jwst.datamodels import dqflags, RampModel, SaturationModel from jwst.pipeline import Detector1Pipeline, Image2Pipeline, Coron3Pipeline from .fnoise_clean import kTCSubtractStep, OneOverfStep -from .expjumpramp import ExperimentalJumpRampStep -from webbpsf_ext import robust -from skimage.metrics import structural_similarity +from webbpsf_ext import robust import logging log = logging.getLogger(__name__) @@ -72,10 +68,8 @@ def __init__(self, """ - self.step_defs['mask_groups'] = MaskGroupsStep self.step_defs['subtract_ktc'] = kTCSubtractStep self.step_defs['subtract_1overf'] = OneOverfStep - self.step_defs['experimental_jumpramp'] = ExperimentalJumpRampStep # Initialize Detector1Pipeline class. super(Coron1Pipeline_spaceKLIP, self).__init__(**kwargs) @@ -127,17 +121,17 @@ def process(self, input = self.run_step(self.group_scale, input) input = self.run_step(self.dq_init, input) input = self.run_step(self.saturation, input) - #input = self.run_step(self.ipc, input) Not run for MIRI + input = self.run_step(self.ipc, input) input = self.run_step(self.firstframe, input) input = self.run_step(self.lastframe, input) input = self.run_step(self.reset, input) input = self.run_step(self.linearity, input) input = self.run_step(self.rscd, input) input = self.run_step(self.dark_current, input) - input = self.run_step(self.refpix, input) - #input = self.run_step(self.charge_migration, input) Not run for MIRI + input = self.do_refpix(input) + input = self.run_step(self.charge_migration, input) input = self.run_step(self.jump, input) - input = self.run_step(self.mask_groups, input) + # TODO: Include same / similar subtract_1overf step??? else: input = self.run_step(self.group_scale, input) input = self.run_step(self.dq_init, input) @@ -151,7 +145,6 @@ def process(self, input = self.run_step(self.charge_migration, input) input = self.run_step(self.jump, input) input = self.run_step(self.subtract_ktc, input) - #1overf Only present in NIR data if 'groups' in self.stage_1overf: input = self.run_step(self.subtract_1overf, input) @@ -160,26 +153,19 @@ def process(self, self.save_model(input, suffix='ramp') # Run ramp fitting & gain scale correction. - if not self.experimental_jumpramp.use: - # Use the default ramp fitting procedure - res = self.run_step(self.ramp_fit, input, save_results=False) - rate, rateints = (res, None) if self.ramp_fit.skip else res - else: - # Use the experimental fitting procedure - res = self.run_step(self.experimental_jumpramp, input) - rate, rateints = res - + res = self.run_step(self.ramp_fit, input, save_results=False) + rate, rateints = (res, None) if self.ramp_fit.skip else res if self.rate_int_outliers and rateints is not None: # Flag additional outliers by comparing rateints and refit ramp input = self.apply_rateint_outliers(rateints, input) if input is None: - input, ints_model = rate, rateints + input, ints_model = (rate, rateints) else: res = self.run_step(self.ramp_fit, input, save_results=False) input, ints_model = (res, None) if self.ramp_fit.skip else res else: # input is the rate product, ints_model is the rateints product - input, ints_model = rate, rateints + input, ints_model = (rate, rateints) if input is None: self.ramp_fit.log.info('NoneType returned from ramp fitting. Gain scale correction skipped') @@ -673,14 +659,14 @@ def run_single_file(fitspath, output_dir, steps={}, verbose=False, **kwargs): # Skip dark current for subarray by default, but not full frame skip_dark = kwargs.get('skip_dark', None) if skip_dark is None: - hdr0 = fits.getheader(fitspath, ext=0) + hdr0 = pyfits.getheader(fitspath, ext=0) is_full_frame = 'FULL' in hdr0['SUBARRAY'] skip_dark = False if is_full_frame else True pipeline.dark_current.skip = skip_dark # Determine reference pixel correction parameters based on # instrument aperture name for NIRCam - hdr0 = fits.getheader(fitspath, 0) + hdr0 = pyfits.getheader(fitspath, 0) if hdr0['INSTRUME'] == 'NIRCAM': # Array of reference pixel borders [lower, upper, left, right] nb, nt, nl, nr = nrc_ref_info(hdr0['APERNAME'], orientation='sci') @@ -714,12 +700,6 @@ def run_single_file(fitspath, output_dir, steps={}, verbose=False, **kwargs): for key1 in steps.keys(): for key2 in steps[key1].keys(): setattr(getattr(pipeline, key1), key2, steps[key1][key2]) - - #Override jump & ramp if necessary - if pipeline.experimental_jumpramp.use: - log.info("Experimental jump/ramp fitting selected, regular jump and ramp will be skipped...") - pipeline.jump.skip = True - pipeline.ramp_fit.skip = True # Run Coron1Pipeline. Raise exception on error. # Ensure that pipeline is closed out. @@ -872,12 +852,11 @@ def run_obs(database, else: itervals = range(nkeys) - groupmaskflag = 0 # Set flag for group masking # Loop through concatenations. for i in itervals: key = keys[i] if not quiet: log.info('--> Concatenation ' + key) - + # Loop through FITS files. nfitsfiles = len(database.obs[key]) jtervals = trange(nfitsfiles, desc='FITS files', leave=False) if quiet else range(nfitsfiles) @@ -890,158 +869,17 @@ def run_obs(database, if not quiet: log.info(' --> Coron1Pipeline: skipping non-stage 0 file ' + tail) continue - # Need to do some preparation steps for group masking before running pipeline - if 'mask_groups' in steps: - if ('groups_to_mask' not in steps['mask_groups']) and (groupmaskflag == 0): - groupmaskflag = 1 - steps = prepare_group_masking(steps, database.obs[key], quiet) - - # Also ensure certain files are skipped if necessary - if 'types' in steps['mask_groups']: - if database.obs[key]['TYPE'][j] not in steps['mask_groups']['types']: - steps['mask_groups']['skip'] = True - else: - steps['mask_groups']['skip'] = False - # Get expected output file name outfile_name = tail.replace('uncal.fits', 'rateints.fits') fitsout_path = os.path.join(output_dir, outfile_name) # Skip if file already exists and overwrite is False. if os.path.isfile(fitsout_path) and not overwrite: - if not quiet: log.info(' --> Coron1Pipeline: skipping already processed file ' - + tail) + if not quiet: log.info(' --> Coron1Pipeline: skipping already processed file ' + tail) else: if not quiet: log.info(' --> Coron1Pipeline: processing ' + tail) _ = run_single_file(fitspath, output_dir, steps=steps, verbose=verbose, **kwargs) - - if (j == jtervals[-1]) and (groupmaskflag == 1): - '''This is the last file for this concatenation, and the groupmaskflag has been - set. This means we need to reset the groups_to_mask back to original state, - which was that it didn't exist. ''' - groupmaskflag = 0 - del steps['mask_groups']['groups_to_mask'] - # Update spaceKLIP database. database.update_obs(key, j, fitsout_path) - -def prepare_group_masking(steps, observations, quiet=False): - - if 'groups_to_mask' not in steps['mask_groups']: - '''First time in the file loop, or groups_to_mask has been preset, - run the optimisation and set groups to mask. ''' - if not quiet: - log.info(' --> Coron1Pipeline: Optimizing number of groups to mask in ramp,' - ' this make take a few minutes.') - - if 'cropwidth' not in steps['mask_groups']: - steps['mask_groups']['cropwidth'] = 20 - if 'edgewidth' not in steps['mask_groups']: - steps['mask_groups']['edgewidth'] = 10 - - # Get crop width, part of image we care about - crop = steps['mask_groups']['cropwidth'] - edge = steps['mask_groups']['edgewidth'] - - # Get our cropped science frames and reference cubes - sci_frames = [] - ref_cubes = [] - nfitsfiles = len(observations) - for j in range(nfitsfiles): - if observations['TYPE'][j] == 'SCI': - with fits.open(os.path.abspath(observations['FITSFILE'][j])) as hdul: - sci_frame = hdul['SCI'].data[:, -1, :, :].astype(float) - - # Subtract a median so we focus on brightest pixels - sci_frame -= np.nanmedian(sci_frame, axis=(1,2), keepdims=True) - - # Crop around CRPIX - crpix_x, crpix_y = hdul["SCI"].header["CRPIX1"], hdul["SCI"].header["CRPIX2"] - xlo = int(crpix_x) - crop - xhi = int(crpix_x) + crop - ylo = int(crpix_y) - crop - yhi = int(crpix_y) + crop - sci_frame = sci_frame[:, ylo:yhi, xlo:xhi] - - # Now going to set the core to 0, so we focus less on the highly variable - # PSF core - sci_frame[:, edge:-edge, edge:-edge] = np.nan - - sci_frames.append(sci_frame) - elif observations['TYPE'][j] == 'REF': - with fits.open(os.path.abspath(observations['FITSFILE'][j])) as hdul: - ref_cube = hdul['SCI'].data.astype(float) - - # Subtract a median so we focus on brightest pixels - ref_cube -= np.nanmedian(ref_cube, axis=(2,3), keepdims=True) - - # Crop around CRPIX - crpix_x, crpix_y = hdul["SCI"].header["CRPIX1"], hdul["SCI"].header["CRPIX2"] - xlo = int(crpix_x) - crop - xhi = int(crpix_x) + crop - ylo = int(crpix_y) - crop - yhi = int(crpix_y) + crop - ref_cube = ref_cube[:, :, ylo:yhi, xlo:xhi] - - # Now going to set the core to 0, so we focus less on the highly variable - # PSF core - ref_cube[:, :, edge:-edge, edge:-edge] = np.nan - - ref_cubes.append(ref_cube) - - # Want to check against every integration of every science dataset to find whichever - # matches the best, then use that for the scaling. - max_grp_to_use = [] - for sci_i, sci_frame in enumerate(sci_frames): - for int_i, sci_last_group in enumerate(sci_frame): - # Compare every reference group to this integration - best_diff = np.inf - for ref_cube in ref_cubes: - this_cube_diffs = [] - for ref_int in ref_cube: - this_int_diffs = [] - for ref_group in ref_int: - diff = np.abs(np.nansum(ref_group)-np.nansum(sci_last_group)) - this_int_diffs.append(diff) - this_cube_diffs.append(this_int_diffs) - - # Is the median of these diffs better that other reference cubes? - if np.nanmin(this_cube_diffs) < best_diff: - # If yes, this reference cube is a better match to the science - best_diff = np.nanmin(this_cube_diffs) - best_maxgrp = np.median(np.argmin(this_cube_diffs, axis=1)) - max_grp_to_use.append(best_maxgrp) - - # Assemble array of groups to mask, starting one above the max group - final_max_grp_to_use = int(np.nanmedian(max_grp_to_use)) - groups_to_mask = np.arange(final_max_grp_to_use+1, ref_cubes[0].shape[1]) - - # Assign to steps so this stage doesn't get repeated. - steps['mask_groups']['groups_to_mask'] = groups_to_mask - - return steps - - -class MaskGroupsStep(Step): - """ - Mask particular groups prior to ramp fitting - """ - class_alias = "maskgroups" - - spec = """ - """ - - def process(self, input): - """Mask particular groups prior to ramp fitting""" - with datamodels.open(input) as input_model: - datamodel = input_model.copy() - - # Set particular groups to DO_NOT_USE - datamodel.groupdq[:,self.groups_to_mask,:,:] = 1 - - return datamodel - - - diff --git a/spaceKLIP/expjumpramp.py b/spaceKLIP/expjumpramp.py deleted file mode 100644 index 41e57b7d..00000000 --- a/spaceKLIP/expjumpramp.py +++ /dev/null @@ -1,1301 +0,0 @@ -from __future__ import division - -from jwst.stpipe import Step -from jwst import datamodels -from jwst.datamodels import dqflags, RampModel -from stcal.ramp_fitting.utils import set_if_total_ramp -from jwst.ramp_fitting.ramp_fit_step import create_image_model, create_integration_model - -import astropy.io.fits as fits -import numpy as np -from scipy import special -import warnings - -from multiprocessing import Pool - -class ExperimentalJumpRampStep(Step): - """ - Experimental step based on Tim Brandt's jump detection - and ramp fitting algorithm. - """ - class_alias = "experimental_jumpramp" - - spec = """ - use = boolean(default=False) # Use regular pipeline by default - nproc = integer(default=4) # Number of processes to use - noise_scale = float(default=1) #Scale factor for noise estimate - nan_dnu = boolean(default=True) #Set do not use pixels to NaN? - """ - - def process(self, input): - """ Perform jump detection and ramp fitting - - Parameters - ---------- - input : JWST pipeline input - Appropriate input to JWST pipeline step, e.g. fits file name or preloaded - JWST datamodel. - - Returns - ------- - img_model : JWST datamodel - JWST datamodel equivalent to a "rate.fits" file. - int_model : JWST datamodel - JWST datamodel equivalent to a "rateints.fits" file. - - """ - with datamodels.RampModel(input) as input_model: - # Run the jump detection and ramp fitting - ramp_results = self.produce_ramp_images(input_model, self.noise_scale) - - # Update DQ arrays with jumps - input_model = self.update_groupdq(input_model, ramp_results['flagged_diffs']) - - # Convert DQ arrays to 3D - dq_ints = self.create_dq_ints(input_model) - - # Want to set do not use pixels to NaN - if self.nan_dnu: - ramp_results = self.set_dnu_to_nan(ramp_results, dq_ints) - - # Define integration info - integ_info = (ramp_results['sci'], - dq_ints, - ramp_results['var_poisson'], - ramp_results['var_rdnoise'], - ramp_results['err']) - - # Create integration model - int_model = create_integration_model(input_model, integ_info, input_model.int_times) - int_model.meta.bunit_data = 'DN/s' - int_model.meta.bunit_err = 'DN/s' - int_model.meta.cal_step.jump = 'COMPLETE' - int_model.meta.cal_step.ramp_fit = 'COMPLETE' - - # Define image info - im_sci, im_vpo, im_vrd, im_err = self.weighted_average(ramp_results['sci'], - ramp_results['var_poisson'], - ramp_results['var_rdnoise'], - ramp_results['err']) - image_info = (im_sci, - input_model.pixeldq, - im_vpo, - im_vrd, - im_err) - - # Create image model - img_model = create_image_model(input_model, image_info) - img_model.meta.bunit_data = 'DN/s' - img_model.meta.bunit_err = 'DN/s' - img_model.meta.cal_step.jump = 'COMPLETE' - img_model.meta.cal_step.ramp_fit = 'COMPLETE' - - return img_model, int_model - - def produce_ramp_images(self, datamodel, noise_scale): - """ - Function to assemble relevant data and run jump detection and ramp fitting. - - Parameters - ---------- - datamodel : JWST datamodel - JWST ramp datamodel - noise_scale : float - Scale factor to apply to the readnoise - - Returns - ------- - results : dict - Dictionary containing equivalents of the integration level 'SCI', 'ERR', - 'VAR_POISSON', and 'VAR_RDNOISE' arrays. - - """ - - # Assemble covariance and differenced frames. - C = self.create_covar(datamodel) - diffs = self.create_diffs(datamodel, C) - - # Grab reference files and create uncertainties - gain = self.get_subarray_ref(datamodel, 'gain') - rn = self.get_subarray_ref(datamodel, 'readnoise') - sig = noise_scale*gain*rn/2**0.5 - - # Need to multiply data by gain for this method, as it works in electrons - # Must undo this at the end!!! - diffs *= gain - - # Get some values - ncols = datamodel.meta.subarray.xsize - nrows = datamodel.meta.subarray.ysize - - if '-seg' in datamodel.meta.filename: - # This is a multi-segment dataset - int_s = datamodel.meta.exposure.integration_start - int_e = datamodel.meta.exposure.integration_end - nints = (int_e - int_s) + 1 - else: - # This is a single file dataset - nints = datamodel.meta.exposure.nints - - # Define arrays to save the results to - rate_ints = np.empty([nints, nrows, ncols]) - uncert_ints = np.empty([nints, nrows, ncols]) - var_po_ints = np.empty([nints, nrows, ncols]) - var_rd_ints = np.empty([nints, nrows, ncols]) - - # Define array to flag which diff frames to use - all_diffs2use = self.create_alldiffs2use(datamodel, diffs.shape) - - # Loop over each integration, run one ramp at a time - for int_i in range(nints): - print('--> Fitting Integration {}/{}...'.format(int_i+1, nints)) - - # Get the diffs for this integration - int_diffs = diffs[int_i] - - # Get array to save which differences we will use - int_diffs2use = all_diffs2use[int_i] - - # Now loop over each column using multiprocessing - colrange = range(ncols) - with Pool(processes=self.nproc) as pool: - results = pool.map(jumpramp_column_helper, [(i, - int_diffs, - C, - sig, - int_diffs2use) for i in colrange]) - pool.close() - pool.join() - - # Unpack results into the arrays - for i, res in enumerate(results): - res_rate, res_uncert, res_poisson, res_rdnoise, int_diffs2use = res - rate_ints[int_i, :, i] = res_rate - uncert_ints[int_i, :, i] = res_uncert - var_po_ints[int_i, :, i] = res_poisson - var_rd_ints[int_i, :, i] = res_rdnoise - all_diffs2use[int_i, :, :, i] = int_diffs2use - - # Explictly remove the gain scaling - rate_ints /= gain - uncert_ints /= gain - var_po_ints /= gain**2 - var_rd_ints /= gain**2 - - results = {'sci':rate_ints, - 'err':uncert_ints, - 'var_poisson':var_po_ints, - 'var_rdnoise':var_rd_ints, - 'flagged_diffs':all_diffs2use} - - return results - - def create_covar(self, datamodel): - """ Create the covariance matrix - - Parameters - ---------- - datamodel : JWST datamodel - JWST ramp datamodel - - Returns - ------- - C : Covar() object - Covariance matrix object - """ - - # Get readtimes array - dt = datamodel.meta.exposure.frame_time - ngroups = datamodel.meta.exposure.ngroups - readtimes = dt*(np.arange(ngroups)+1) - - # Produce covariance - C = Covar(readtimes) - - return C - - def create_diffs(self, datamodel, C): - """ Create the differenced frames - - Parameters - ---------- - datamodel : JWST datamodel - JWST ramp datamodel - C : Covar() object - Covariance matrix object - - Returns - ------- - diffs : ndarray - Scaled differenced frames of a given set of ramp images. - - """ - im = datamodel.data - - # Difference images - diffs = im[:,1:,:,:]-im[:,:-1,:,:] - - # Scale using covariance - diffs /= C.delta_t[np.newaxis, :, np.newaxis, np.newaxis] - - return diffs - - def get_subarray_ref(self, datamodel, filename, trim=True): - """ Get a reference file and extract on subarray - - Parameters - ---------- - datamodel : JWST datamodel - JWST ramp datamodel - filename : str - File name string for the given reference file, e.g. 'gain' or 'readnoise' - trim : bool - Boolean for trimming reference file to the subarray of the provided datamodel. - - Returns - ------- - subref : ndarray - Loaded reference file array. - - """ - - file = self.get_reference_file(datamodel, filename) - with fits.open(file) as hdul: - full = hdul['SCI'].data - - if trim == True: - # Trim full array to subarray for this data - xstrt = datamodel.meta.subarray.xstart-1 #SUBSTRT1 - ystrt = datamodel.meta.subarray.ystart-1 #SUBSTRT2 - xsize = datamodel.meta.subarray.xsize #SUBSIZE1 - ysize = datamodel.meta.subarray.ysize #SUBSIZE2 - subref = full[ystrt:ystrt+ysize, xstrt:xstrt+xsize] - else: - subref = full - - return subref - - def create_alldiffs2use(self, datamodel, shape): - """ - - Parameters - ---------- - datamodel : JWST datamodel - JWST ramp datamodel - shape : list - Shape of the corresponding diffs array - - Returns - ------- - all_diffs2use : ndarray - Array flagging which diffs should be used in the jump detection and ramp fitting - - """ - - # Make initial diffs2use - all_diffs2use = np.ones(shape, dtype=np.uint8) - - # Want to preflag any do not use pixels - dnu = dqflags.pixel['DO_NOT_USE'] - - grp_dq = datamodel.groupdq.copy() # Don't adjust actual groupdq - dnu_loc = np.bitwise_and(grp_dq, dnu) # Locations of do not use pixels per group - dnu_mask = np.where(dnu_loc > 0) - good_mask = np.where(dnu_loc == 0) - - grp_dq[dnu_mask] = 0 # These are bad - grp_dq[good_mask] = 1 # These are good - - # Add rather than subtract diffs - dq_diffs = grp_dq[:,1:,:,:]+grp_dq[:,:-1,:,:] - - # By adding, any diffs with a value <2 have at least one DNU pixel, so shouldn't be used. - diff_mask = np.where(dq_diffs < 2) - - # Set the diff mask diffs to zero so that aren't used. - all_diffs2use[diff_mask] = 0 - - return all_diffs2use - - def update_groupdq(self, datamodel, flagged_diffs): - """ Update group DQ based on detected jumps - - Parameters - ---------- - datamodel : JWST datamodel - JWST ramp datamodel - flagged_diffs : ndarray - Flagged array of size the same as differenced frames, with jumps marked. - - Returns - ------- - updated_datamodel : JWST datamodel - JWST ramp datamodel with updated groupdq - - """ - - # Create array to hold the jump information and define jump flag - dq_grp = datamodel.groupdq.copy() - jump = dqflags.pixel['JUMP_DET'] - - # Loop over each integration and create jump array - for int_i, integ in enumerate(flagged_diffs): - for i, diff in enumerate(integ): - # Locate the jumps in the differenced frames, Tim's code sets to 0. - jump_check = np.where(diff == 0) - - # If jump is in differenced frame N, then flag as a jump in real frames N and N+1 - this_dq_n = dq_grp[int_i, i] - this_dq_np1 = dq_grp[int_i, i+1] - - this_dq_n[jump_check] = np.bitwise_or(this_dq_n[jump_check], jump) - this_dq_np1[jump_check] = np.bitwise_or(this_dq_np1[jump_check], jump) - - dq_grp[int_i, i] = this_dq_n - dq_grp[int_i, i+1] = this_dq_np1 - - # Apply jump array to groupdq - datamodel.groupdq = dq_grp - - return datamodel - - def create_dq_ints(self, datamodel): - """ Turn the 2D and 4D pixel DQ to 3D DQ - - Parameters - ---------- - datamodel : JWST datamodel - JWST ramp datamodel - - Returns - ------- - dqints : ndarray - 3D data quality array, corresponding to each integration - - """ - - # Get 2D/4D arrays - dq_pix_base = datamodel.pixeldq - dq_grp = datamodel.groupdq - - # Define some DQ flags - sat = dqflags.pixel['SATURATED'] - jump = dqflags.pixel['JUMP_DET'] - dnu = dqflags.pixel['DO_NOT_USE'] - - # Define array to save values to - ncols = datamodel.meta.subarray.xsize - nrows = datamodel.meta.subarray.ysize - if '-seg' in datamodel.meta.filename: - # This is a multi-segment dataset - int_s = datamodel.meta.exposure.integration_start - int_e = datamodel.meta.exposure.integration_end - nints = (int_e - int_s) + 1 - else: - # This is a single file dataset - nints = datamodel.meta.exposure.nints - dq_ints = np.empty([nints, nrows, ncols], dtype=np.uint32) - - # Loop over integrations - for int_i in range(nints): - # Make copy of Pixel DQ - dq_pix = dq_pix_base.copy() - - # This integration DQ - dq_grp_int = dq_grp[int_i] - - # Check Saturated and Do Not Use - set_if_total_ramp(dq_pix, dq_grp_int, sat | dnu, dnu) - - # Assume complete saturation if group 0 saturated - dq_grp0_sat = np.bitwise_and(dq_grp_int[0], sat) - dq_pix[dq_grp0_sat != 0] = np.bitwise_or(dq_pix[dq_grp0_sat != 0], sat | dnu) - - # If jump occurs mark the appropriate flag. - jump_loc = np.bitwise_and(dq_grp_int, jump) - jump_check = np.where(jump_loc.sum(axis=0) > 0) - dq_pix[jump_check] = np.bitwise_or(dq_pix[jump_check], jump) - - # Save to main array - dq_ints[int_i] = dq_pix - - return dq_ints - - def set_dnu_to_nan(self, ramp_results, dq_ints): - ''' Set any do not use pixels to NaNs - - Parameters - ---------- - results : dict - Dictionary containing equivalents of the integration level 'SCI', 'ERR', - 'VAR_POISSON', and 'VAR_RDNOISE' arrays. - dq_ints : ndarray - 3D data quality array - - Returns - ------- - results : dict - Dictionary containing equivalents of the integration level 'SCI', 'ERR', - 'VAR_POISSON', and 'VAR_RDNOISE' arrays. - - ''' - - # Define do not use pixel check - dnu = dqflags.pixel['DO_NOT_USE'] - dnu_loc = np.bitwise_and(dq_ints, dnu) - dnu_check = np.where(dnu_loc > 0) - - # Assign NaNs to arrays - ramp_results['sci'][dnu_check] = np.nan - ramp_results['var_poisson'][dnu_check] = np.nan - ramp_results['var_rdnoise'][dnu_check] = np.nan - ramp_results['err'][dnu_check] = np.nan - - return ramp_results - - def weighted_average(self, sci, var_poisson, var_rdnoise, err): - ''' Compute a weighted average of a 3D multi-integration dataset - - Parameters - ---------- - sci : ndarray - Array of science data - var_poisson : ndarray - Array of poission variance - var_rdnoise : ndarray - Array of readnoise variance - err : ndarray - Array of uncertainties - - Returns - ------- - w_av_sci : ndarray - Weighted average array of science data - w_av_vpo : ndarray - Weighted average array of poisson variance - w_av_vrd : ndarray - Weighted average array of readnoise variance - w_av_err : ndarray - Weighted average array of uncertainties - - ''' - # Remove NaNs - sci = np.nan_to_num(sci) - var_poisson = np.nan_to_num(var_poisson) - var_rdnoise = np.nan_to_num(var_rdnoise) - err = np.nan_to_num(err) - - # Get normalised weights using variances and normalize - variance = var_poisson + var_rdnoise - weights = 1 / variance - weights /= np.sum(weights, axis=0) - - # Produce averaged results - w_av_sci = np.average(sci, weights=weights, axis=0) - w_av_err = np.sqrt(np.average(err**2, weights=weights, axis=0)) - w_av_vpo = np.average(var_poisson, weights=weights, axis=0) - w_av_vrd = np.average(var_rdnoise, weights=weights, axis=0) - - return w_av_sci, w_av_vpo, w_av_vrd, w_av_err - - -def jumpramp_column_helper(args): - """ Function to unpack parameters in multiprocessing - - Parameters - ---------- - args : list - List of arguments to be passed to the jumpramp calculation - - Returns - ------- - result.countrate : ndarray - Array of countrates from ramp fitting. - result.uncert : ndarray - Array of uncertainties from ramp fitting. - result.var_poisson : ndarray - Array of poisson variance from ramp fitting. - result.var_rdnoise : ndarray - Array of readnoise from ramp fitting. - int_diffs2use : ndarray - The integration level differenced frames to be used during ramp fitting. - - """ - - i, int_diffs, C, sig, int_diffs2use = args - result, int_diffs2use = jumpramp_column(i, int_diffs, C, sig, int_diffs2use) - return result.countrate, result.uncert, result.var_poisson, result.var_rdnoise, int_diffs2use - -def jumpramp_column(i, int_diffs, C, sig, int_diffs2use): - """ Function to run of jump/ramp on specific column of pixels, assumes an input 4D array of - data at the group level. - - Parameters - ---------- - i : int - Index of the integration to run ramp fitting over. - int_diffs : ndarray - Array of differenced frames - C : Covar() object - Covariance matrix object - sig : ndarray - Uncertainty on the data - int_diffs2use : ndarray - The integration level differenced frames to be used during ramp fitting. - - Returns - ------- - result : Ramp_Result() object - Output object from ramp fitting algorithm - - int_diffs2use : ndarray - Updated integration level differenced frames to be used during ramp fitting. - """ - - # Run jump detection - int_diffs2use, countrates = mask_jumps(int_diffs[:, :, i], - C, - sig[:, i], - diffs2use=int_diffs2use[:, :, i]) - - # Run ramp fitting - result = fit_ramps(int_diffs[:, :, i], - C, - sig[:, i], - diffs2use=int_diffs2use, - countrateguess=countrates*(countrates > 0)) - - return result, int_diffs2use - - -class Covar: - - """ - - class Covar holding read and photon noise components of alpha and - beta and the time intervals between the resultant midpoints - - """ - - def __init__(self, read_times, pedestal=False): - - """ - - Compute alpha and beta, the diagonal and off-diagonal elements of - the covariance matrix of the resultant differences, and the time - intervals between the resultant midpoints. - - Arguments: - 1. readtimes [list of values or lists for the times of reads. If - a list of lists, times for reads that are averaged - together to produce a resultant.] - - Optional arguments: - 2. pedestal [boolean: does the covariance matrix include the terms - for the first resultant? This is needed if fitting - for the pedestal (i.e. the reset value). Default - False. ] - - """ - - mean_t = [] # mean time of the resultant as defined in the paper - tau = [] # variance-weighted mean time of the resultant - N = [] # Number of reads per resultant - - for times in read_times: - mean_t += [np.mean(times)] - - if hasattr(times, "__len__"): - N += [len(times)] - k = np.arange(1, N[-1] + 1) - tau += [1/N[-1]**2*np.sum((2*N[-1] + 1 - 2*k)*np.array(times))] - else: - tau += [times] - N += [1] - - mean_t = np.array(mean_t) - tau = np.array(tau) - N = np.array(N) - - delta_t = mean_t[1:] - mean_t[:-1] - - self.pedestal = pedestal - self.delta_t = delta_t - self.mean_t = mean_t - self.tau = tau - self.Nreads = N - - self.alpha_readnoise = (1/N[:-1] + 1/N[1:])/delta_t**2 - self.beta_readnoise = -1/N[1:-1]/(delta_t[1:]*delta_t[:-1]) - - self.alpha_phnoise = (tau[:-1] + tau[1:] - 2*mean_t[:-1])/delta_t**2 - self.beta_phnoise = (mean_t[1:-1] - tau[1:-1])/(delta_t[1:]*delta_t[:-1]) - - # If we want the reset value we need to include the first - # resultant. These are the components of the variance and - # covariance for the first resultant. - - if pedestal: - self.alpha_readnoise = np.array([1/N[0]/mean_t[0]**2] - + list(self.alpha_readnoise)) - self.beta_readnoise = np.array([-1/N[0]/mean_t[0]/delta_t[0]] - + list(self.beta_readnoise)) - self.alpha_phnoise = np.array([tau[0]/mean_t[0]**2] - + list(self.alpha_phnoise)) - self.beta_phnoise = np.array([(mean_t[0] - tau[0])/mean_t[0]/delta_t[0]] - + list(self.beta_phnoise)) - - def calc_bias(self, countrates, sig, cvec, da=1e-7): - - """ - Calculate the bias in the best-fit count rate from estimating the - covariance matrix. This calculation is derived in the paper. - - Arguments: - 1. countrates [array of count rates at which the bias is desired] - 2. sig [float, single read noise] - 3. cvec [weight vector on resultant differences for initial - estimation of count rate for the covariance matrix. - Will be renormalized inside this function.] - Optional argument: - 4. da [float, fraction of the count rate plus sig**2 to use for finite - difference estimate of the derivative. Default 1e-7.] - - Returns: - 1. bias [array, bias of the best-fit count rate from using cvec - plus the observed resultants to estimate the covariance - matrix] - - """ - - if self.pedestal: - raise ValueError("Cannot compute bias with a Covar class that includes a pedestal fit.") - - alpha = countrates[np.newaxis, :]*self.alpha_phnoise[:, np.newaxis] - alpha += sig**2*self.alpha_readnoise[:, np.newaxis] - beta = countrates[np.newaxis, :]*self.beta_phnoise[:, np.newaxis] - beta += sig**2*self.beta_readnoise[:, np.newaxis] - - # we only want the weights; it doesn't matter what the count rates are. - n = alpha.shape[0] - z = np.zeros((len(cvec), len(countrates))) - result_low_a = fit_ramps(z, self, sig, countrateguess=countrates) - - # try to avoid problems with roundoff error - da_incr = da*(countrates[np.newaxis, :] + sig**2) - - dalpha = da_incr*self.alpha_phnoise[:, np.newaxis] - dbeta = da_incr*self.beta_phnoise[:, np.newaxis] - result_high_a = fit_ramps(z, self, sig, countrateguess=countrates+da_incr) - # finite difference approximation to dw/da - - dw_da = (result_high_a.weights - result_low_a.weights)/da_incr - - bias = np.zeros(len(countrates)) - c = cvec/np.sum(cvec) - - for i in range(len(countrates)): - - C = np.zeros((n, n)) - for j in range(n): - C[j, j] = alpha[j, i] - for j in range(n - 1): - C[j + 1, j] = C[j, j + 1] = beta[j, i] - - bias[i] = np.linalg.multi_dot([c[np.newaxis, :], C, dw_da[:, i]]) - - sig_a = np.sqrt(np.linalg.multi_dot([c[np.newaxis, :], C, c[:, np.newaxis]])) - bias[i] *= 0.5*(1 + special.erf(countrates[i]/sig_a/2**0.5)) - - return bias - - -class Ramp_Result: - - def __init__(self): - - self.countrate = None - self.chisq = None - self.uncert = None - self.var_poisson = None - self.var_rdnoise = None - self.weights = None - self.pedestal = None - self.uncert_pedestal = None - self.covar_countrate_pedestal = None - - self.countrate_twoomit = None - self.chisq_twoomit = None - self.uncert_twoomit = None - - self.countrate_oneomit = None - self.jumpval_oneomit = None - self.jumpsig_oneomit = None - self.chisq_oneomit = None - self.uncert_oneomit = None - - def fill_masked_reads(self, diffs2use): - - """ - - Replace countrates, uncertainties, and chi squared values that - are NaN because resultant differences were doubly omitted. - For these cases, revert to the corresponding values in with - fewer omitted resultant differences to get the correct values - without double-coundint omissions. - - Arguments: - 1. diffs2use [a 2D array matching self.countrate_oneomit in - shape with zero for resultant differences that - were masked and one for differences that were - not masked] - - This function replaces the relevant entries of - self.countrate_twoomit, self.chisq_twoomit, - self.uncert_twoomit, self.countrate_oneomit, and - self.chisq_oneomit in place. It does not return a value. - - """ - - # replace entries that would be nan (from trying to - # doubly exclude read differences) with the global fits. - - omit = diffs2use == 0 - ones = np.ones(diffs2use.shape) - - self.countrate_oneomit[omit] = (self.countrate*ones)[omit] - self.chisq_oneomit[omit] = (self.chisq*ones)[omit] - self.uncert_oneomit[omit] = (self.uncert*ones)[omit] - - omit = diffs2use[1:] == 0 - - self.countrate_twoomit[omit] = (self.countrate_oneomit[:-1])[omit] - self.chisq_twoomit[omit] = (self.chisq_oneomit[:-1])[omit] - self.uncert_twoomit[omit] = (self.uncert_oneomit[:-1])[omit] - - omit = diffs2use[:-1] == 0 - - self.countrate_twoomit[omit] = (self.countrate_oneomit[1:])[omit] - self.chisq_twoomit[omit] = (self.chisq_oneomit[1:])[omit] - self.uncert_twoomit[omit] = (self.uncert_oneomit[1:])[omit] - - -def fit_ramps(diffs, Cov, sig, countrateguess=None, diffs2use=None, - detect_jumps=False, resetval=0, resetsig=np.inf, rescale=True, - dn_scale=10): - - """ - Function fit_ramps. Fits ramps to read differences using the - covariance matrix for the read differences as given by the diagonal - elements and the off-diagonal elements. - - Arguments: - 1. diffs [resultant differences, shape (ndiffs, npix)] - 2. Cov [class Covar, holds the covariance matrix information] - 3. sig [read noise, 1D array, shape (npix)] - - Optional Arguments: - 4. countrateguess [array of shape (npix): count rates to be used - to estimate the covariance matrix. Default None, - in which case the average difference will be used, - replacing negative means with zeros.] - 5. diffs2use [shape (ndiffs, npix), boolean mask of whether to use - each resultant difference for each pixel. Default - None] - 6. detect_jumps [run the jump detection machinery leaving out - single and consecutive resultant differences. - Default False] - 7. resetval [float or array of shape (npix): priors on the reset - values. Default 0. Irrelevant unless - Cov.pedestal is True.] - 8. resetsig [float or array of shape (npix): uncertainties on the - reset values. Default np.inf, i.e., reset values - have flat priors. Irrelevant unless Cov.pedestal - is True.] - 9. rescale [boolean, scale the covariance matrix internally to avoid - possible overflow/underflow problems for long ramps. - Slightly increases computational cost. Default - True. ] - - Returns: - Class Ramp_Result holding lots of information - - """ - - if diffs2use is None: - diffs2use = np.ones(diffs.shape, np.uint8) - - if countrateguess is None: - # initial guess for count rate is the average of the unmasked - # resultant differences unless otherwise specified. - if Cov.pedestal: - countrateguess = np.sum((diffs*diffs2use)[1:], axis=0)/np.sum(diffs2use[1:], axis=0) - else: - countrateguess = np.sum((diffs*diffs2use), axis=0)/np.sum(diffs2use, axis=0) - countrateguess *= countrateguess > 0 - - # Elements of the covariance matrix - - alpha_phnoise = countrateguess*Cov.alpha_phnoise[:, np.newaxis] - alpha_readnoise = sig**2*Cov.alpha_readnoise[:, np.newaxis] - alpha = alpha_phnoise + alpha_readnoise - beta_phnoise = countrateguess*Cov.beta_phnoise[:, np.newaxis] - beta_readnoise = sig**2*Cov.beta_readnoise[:, np.newaxis] - beta = beta_phnoise + beta_readnoise - - # Mask resultant differences that should be ignored. This is half - # of what we need to do to mask these resultant differences; the - # rest comes later. - - d = diffs*diffs2use - beta = beta*diffs2use[1:]*diffs2use[:-1] - ndiffs, npix = alpha.shape - - # Rescale the covariance matrix to a determinant of 1 to - # avoid possible overflow/underflow. The uncertainty and chi - # squared value will need to be scaled back later. Note that - # theta[-1] is the determinant of the covariance matrix. - # - # The method below uses the fact that if all alpha and beta - # are multiplied by f, theta[i] is multiplied by f**i. Keep - # a running track of these factors to construct the scale at - # the end, and keep scaling throughout so that we never risk - # overflow or underflow. - - if rescale: - theta = np.ones((ndiffs + 1, npix)) - theta[1] = alpha[0] - - scale = theta[0]*1 - for i in range(2, ndiffs + 1): - - theta[i] = alpha[i - 1]/scale*theta[i - 1] - beta[i - 2]**2/scale**2*theta[i - 2] - - # Scaling every ten steps is safe for alpha up to 1e20 - # or so and incurs a negligible computational cost for - # the fractional power. - - if i % int(dn_scale) == 0 or i == ndiffs: - f = theta[i]**(1/i) - scale *= f - theta[i - 1] /= theta[i]/f - theta[i - 2] /= theta[i]/f**2 - theta[i] = 1 - else: - scale = 1 - - alpha /= scale - beta /= scale - - # All definitions and formulas here are in the paper. - - theta = np.ones((ndiffs + 1, npix)) - theta[1] = alpha[0] - for i in range(2, ndiffs + 1): - theta[i] = alpha[i - 1]*theta[i - 1] - beta[i - 2]**2*theta[i - 2] - - #print(np.mean(theta[-1]), np.std(theta[-1])) - - phi = np.ones((ndiffs + 1, npix)) - phi[ndiffs - 1] = alpha[ndiffs - 1] - for i in range(ndiffs - 2, -1, -1): - phi[i] = alpha[i]*phi[i + 1] - beta[i]**2*phi[i + 2] - - sgn = np.ones((ndiffs, npix)) - sgn[::2] = -1 - - Phi = np.zeros((ndiffs, npix)) - for i in range(ndiffs - 2, -1, -1): - Phi[i] = Phi[i + 1]*beta[i] + sgn[i + 1]*beta[i]*phi[i + 2] - - # This one is defined later in the paper and is used for jump - # detection and pedestal fitting. - - PhiD = np.zeros((ndiffs, npix)) - for i in range(ndiffs - 2, -1, -1): - PhiD[i] = (PhiD[i + 1] + sgn[i + 1]*d[i + 1]*phi[i + 2])*beta[i] - - Theta = np.zeros((ndiffs, npix)) - Theta[0] = -theta[0] - for i in range(1, ndiffs): - Theta[i] = Theta[i - 1]*beta[i - 1] + sgn[i]*theta[i] - - ThetaD = np.zeros((ndiffs + 1, npix)) - ThetaD[1] = -d[0]*theta[0] - for i in range(1, ndiffs): - ThetaD[i + 1] = beta[i - 1]*ThetaD[i] + sgn[i]*d[i]*theta[i] - - beta_extended = np.ones((ndiffs, npix)) - beta_extended[1:] = beta - - # C' and B' in the paper - - dC = sgn/theta[ndiffs]*(phi[1:]*Theta + theta[:-1]*Phi) - dC *= diffs2use - - dB = sgn/theta[ndiffs]*(phi[1:]*ThetaD[1:] + theta[:-1]*PhiD) - - # {\cal A}, {\cal B}, {\cal C} in the paper - - A = 2*np.sum(d*sgn/theta[-1]*beta_extended*phi[1:]*ThetaD[:-1], axis=0) - A += np.sum(d**2*theta[:-1]*phi[1:]/theta[ndiffs], axis=0) - - B = np.sum(d*dC, axis=0) - C = np.sum(dC, axis=0) - - r = Ramp_Result() - - # Finally, save the best-fit count rate, chi squared, uncertainty - # in the count rate, and the weights used to combine the - # resultants. - - if not Cov.pedestal: - r.countrate = B/C - r.chisq = (A - B**2/C)/scale - r.uncert = np.sqrt(scale/C) - r.weights = dC/C - r.var_poisson = np.sum(r.weights**2*alpha_phnoise, axis=0) - r.var_poisson += 2*np.sum(r.weights[1:]*r.weights[:-1]*beta_phnoise, axis=0) - r.var_rdnoise = np.sum(r.weights**2*alpha_readnoise, axis=0) - r.var_rdnoise += 2*np.sum(r.weights[1:]*r.weights[:-1]*beta_readnoise, axis=0) - - # If we are computing the pedestal, then we use the other formulas - # in the paper. - - else: - dt = Cov.mean_t[0] - Cinv_11 = theta[0]*phi[1]/theta[ndiffs] - - # Calculate the pedestal and slope using the equations in the paper. - # Do not compute weights for this case. - - b = dB[0]*C*dt - B*dC[0]*dt + dt**2*C*resetval/resetsig**2 - b /= C*Cinv_11 - dC[0]**2 + dt**2*C/resetsig**2 - a = B/C - b*dC[0]/C/dt - r.pedestal = b - r.countrate = a - r.chisq = A + a**2*C + b**2/dt**2*Cinv_11 - r.chisq += - 2*b/dt*dB[0] - 2*a*B + 2*a*b/dt*dC[0] - r.chisq /= scale - - # elements of the inverse covariance matrix - M = [C, dC[0]/dt, Cinv_11/dt**2 + 1/resetsig**2] - detM = M[0]*M[-1] - M[1]**2 - r.uncert = np.sqrt(scale*M[-1]/detM) - r.uncert_pedestal = np.sqrt(scale*M[0]/detM) - r.covar_countrate_pedestal = -scale*M[1]/detM - - # The code below computes the best chi squared, best-fit slope, - # and its uncertainty leaving out each resultant difference in - # turn. There are ndiffs possible differences that can be - # omitted. - # - # Then do it omitting two consecutive reads. There are ndiffs-1 - # possible pairs of adjacent reads that can be omitted. - # - # This approach would need to be modified if also fitting the - # pedestal, so that condition currently triggers an error. The - # modifications would make the equations significantly more - # complicated; the matrix equations to be solved by hand would be - # larger. - - if detect_jumps: - - # The algorithms below do not work if we are computing the - # pedestal here. - - if Cov.pedestal: - raise ValueError("Cannot use jump detection algorithm when fitting pedestals.") - - # Diagonal elements of the inverse covariance matrix - - Cinv_diag = theta[:-1]*phi[1:]/theta[ndiffs] - Cinv_diag *= diffs2use - - # Off-diagonal elements of the inverse covariance matrix - # one spot above and below for the case of two adjacent - # differences to be masked - - Cinv_offdiag = -beta*theta[:-2]*phi[2:]/theta[ndiffs] - - # Equations in the paper: best-fit a, b - # - # Catch warnings in case there are masked resultant - # differences, since these will be overwritten later. No need - # to warn about division by zero here. - - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - a = (Cinv_diag*B - dB*dC)/(C*Cinv_diag - dC**2) - b = (dB - a*dC)/Cinv_diag - - r.countrate_oneomit = a - r.jumpval_oneomit = b - - # Use the best-fit a, b to get chi squared - - r.chisq_oneomit = A + a**2*C - 2*a*B + b**2*Cinv_diag - 2*b*dB + 2*a*b*dC - # invert the covariance matrix of a, b to get the uncertainty on a - r.uncert_oneomit = np.sqrt(Cinv_diag/(C*Cinv_diag - dC**2)) - r.jumpsig_oneomit = np.sqrt(C/(C*Cinv_diag - dC**2)) - - r.chisq_oneomit /= scale - r.uncert_oneomit *= np.sqrt(scale) - r.jumpsig_oneomit *= np.sqrt(scale) - - # Now for two omissions in a row. This is more work. Again, - # all equations are in the paper. I first define three - # factors that will be used more than once to save a bit of - # computational effort. - - cpj_fac = dC[:-1]**2 - C*Cinv_diag[:-1] - cjck_fac = dC[:-1]*dC[1:] - C*Cinv_offdiag - bcpj_fac = B*dC[:-1] - dB[:-1]*C - - with warnings.catch_warnings(): - warnings.simplefilter('ignore') - # best-fit a, b, c - c = (bcpj_fac/cpj_fac - (B*dC[1:] - dB[1:]*C)/cjck_fac) - c /= cjck_fac/cpj_fac - (dC[1:]**2 - C*Cinv_diag[1:])/cjck_fac - b = (bcpj_fac - c*cjck_fac)/cpj_fac - a = (B - b*dC[:-1] - c*dC[1:])/C - r.countrate_twoomit = a - - # best-fit chi squared - r.chisq_twoomit = A + a**2*C + b**2*Cinv_diag[:-1] + c**2*Cinv_diag[1:] - r.chisq_twoomit -= 2*a*B + 2*b*dB[:-1] + 2*c*dB[1:] - r.chisq_twoomit += 2*a*b*dC[:-1] + 2*a*c*dC[1:] + 2*b*c*Cinv_offdiag - r.chisq_twoomit /= scale - - # uncertainty on the slope from inverting the (a, b, c) - # covariance matrix - fac = Cinv_diag[1:]*Cinv_diag[:-1] - Cinv_offdiag**2 - term2 = dC[:-1]*(dC[:-1]*Cinv_diag[1:] - Cinv_offdiag*dC[1:]) - term3 = dC[1:]*(dC[:-1]*Cinv_offdiag - Cinv_diag[:-1]*dC[1:]) - r.uncert_twoomit = np.sqrt(fac/(C*fac - term2 + term3)) - r.uncert_twoomit *= np.sqrt(scale) - - r.fill_masked_reads(diffs2use) - - return r - - -def mask_jumps(diffs, Cov, sig, threshold_oneomit=20.25, - threshold_twoomit=23.8, diffs2use=None): - - """ - - Function mask_jumps implements a likelihood-based, iterative jump - detection algorithm. - - Arguments: - 1. diffs [resultant differences] - 2. Cov [class Covar, holds the covariance matrix information. Must - be based on differences alone (i.e. without the pedestal)] - 3. sig [read noise, 1D array] - Optional arguments: - 4. threshold_oneomit [float, minimum chisq improvement to exclude - a single resultant difference. Default 20.25, - i.e., 4.5 sigma] - 5. threshold_twoomit [float, minimum chisq improvement to exclude - two sequential resultant differences. - Default 23.8, i.e., 4.5 sigma] - 6. diffs2use [a 2D array of the same shape as d, one for resultant - differences that appear ok and zero for resultant - differences flagged as contaminated. These flagged - differences will be ignored throughout jump detection, - which will only flag additional differences and - overwrite the data in this array. Default None] - - Returns: - 1. diffs2use [a 2D array of the same shape as d, one for resultant - differences that appear ok and zero for resultant - differences flagged as contaminated.] - 2. countrates [a 1D array of the count rates after masking the pixels - and resultants in diffs2use.] - - """ - - if Cov.pedestal: - raise ValueError("Cannot mask jumps with a Covar class that includes a pedestal fit.") - - # Force a copy of the input array for more efficient memory access. - - d = diffs*1 - - # We can use one-omit searches only where the reads immediately - # preceding and following have just one read. If a readout - # pattern has more than one read per resultant but significant - # gaps between resultants then a one-omit search might still be a - # good idea even with multiple-read resultants. - - oneomit_ok = Cov.Nreads[1:]*Cov.Nreads[:-1] >= 1 - oneomit_ok[0] = oneomit_ok[-1] = True - - # Other than that, we need to omit two. If a resultant has more - # than two reads, we need to omit both differences containing it - # (one pair of omissions in the differences). - - twoomit_ok = Cov.Nreads[1:-1] > 1 - - # This is the array to return: one for resultant differences to - # use, zero for resultant differences to ignore. - - if diffs2use is None: - diffs2use = np.ones(d.shape, np.uint8) - - # We need to estimate the covariance matrix. I'll use the median - # here for now to limit problems with the count rate in reads with - # jumps (which is what we are looking for) since we'll be using - # likelihoods and chi squared; getting the covariance matrix - # reasonably close to correct is important. - - countrateguess = np.median(d, axis=0)[np.newaxis, :] - countrateguess *= countrateguess > 0 - - # boolean arrays to be used later - recheck = np.ones(d.shape[1]) == 1 - dropped = np.ones(d.shape[1]) == 0 - - for j in range(d.shape[0]): - - # No need for indexing on the first pass. - if j == 0: - result = fit_ramps(d, Cov, sig, countrateguess=countrateguess, - diffs2use=diffs2use, detect_jumps=True) - # Also save the count rates so that we can use them later - # for debiasing. - countrate = result.countrate*1. - else: - result = fit_ramps(d[:, recheck], Cov, sig[recheck], - countrateguess=countrateguess[:, recheck], - diffs2use=diffs2use[:, recheck], - detect_jumps=True) - - # Chi squared improvements - - dchisq_two = result.chisq - result.chisq_twoomit - dchisq_one = result.chisq - result.chisq_oneomit - - # We want the largest chi squared difference - - best_dchisq_one = np.amax(dchisq_one*oneomit_ok[:, np.newaxis], axis=0) - best_dchisq_two = np.amax(dchisq_two*twoomit_ok[:, np.newaxis], axis=0) - - # Is the best improvement from dropping one resultant - # difference or two? Two drops will always offer more - # improvement than one so penalize them by the respective - # thresholds. Then find the chi squared improvement - # corresponding to dropping either one or two reads, whichever - # is better, if either exceeded the threshold. - - onedropbetter = (best_dchisq_one - threshold_oneomit > best_dchisq_two - threshold_twoomit) - - best_dchisq = best_dchisq_one*(best_dchisq_one > threshold_oneomit)*onedropbetter - best_dchisq += best_dchisq_two*(best_dchisq_two > threshold_twoomit)*(~onedropbetter) - - # If nothing exceeded the threshold set the improvement to - # NaN so that dchisq==best_dchisq is guaranteed to be False. - - best_dchisq[best_dchisq == 0] = np.nan - - # Now make the masks for which resultant difference(s) to - # drop, count the number of ramps affected, and drop them. - # If no ramps were affected break the loop. - - dropone = dchisq_one == best_dchisq - droptwo = dchisq_two == best_dchisq - - drop = np.any([np.sum(dropone, axis=0), - np.sum(droptwo, axis=0)], axis=0) - - if np.sum(drop) == 0: - break - - # Store the updated counts with omitted reads - - new_cts = np.zeros(np.sum(recheck)) - i_d1 = np.sum(dropone, axis=0) > 0 - new_cts[i_d1] = np.sum(result.countrate_oneomit*dropone, axis=0)[i_d1] - i_d2 = np.sum(droptwo, axis=0) > 0 - new_cts[i_d2] = np.sum(result.countrate_twoomit*droptwo, axis=0)[i_d2] - - # zero out count rates with drops and add their new values back in - - countrate[recheck] *= drop == 0 - countrate[recheck] += new_cts - - # Drop the read (set diffs2use=0) if the boolean array is True. - - diffs2use[:, recheck] *= ~dropone - diffs2use[:-1, recheck] *= ~droptwo - diffs2use[1:, recheck] *= ~droptwo - - # No need to repeat this on the entire ramp, only re-search - # ramps that had a resultant difference dropped this time. - - dropped[:] = False - dropped[recheck] = drop - recheck[:] = dropped - - # Do not try to search for bad resultants if we have already - # given up on all but one, two, or three resultant differences - # in the ramp. If there are only two left we have no way of - # choosing which one is "good". If there are three left we - # run into trouble in case we need to discard two. - - recheck[np.sum(diffs2use, axis=0) <= 3] = False - - return diffs2use, countrate - - - -def getramps(countrate, sig, readtimes, nramps=10): - - """ - Function getramps: make some synthetic ramps - Arguments: - 1. countrate [electrons/time] - 2. sig [single read read noise, electrons] - 3. readtimes [list of values or lists for the times of reads. If - a list of lists, times for reads that are averaged - together to produce a resultant.] - Optional Arguments: - 4. nramps [number of ramps desired, default 10] - - Returns: - 1. counts [electrons in each read, shape (nreads, nramps)] - """ - - t_last = 0 - counts = np.zeros((len(readtimes), nramps)) # resultant values to return - counts_total = np.zeros(nramps) # Running total of electrons - - for k in range(len(readtimes)): - t = readtimes[k] - counts_last = counts_total*1. - counts_resultant = counts_last*0 - - if hasattr(t, "__len__"): - for i in range(len(t)): - lam = countrate*(t[i] - t_last) # expected number of electrons - counts_total += np.random.poisson(lam, nramps) - counts[k] += counts_total/len(t) - t_last = t[i] - - # add read noise - counts[k] += np.random.randn(nramps)*sig/np.sqrt(len(t)) - - else: - if t_last is None: - t_last = t - lam = countrate*(t - t_last) # expected number of electrons - counts_total += np.random.poisson(lam, nramps) - counts[k] = counts_total - t_last = t - - # add read noise - counts[k] += np.random.randn(nramps)*sig - - return counts - diff --git a/spaceKLIP/imagetools.py b/spaceKLIP/imagetools.py index 52698a7b..e94a8106 100644 --- a/spaceKLIP/imagetools.py +++ b/spaceKLIP/imagetools.py @@ -588,138 +588,6 @@ def subtract_median(self, # Update spaceKLIP database. self.database.update_obs(key, j, fitsfile, maskfile) - def subtract_background_godoy(self, - types=['SCI', 'REF'], - subdir='bgsub'): - - """ - Subtract the corresponding background observations from the SCI and REF - data in the spaceKLIP database using a method developed by Nico Godoy. - - Parameters - ---------- - types : list of str - File types to run the subtraction over. - - subdir : str, optional - Name of the directory where the data products shall be saved. The - default is 'bgsub'. - - Returns - ------- - None. - - """ - - # Set output directory. - output_dir = os.path.join(self.database.output_dir, subdir) - if not os.path.exists(output_dir): - os.makedirs(output_dir) - - # Loop through concatenations. - for i, key in enumerate(self.database.obs.keys()): - - # Load in bunch of stuff - # Find science, reference, and background files. - ww_sci = np.where(self.database.obs[key]['TYPE'] == 'SCI')[0] - ww_ref = np.where(self.database.obs[key]['TYPE'] == 'REF')[0] - ww_sci_bg = np.where(self.database.obs[key]['TYPE'] == 'SCI_BG')[0] - ww_ref_bg = np.where(self.database.obs[key]['TYPE'] == 'REF_BG')[0] - - # Loop over science and reference files - for typ in types: - if typ == 'SCI': - ww, ww_bg = ww_sci, ww_sci_bg - elif typ == 'REF': - ww, ww_bg = ww_ref, ww_ref_bg - - # Gather background files. - if len(ww_bg) == 0: - raise UserWarning('Could not find any background files.') - else: - bg_data, bg_erro, bg_pxdq = [], [], [] - for j in ww_bg: - # Read background file. - fitsfile = self.database.obs[key]['FITSFILE'][j] - data, erro, pxdq, head_pri, head_sci, is2d, imshifts, maskoffs = ut.read_obs(fitsfile) - - # Compute median science background. - bg_data += [data] - bg_erro += [erro] - bg_pxdq += [pxdq] - bg_data, bg_erro, bg_pxdq = np.array(bg_data), np.array(bg_erro), np.array(bg_pxdq) - - # If multiple files, take the median. Otherwise, carry on. - if bg_data.ndim == 4: - bg_data = np.nanmedian(bg_data, axis=0) - - # Loop over individual files - for j in ww: - # Read FITS file. - fitsfile = self.database.obs[key]['FITSFILE'][j] - data, erro, pxdq, head_pri, head_sci, is2d, imshifts, maskoffs = ut.read_obs(fitsfile) - - # Subtract the background per frame - data -= bg_data - - # Loop over integrations - data_bg_sub = np.empty_like(data) - for k in range(data.shape[0]): - # Subtract median of corresponding background frame from the frame - bg_submed = bg_data[k,:,:] - np.nanmedian(bg_data[k,:,:]) - # Do the same for the data (that's already background subtracted) - data_submed = data[k,:,:] - np.nanmedian(data[k,:,:]) - - # Specify sections for initial guess - # sect1 = data_submed[108:118,12:62]/bg_submed[108:118,12:62] - # sect2 = data_submed[93:106,152:207]/bg_submed[93:106,152:207] - sect1 = data_submed[112:118,4:10]/bg_submed[112:118,4:10] - sect2 = data_submed[95:101,207:212]/bg_submed[95:101,207:212] - - # Reshape into 1d arrays and concatenate - s1 = sect1.reshape(1,sect1.shape[0]*sect1.shape[1]) - s2 = sect2.reshape(1,sect2.shape[0]*sect2.shape[1]) - s12 = np.concatenate((s1[0,:],s2[0,:])) - - # Take median of concatenated array - cte = np.nanmedian(s12) - - # Use filter to determine mask for estimating BG scaling - # at the moment only have it working for F1140C. - filt = self.database.obs[key]['FILTER'][j] - if filt not in ['F1065C', 'F1140C', 'F1550C']: - raise NotImplementedError('Godoy subtraction is only supported for MIRI FQPMs at this time!') - else: - bgmaskbase = os.path.split(os.path.abspath(__file__))[0] - bgmaskfile = os.path.join(bgmaskbase, 'resources/miri_bg_masks/godoy_mask_{}.fits'.format(filt.lower())) - - # Run minimisation function, 'res' will tell us if there is any residual - # background that wasn't removed in the initial attempt. I.e. do we - # need to subtract a little bit more or less? - res = minimize(ut.bg_minimize, - x0=cte*100, - args=(data_submed, bg_submed, bgmaskfile), - method='L-BFGS-B', - tol=1e-7) - - # Extract scale factor for the background from res - scale = res.x/100 - - # Scale the background, and now subtract this correction from the original - # background subtracted data - data_improved_bgsub = data_submed - bg_submed*scale - - # Subtract median of residual frame to remove any residual median offset - data_bg_sub[k] = data_improved_bgsub - np.nanmedian(data_improved_bgsub) - - # Write FITS file and PSF mask. - fitsfile = ut.write_obs(fitsfile, output_dir, data_bg_sub, erro, pxdq, head_pri, head_sci, is2d, imshifts, maskoffs) - - # Update spaceKLIP database. - self.database.update_obs(key, j, fitsfile) - - pass - def subtract_background(self, nints_per_med=None, subdir='bgsub'): @@ -2810,4 +2678,3 @@ def create_rec_mask(h, w, center=None, z=None): log.info(f" Plot saved in {output_file}") plt.close(fig) - diff --git a/spaceKLIP/pyklippipeline.py b/spaceKLIP/pyklippipeline.py index bac3de6b..78e37b31 100644 --- a/spaceKLIP/pyklippipeline.py +++ b/spaceKLIP/pyklippipeline.py @@ -1,4 +1,3 @@ - from __future__ import division import matplotlib diff --git a/spaceKLIP/resources/miri_bg_masks/godoy_mask_f1065c.fits b/spaceKLIP/resources/miri_bg_masks/godoy_mask_f1065c.fits deleted file mode 100644 index f95c11a2b3b8073975e55644b15d3f0b45d7c6c8..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 184320 zcmeI3&8jU|a)k9L=LL-2G9DUZV>Jlura>BLdTptp2a=2=4|-%jnjg(dH6JH#-dumR zt9I4iI|P@^jEszxl`Gdd`axfP`R{-Hw=dq_{_(B&>+83_ef!TJfB4gnKYjo6+b=)9 z{r=0Zzk2)mmk&RD`{Ac=3;wUrzkK=iAHV#j)cc#i|NK|I4?q3%;V*CBe)#3X+b@6l zukUJ=ck~bc{LPnN{Ti12h29VU{r!J@_tV>te?s}s-+%YR&wu{$ho9d?#QRO}cctFv zzxxLvPP)J8ecsUfyGaqDeeuoLU;N=;OTD+ZuZ1}2{w6lbHO&+71Uvyxz!UHUJONK& zJb{1u-_JhtpW|lmKmYgHXZ+DQgN3bd#lnZ5amOFe&+~Y#c~tKL%x~+pcEN9d-}}z* zC+55QwtOf1-ihb=yz^|mE#95){(0}-cYFNZR=>|Y_v7d88GFXw0q=l!z&qd_@D6wf zyaV0=?|^r}JK!Dg4tNK=1KxrA?129r=zZqwefFHa1Kt7efOo(<;2rP|cn7=#-U07` zcfdQ~9q@6$89&J6Fp*VTFXt!>#4 zbobEOa=vS>v(s0(3wH8e{Ty{3x_(>U@q6v^Gn~V_>3s_vzb(J5-kw}B$L^ld+umpI zt#`&ueCG4~47tKS?$rZ2VZEcC!Pz^&KDPQ8^);D_Dv zZCZ~;g>BJ(wHlq<-rvsM?5+C@dOG)&IpsQyV=g(~GOJ1hY>VEj)*Lyr(07=sJNVgBuQ=GMtuz9g)uDctZ{GZDnbpdDQ1gto+BL9hMis_fsEZh1`4vX(idD^( zx>jGS8hi7Z%ig+Y>D_WaXZ3uJlHENvmy?gG2*bw!zWLSt@v2`+PloDcsh6Ghd<)yylRbn&dNM+&L>aY%EPmr+au5E{b#k# zLidZE>*bZBwOSLt=BQCD^K(s=uQ1i*9C$fLOy1USXRjS-?G||0>MYGGNAXczat#S9^9(8ia*j+J^Wk)t`-g$aoE5# zwtwtAxC5-t8?3q$a>`%vsmJEp>dfH4$5KZ;s)Ijc(QI*ZzIAr1UU9V4MxI!?>XXN6 zO<-`zD@M-teNKDMLT8STJ~0+NELd2@v#-3&9p$oCK9_oUIR>XL$2~qZ!J#hma|{nV z%EJo_S1tL#@UvAuJjLKw7&+T|zB)5zyE?0$#wa&)mTF`^xr~ou>cT2F$K*1u^*Kk4 zTIGX-wc@~MwBezLRh%;_tNd0CoUoOKv3f4Eb8GPwZ_QpgtM_nA?dUADd%VmcR@WHpk$zmin;NfDvb_+?*puL%9_O-}Zg}tFzD<=h|K_bMGCq4|*?N{?o_eNZ{26ES*;7*N+V;{`CB-- zrWMbeSJhT^Ev!|qxYo9C<=DUPyVE-fpZP1EG2)R`4g67!=Bb&nGj%hZ*|TPNEspk) z{Z>tdQ=4<(w(t9V#w_x^L`)y#h+9@`G6#){)5vL|8Kf~$OM{ZT4;)7LobiU@P zuP}19?{{|ZS?Imsuk=6-Hfr(1R$7|h)6ZvD+B2N04s4}fzn%K9idVjJTfWx8QOy2$ z&PU8h?*d&(p2xaViR9(R-rtrn-3i<(is(xoOWe=-kFg;zOrRxSFgbzAwIhuy04 zN_EVtXLTMu-fV6!z0jFGGxM}o`%%2AQEn@r@rbo<51%vk(6Yw<_Y*Vaj(V2aXAbqP zeCBJsr6$*Fyz1xP=?A>h>v719>flvYYpZ-<%E51OjUV%7R_Dq+t=N3VZ0*)M^rPPvbymDX7*V$D5E%T`K6hm9Fkxz}Xi={OlvBQ z7MAfGqZ>GX^4Z{LTY9}}e^;81&Px5QJ*q~_SMAIpr+8o*k9;(!Bc?9LfrC*4mz>5e z|CUT0v$Ubpl`?cO$e*;{i(Z{_?RcgsvyYxKKTYt29Io$9^d&(0fu*insg zXZ+Bc;a!DO^*fr$9{K$JduXag){j{8&+hB@SKq36_2DWOJ~+0O&v+H15w^ zM?Q-f+vMtFuu4p1#h%(j{Nzu7=So_Ne~Q zS37%_)~x!Rv08Tze)N{lsj$pBB4$o%d*>)`)(^cFw|CY|uBCe>p8Gz!m)?9^XXhS! zacK1NSI?c*e{1%NGdsoh^wHnb?%Cg?X~i6K&&(E_c*bt2p@$j$y*Sb5URvnPetYH- zpKbY$J?+`wbk_MgLwcw%a%z=l|NWiDjC%9v%}{eaf26PZN=NfEex>mWSe}i3Gdac7 zCm$AEF~wVX_~gTa+19g|nUT&D{~q(?tRrmne}wK9O=j7uC9f8)V)*1*7M@z=Dqi_; ztK8MG=>6(F9^D`Jr}l_5RsRe>e^b>%kCSr|R~m{{HDK_F_p+h8m)2IDEA-IneWo|= zW`^542b~sY51ca|!N&Kt`qWwUIOLRH`S8M4`A5X~rhE0ydD-6U{Ck|!`?j2WlsL~V z_Z@k~GJh1`QajRi{+ORvdS~I^GUu#S#UuTyrovRG%I}5gIpJ5jnKOz*W5lWQoL%LN zZJX(p=EnW5+Ldou)!v)m>$P=OuO{cv$hhOje;>Ep;mTgX#m<}y_Y8Y$e>3&L(Qlir zG4Gsrr+20PnfIZ6T+v^d%}W0x-AZdzYj72_TmUkziCf>)t&XL^P{%HSFzb1 z%-jEMaaU)~mLBeDtMA=6d$pPCc>nj)nyKDPi_^oY&m8b8*j8^DkNB2; zD?P=wXluRw`@M@jX07)24|h` zbKieojb;&d)#GrNBiqx6+^GH@`7J$;YMk%y&m+&L+R2sgp}wjihgA+3R&&{hbGFTW zRq>V{E1emPvu3#Fob!gcqp+)Y(5h9urQ74=+}ZKWTvmFV(aW4lKXV+f@9>tHaTl|- z!Jp+Sk8^rf`Qi7tSI+SrvpvGUlBM_D6L@7a*YRI>cxDbd&l!H^XnpwPGycdi-=tQ* zRq?IeZ~~IAgW$p8VE*%=QjX?P~wMTFv2SeQ>pK6g$(WAGO+#;<82qfhKeA08Gyv08F1t6Y5Lf?+imKKRwG_Px^UYMt`0;K#S# zyN@%@Y<(kl-OA74Xk?pfSMxo7&b>;UXF7`WT$E2vv60`Z34Bz0bp8rG=r^`zrJR;e zE#IUue$?VuR<*1AGh+J9d-;3Jsp_Sro4LWu*i4N6h&!V{!x^2kcYf4UE8oK-=RfTk z9C3&Do_BnMp1I|lX6oUL^n+J%T6o1h1Izv8IQJC!EB(28f1`Qk{&Rffto5wcIm+4M z_rAZE)3auIT5}a&dsn~eFY@|FIf~!Q*Ym@_*STBHndxE$0i$4nmm zJ8tctIbcH%-}bgy zD7Mv)^S0Ky?tOOLTQw7MEnC$nrgeemIAfZBbzkS$nv-(!yfioK@1Ze!)>V1TU~7NW zXPfnMo;dSX;;cce7CpsY#pl~wmO11+xZ2E_iKBjo_sATbQT|@dt$d@_sh>5-1&)u+ z+$wHyx8l@2yRBxqg5~+F)V>mzIqmt)XksgEa$&3d2-AAiP`uL6{Hy!be2i_|$E(aO z@APV4xt=(5!msq;RF=9b_v~0TE3o|S;8^m+Sn#l5m4#R3jP+-i=ff;VGtYGztA_f@ zQC~SNKlrWOz3_YXp!31k@2s%MgNIctC zTGR5k)U>n|M^~|X`EfR9(MPuAv04KbmO91MC%=b1N_UUmUR9%YCly{{%uRFQxA4rX z;#+Fb3k`hi493}++roR=k)sZtT5@|?wD;1p&ino4SukttO1`S;VdTQ*b0XIK$Y14W zJ1`|@iPw`OP&}Do-KK7D-RA28#q3;!pP0q z)_GCCmwwetPn({w^g?WVjdvCIFmpG0|0{Q=bu;Hw+LgvuOzSvHvF!Iak<0uoai!Bc zmownlUgs^{>n_*BiFN49lcdKw3K82(>=;F ziMeI$5#s8;!A4fUPr2YTm~vM9IHSj}dKmG0`I#JNVsj0+?RDQuzO#r|@2jOz;VmrA zY~?DP9wVoI#T{Xz`cb~p-%^84)T*y>#yKb3SzfMl{Ju|HXTj_%yQ-NrsLyySUgd`l zHELV>J$zI@%8#_kqZ2WHSn}8wmT}dneDv~p71!9ucXI1IboT12F1TWq4-X56cw~dG zoE&Q|eDo`uxf*AF)>1BhRikzAW^C4qc%=?*rDLqW7iY~U?ku0pJz{fSV>qg#*vJnp z<>Xj%;pek8#_!d@39EXuHMu_55yK-63x?Gke&q16@NGYx;iJujI}N+iJ9SkJT5wvH z+Frhu3!M6FwT4{g!y!LosmWGr^Eu>3HRR!6$&Pw!ee;#R;gNc&=GnRfaJBhcH5Wel zjOQ4hTIB>zTyYd@`P7G1P77B~ZXHOJh- zG6TibCqJ@U)AB0~d0=YED;7TZt+u7x!sEPNj(*?} z+unX>%rq?XRIKvh^{m$Qc)grzgdbeRTD5Rn`Bl7irsDW$6Jx<|wdi59-Ux5yvyS8U z_po&q%sMO@U^!L}xg#xoVB;)&ESUD8ImNg7TlA<2yY<{w-7VTj^wa9o_@{TD=g(ap zIjekToW;4=5k`N?2g|nNR9LQCiTBo9IosS%_xg%6i5;OY_C+yrTP>Qo9vmBae5_*B zYYq$^F`TgEG6vW7{r)|Rk9Tlo9=c!Zhla-Cujux;=)yr8yk~P?XjB$may{M5rN(x@ zdw=9vj`oAy;v0!88yv-{$##a5YjTW6;P@>|EqUU~21jvfz%>R_zWtBq=Xtc}!7k~H zz_Nt{POO%kTFuY;)CV6Qn=v#akKeMYp*qTkGoz<<_Q&&or5UkD)i;xy-7)27ztZnv z)CUJ2o3WnGOfGb1c;;6)$70b9(q58a@3-&pq;duF!j{2XHKVl4J0e4gQv#*5osZvC3WX zZ>?WB>k+w+FrTVNu&`I1`<5Q5;hS;{zGAZmF}9_3YrduXD)?xASC~WenERmbDwnxu z#`HPU16=*43Y*bV4jN#@Sn$dQr^1bm=Q1-XML$M>iWy&U(>2yT@} zgB&py+_s*d=W^~`ZrxMr)yAFXc#nJTY4&n^oNbMB_w*6x&raUbI~rKDatwZinB6V*S8KCBtCj1v#Q7}x z%lH;;*Z7UxYX@|1ud=)Jbgvos9KW|Cc7Q#`9&v^>y}Pft1H7rQM|?l}M`%}?M}0?i zma$`Y=9%Vt)em; z;&Jc7zc(KD?#_C5?2Lad_jkwh>Cfk^{k`hB`mVe?-kq!N>7(`KcjdSAXgm8T{p#C# zl%8Ly&+6|ZD(sbJ;&1iugIjl?^>YvWR(;odM_+XZctiR36!@#o#oz6p3AgS5Z^&=U zyL9U=x&GrF^4t1&9^MbnW6vGneW_-JY2Lnf&^s{Rf$AMqm}fbjrDx9`_bDE=eF$!?<%Zk=oxqi z_S*s8mS?)(Og&4_GIqefue=BTS>PS;4tNK=1Kt7efOo(<;2rP|cn7=#-U07`cfdQ~ z9q+dsV(|9$=T_iz9C;}3uN@u%;9e*5Lew?BOO z^;d5{|MKC7Z$JF>ZGrza^sirj{pT;gDf#~HA3y(<@54_&efZnkw;z7_@b=5!{_DG% zo5NFZzbQ`+t&h|bblk8#2V%ecmv*mH{cC;1Kxl)Fy6pF z|Lf_@Dp#bXO6p;A~;4-+S{7wZ5$!_wM)Gy?6QC zyY#p9-0$n}_qzA>w{_j_kI+5uc?7pPmutoEv3zHD)EVy>_Kx$t>P&u*$JL{LQ~0l{ zW8eF6b$rLS<9kwj$2mQdb^mA@&~GaIN1NaC`pv8XYrq~|RH8J+GuSH)Y!*7SuZ?3#MlcyC8*prz5R znpIxq()iJLGynGQ$M3anhq|}D^?xL{znOLKy=HH%cU;5uC^dTLr88@-`tVdSa!YK* zLu?6;vFCkDtzEy$yKzr!YagvRb=d~KinVw~apf9eTlqaWTHKCb|IOa*4ODO8ZK zYrSg~uVTHo>EUyZGq%UOGzX{Jr5s1&W45DvnE(16s|NMWHNNznlqYI4m!4bLB|ViR zb8Bu1r?J;-^U>cRYmxhCoLJbAR;i_Spwo z!?IV-VSBO2DHgw~feVaU#p7?W^imG(SLMNi<81q`0q@S~Nu^gw4H{gt| z@L>0{+Dkc!V>M@gUNhd+r@Lxtia3_t=@}mjCM>bA8qYcMiYL~q%emgU$E`UyTK&n} ze)^`~$)`K1w8fr!xe;bmTji@5eZr3NidA3h!q=Soiqp9L_RhSU_HIUc8s%fRV06#K z!fG79%2zS)VMqCl<*aq#6VFz0nj@|;`|Uk$*-h(xdT&6kRXaP+79MJIuT~t~C=X68 zHCWBz=UhI6QyUmPXCDr=f#IC-oXKOWcrQloRc-8T>)qMC+ThePzV;#(R^u%ne9jZc z=6vf69;>)2wj{4u_4nW%acIBr!K(!mG5l;5uQ-iUQ^m-=s@-EJ%7M>5IjbcXmbk_$ z9~?Ggh=bV*>-e%E6IyxrQ_N%FbZ44m>vJ^qhTa!jdD#IV>^7Y5vuG z-G}cfd8m(Rp|wpCNb zt5}cwnR47c+twX>xL&TswS=@-{>U?z%Vqs_J(;EDlTClm7o?Dn1{uzw& zkNg(js4nuY+A3c0RXE3%?e)xk=^N6wijiwst-)_$M{(-I23F7bnS*$?8JEu)KRW9? zaE)rngYQ|eGdA;3ljA-g`~9(ZEc@mhmRMNg8teJs)%NhMoYs$gIA(EyA7M40@w5GA zIN-`wYe&A;Dpuq6Ki2gT`_SFv(K1AJKJ)p+GA zcI0b4SdH1Qb-4Fl*uC1NJ+wA_>eVNvHrFjZ2Uk@?k1D6JS)A6WPkl8<6^q(djmAel zGgLiBwR`7vmc6+jkx%!5pKa@2iVKXMdw#{&!;fK5;C0 z&W~uTdp+Vk_L!OX9NFB99co=*^gQw_&QVP(zqEJ7qZ~_SQ*4V@xxr}No<96x_vAI# zvVH11^nCW3F%c)GJ~3=72X8A+{n3^=q7S~>UjG)Z6>s%WeD<|}&gofkYx*dc{aOpN zd&zIH=A3&i<8k>nzht2iM0p{+OJKTNFqvqrOExV;}i!;aRk>hc0 zJ2fM#HTa{NSIy7v2b|h1J5wI~ zEqq!F-^glB_#=Ge_r}fiY+;tnGQviF=1-nwrskC~Xy)u4=sWG%TX}mvJ$sm(i<%s} z;=E*z(LKo@ZS$Uv?q~i_e=m81t-FC|Pa8As$#)i4HTTy2F0IdZU0QGar{Buboosmn ztsV7M>?`=Aee9X2l#YJ3g z)arSOpZn(X5-wvueutTN$@{Kq6r(%VxcXo;*7BENTO4P_TeC)d$qX$#{8e6aBY$?5 z5s&j9e@|=oq3_Y+2u#jt9X=S&BMUz42zzy&S#^%3bGPsncl04XvO1I2<+$cXoX*eh z|A<}Xx6>Xur{}=sv)1OkVz2b`T=~5BYhcTwr{j2_yTej3YuU~=rpbJeek9XZ1cH@m}*YQP;~bMICU z#kDX;$76=0dj!AU6@N)zYSvnMt*t++uC=SP`sY6EZDgq%;kM@KQ`z9I;u_o9f42UL zp0l%MoX$GqNBs=$5ip~>%;ZNH_K)wyIHFm4DWl-aCBmsc~X0t625HYs`N8_lA98byvh@ ztk#|3&&>3wn7qHG_n=BJpL#vt#pnKandhoG*n`@cS@Jt;jdI<|*BQgt*>W!D6{~(#mup)X zjR&WmNB&a$h;K`NW*=MlTeb9UmT8~>j!-^1N0r!$AIvuUmsS4_{}THoV#ZrS&E_P+G(xA@Gxg=1?iGhlmsU}h{l zGuR5VwZ7tX?B3rM@3Z9X!B#f-s`!yH?$ENm-mRJ(SB{+X`7tye@1*aj{*iAOeP;M| zo~!-6=3T=xQ)A4Wrk;6ESAQd;{!4d1;%ntweJXql>-aU_6K`~v{WH&9?V3HEb(bx( zt*LcQ{mth0-fK^p_gT(c_OoUO+HX{^*d>1OF2Nf+x3AHA9g$!2Hb=~{hI2Ig8qTO0 z@jCxlJL>76_m;i7nOu*%m7C$c3x+%BJj9N)+7JK8e&(*fcXnUdnVnU2GqD+N#msPC zjnR3p_ItaydbVfpue-1A8u==&G36NfTG#VOJXhv*Zk;ZYvEgy$2js9l5ZhdOXTIowdR{R_o%>12Do1Pvld(}l?3FCLx|@wRxMHTO`u6UsdhFqPIpx^m z(_;(Av(&N+wR$&=hp##Ht6Z+F&dPzWxGJXk89z8?a3jnfvy6Id$&c>Y%D3)C-1bx5 zZ+)k$cFNrJ=e*QXAGTKuu7`=-49|$;>U?#FD#o1b=B}1G!=5#x_FUtiHRnv9HS=pt z^w4wmM|tqf$~i2tSq`ne%GX($zcmZ_%I5o5n4@Dl_tCwSU;9>mV5@k$3H@7iPbhi=#}=aF5CAMw3vo;_xp_mT6@cGmko+ugX2 z-Z6aQVXJe-a|X}(XnS_Qz3!v?#cy$wxBXa`ciKnm?Q*;~bH?hKI*lFapWScoyD-1n zsz)oYar@VPch$QQn|&Mb*r*F%>&ZQ$_Rl#r1A^S3ZBNuRHH!>pN<_`ZIS{eDFM$Y%5MaESSA4+xx(ZaXX-hRY{efn6;Cggc}CcbXUGr3%o z&&m;g&TAj_sjFhj15RzLMltHQ@KwBuDHp!t!q+pHu&p>bYHe%noVgc$x8D2A;J{S# zlk4?}I&c+-=9PzwhDwij>J64ySNM=|!F_xsB49_LCISp!-}u8L)@e9kp`&VJQr z4MuBbI2Escz6)Gf#fDGZ_W5?fTkG2r8(FOhU-9HL7JiQBbH4iHYj`t{a^zg=Oie3a#dA$QS2evD^Cc9(S&5 zVKv_JnWvR^{NuM+vk!Kh?H*ikWbWQMYBZnEn!~40w&Z(OYrs`8a$z+eewEW$_O&K+ zan3by);jVPo|x^YZ*}HQcq8oE-qfpwOJhC1!bhxE6FKTC+p4Xw^wM0e*Idufb$X_b z7-zNQ)J9&<_*!SbuZO$rhqnv6btdI%`Sj3y#PpnfYWA|huUy(s<0D_|vOlVge6EQz zbvZ};nXTUWneVyx-E?p2YaflP4@P6z2UFR=YP>qbr?{m)9GSnzqcv4s73=X`k<&Z0 z`h#iZBMwGw5C3Yp*h}lqy0cjvT)8%%!DXws$mMe@mU;3So>n~aIk)y)^&+-spV55u z%sAp|dzhT#to8WhUfr^fuz44n!_S!BnVObu)ha$P`P_=N_=v;9S#9L>-16b*;aiwi zywz(-yv5y$2e;<%tDe1F#Zz&-T8y1xxj*NY1=q8!TH?Lh$OR4zHKQ2#mW8L4pNYev zvuK^>Ge&b-vxnc~ZsjzN-2?sS{M|EIFxL41!C*A36?9AQiOw6h=W+A8N93!{IuH|IrExmiS z@UNY($L-wscW1OC-bZm&O!IsAy&mMS(GNc?{fV{C#BJYq=3U%t7pmvI-X?z+@Df)W zxqQ|bwZ!yX`EY5DSXhnQ_deEYV5B89x0b(E*W)FJ4Gw&)##%jkaph62i1+H8^L}PD zfWES&R^!AfyN4d>75wS;i=grD&vol!T!W!`*-hd5_= z^-K)xjLkLhX)UpgRg8UqQ?nXSJuTH#uH`)2TXDITb69XyjGXQJp1h0u?Lu`tdLPw! z<&Wxf9uDH1!}en2Y~T0fT|DkCRQD}^&HIkNBYIWbRj!IvwWHWxJfl9&KkhfN*V~~r zwQv*Tyu~izs<6bu?lH$!o-H%DcH0|1%KlWdS#S9~y3_8+3us-Rh{cSHLm!jK0I1SO#5lxe*4?>l3pt7yRAmbLyWT0i@q{BqNqSI*ULdd~Bg?=@5E&RC5e#S++D;+~@Io`jNLgrQKbjvcqKv{ z`}Dot_m%*GJjG&>hlhQ?`ReU||M|bZczOA^m*l^9FTZ>F-;W=@`}pbmA6|a^`0}T> z@4kBZ;m7wMzIp%Yn}q)}^dE2E{rT93-v9Xi<;S1?`faZA zihlX;U%&n8=VRJm=zaLV@Bio9PcI+8gYsX#|MtTVfBE>~hgT8ue$o3Z)%*OjeBGEaN{Tws%G&`6HZX zoc-!MW6rRv^{qGecWCcBAU|w1`&O*VF<#HJnhn2&XFfl_YOcmTOa8XJ58emwfOo(< z;2n6x9q_*$z1OSun*ZJOtb3(1$mWn?tKML>&XsyyefJSPy!Qxq4+nK%_h?-)7xz}n zci|Z0XY{Ys>%Hc8s~K^f@Q=DL>a5U*7VI8PXsBgPjxl~VpZg5|Z2jH3AI#&{b-CWr zJwQKu?^cX^44=7~ul(LQ(*VY~z``=#^O4u%t%{=_TWDt(VvlI~zTrQ^9rd20Gf-Z3 z-ZH;TCot8@u;5fUVy(E=*n5szqyO)%b@u#h4rlLqroZfdfw4x$;W+$^rx>_1RyE>0 zn$Nje&Wu);zXd~2oC^$=V~^_Nu54G{HFQ_$agXd*`pi+Eb2NrCbHWGSvYK<`S^GG$ zmbK^iSl>f>MR93YsUf`o#<;A(Bn`vhnj6sjFQ6ITm>Wl>k;_`_tM15|^TQv>0f%#eJ)@nS+YvM3w`ln*YFGMO zwJSZ~tP8BN%niTd0arfGh=&C>YgvQ+5uH7CNAxq!Y>%17BVs(4mfyPj=zCWBE$vpU z(gTiMoGS|++lrgvoVUWyY6P}Id!~74<;)*dpe2veXRg2H8 z!#MlQnekZz{!D(43r)4Xyj#V!uQ<=t=KPtwN8s|^GoSCMagJp^^RxIY#(MBs1DIOo zZ1Gp?aSy7EGjr^>zvr0Qt!IRHVw?5AarK$siZO4*53Lbi#z7w%0}iY4J!fc!-J_$q ztP5W`KIb#Xyx1N+E3TFKRsD<`81t3Gu@Rs1NBE;0*43UzHCjB@S02Y^YwyX`KA!EY zv%~p_IY4{FR`s*oj2pZRZ(*w4;%1j(H&e1z%J!G5;V{QBOj#cvmzIt}h9ks3c z&}eaHV$fqdlYcaq`+`4G6I`{(RlIdZ&K_Io@6n80)!nPNH(&eF^Im<r(TD;#-Ydz>l`oVE7)(9G;g?Y&wn^0PYzwn978$m(UhBeAT93Y+O; zt45UoXsNy3yG+?V7#9=EQ3IVUnmSn%_Eyb{0oXjkWm7#ktBsyt3$_vYN*{jq#az@O$y#frFea78+r<rX2ElZ-^&?| zLyv9b@|pcsjCnIYdL7LHf5wh#2Tyf1H;XHVGjeg>VpnKvsryVhm~q(b-W;zunx{F) z%kU#H?Q8bTy=gATM*LpB#;Y^esyID-g+{N2G@JUug){fx&L zbI;~;pOvpYE05!dF_!t@vCVjlG4~ce*AD*{HLqSnXV#iq#$}9W!8mIwezr5;0e-7T zaKdIh#IqPMaM&W&vdo9CeBfKGYE|bg+Rl0XeX8clIM0jm@DWq|9_Q$#RTsUk;`Q!7 zT2FiCyk4)pT-BP@+oGScx9A#k^DpXUzxRvsGS&_w+c9 zp4BpEug}`ETAZ8tm0rd*_S@fq6*C!m1LkN=j;YUFjqy2ne4g>M+zgx13@-P+LNBYY zHTazIfiZ8T5B;z!bst5G=hK^!a{0{j7{QstKI??fu@(Lin$g!x-&=do8hq}}hPktT zPkT>%wwJ5$tC>{&SNy;uCyPTXi)9%2na%jKSk>pO4)f6~Fg*hgyQMz#D$DvUpYvP1 zTgJGjNAjx~R(GuUt8vs{t$T!C^jhg^OwUywaLhq1;#ar5dujZPXMKa{Q87JN_fVeC z?OEy^x6F$*w&b5JhwlTw>JjlTYS|!@|^WIQV)3<7qON_UW?<{{eC;}KCt{2z(y_Sj`}@q*2%aW4?p8v z1=G7!=juKc|42M~Ia2?s+WbbXJ1}p=&$LGJo|%))EW<`-l;s_TML$PrUsV@xG}1q~ zw#WQk9C-t9_zi&B{`32e*_>e(Gkx%T+P~lTh1s8Bf6my$%ues^A8&!0dLGf|oEbmU zp2=UuBWlvoDYH2BV^xV?S zVwrBnt1#qOapU#(huMYA=EU*LAI-^lz%tz#6)}B-xI;Z&c@Oytvm9J?Xr9HF|uGvlio-&s_M-3yjZK_i@Lq_K5qpdeM6N z9c#a*nby!X!>gWvBM0YuTl5oFdqGTb zoDr{V#RIOk(lNg7d&A6IGX*EJRbIwrEX(6~=4*~}G#+Q=ocY}97k5MdIQK00neF8| zzVAE3%tz)3F54=PIqGLQ%E`{HyjCpJIwKZ$nCUCKvtlhh<@&tuTjKfY+&NbHGdaqQ zbEcc|GAzq&VVTC6vAD}@f35rX^m@6D=leP`7k$T(oEeVlKU;>4Y1 z`x%*oIorMwyv6MNkOw=P4}RpsXA9k)Mb0QM_!$>EYBTOxvFdKEe)YaS=lgi;x!}!O zb%RsIk)Lrm*78|wj>a>s41-=33%rUUH?t#o;Pfo(_qdL~?w!s8b5wlfj9l0%2XVH{Wn6vm z6xTC&VKv6*Q6A^?d}uN+!w~QF2;PWA9oSwDc$LlaN4OQoc>V2Fv*9kP({o$T1PAi6d+tN2-{p}7t!c4FOt-F9z*_m8$E3FlH)R@&lU#;AV!#Mnr`rzKvqR)HwaL?Mf zKkSNLpvAVO%RKlkK4Rd)0uKvJEjWsgo?F_P2I}?lm^b2QwGg z7rn5Z$ zIljWrY5-fIzgGjjusy4$_QgG`53b^T1}C%3W1qQ<_k844JjBBWhtI$=TzPxXS%0dT zaYiuJ*?Vr)+^fAJpZi;(#~Pu{=hZ&& znUndjIetqY^gk`u+JX z;q!f3vEXT3^I1#%EH}%qco}!JMy09wde%IR2IpjcONaRtk2&F2yug)7IsxV z>)daacu#!;;FVRr&-Z&n&;EL|_jh%@H}iMqH}ej72fPFQ9pJxjo>zZfD`te5uh6o_ zZQg9lobk@ETQrPwUD&Hzf*V@&vyN6Z#7fAAuQJg z2Fo$~$Gth`&DQfiZr(@DmftUYuG5M!?`U89I=bE&^v^Q;GtA=3J@7m7Js7{|exncW z?HPRB`<}CZhB+U3M}9+o#|l5P7i<6e-wvL~wdc|LhB)STzVwE6@J!_x+jn34c-)a0ejlZoky`g&pUGZGk&9yi8 zv&OLxXYKDwePNccSL)rpuKz8-9)$IL_szhw+&90i`*+se#XfAkm-Alt$+LW&#jW-N zGe6_^h2Cy83)i{*{yp-{vA<{hwjO!zzT0zr$N0Yd-p}#gerLb)=nnX|!MifLE9TAZ ziodP7d3nD)FYkbNz&qd_@D6wfyaV0=?|^r}JK!Dg4tNK=1Kt7efOo(<;2rP|cn7=# z-U07`cfdQ~9q Date: Thu, 27 Jun 2024 11:20:18 -0400 Subject: [PATCH 06/11] Minor changes --- spaceKLIP/coron1pipeline.py | 2 ++ spaceKLIP/imagetools.py | 6 ++---- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/spaceKLIP/coron1pipeline.py b/spaceKLIP/coron1pipeline.py index f13e171b..c70ffb8f 100644 --- a/spaceKLIP/coron1pipeline.py +++ b/spaceKLIP/coron1pipeline.py @@ -1,3 +1,5 @@ +from __future__ import division + import matplotlib # ============================================================================= diff --git a/spaceKLIP/imagetools.py b/spaceKLIP/imagetools.py index e94a8106..93350113 100644 --- a/spaceKLIP/imagetools.py +++ b/spaceKLIP/imagetools.py @@ -1270,7 +1270,8 @@ def find_bad_pixels_gradient(self, pxdq, key, gradient_kwargs={}): - + print('') + log.info(' --> Warning!: This routine has not been thoroughly tested and requires further development') # Check input. if 'sigma' not in gradient_kwargs.keys(): gradient_kwargs['sigma'] = 0.5 @@ -1338,9 +1339,6 @@ def find_bad_pixels_gradient(self, image[bad_pixels] = np.nan - # plt.imshow(image_to_gradient) - # plt.show() - # Flag DQ array ww[i] = ww[i] | bad_pixels pxdq[i][ww[i]] = 1 From a8886a821362bd8feb49639a8f9f4d0e2d56f76a Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Thu, 27 Jun 2024 11:30:24 -0400 Subject: [PATCH 07/11] Updating scipy.integrate.simps to scipy.integrate.simpson --- spaceKLIP/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spaceKLIP/utils.py b/spaceKLIP/utils.py index 561da6cf..3f454a06 100644 --- a/spaceKLIP/utils.py +++ b/spaceKLIP/utils.py @@ -16,7 +16,7 @@ import importlib import scipy.ndimage.interpolation as sinterp -from scipy.integrate import simps +from scipy.integrate import simpson from scipy.ndimage import fourier_shift, gaussian_filter from scipy.ndimage import shift as spline_shift From 4dcd6c3ef49ebaa2db42b3601bf4a29051b2ebd0 Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Thu, 27 Jun 2024 11:32:43 -0400 Subject: [PATCH 08/11] Updating scipy.integrate.simps to scipy.integrate.simpson part 2 --- spaceKLIP/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/spaceKLIP/utils.py b/spaceKLIP/utils.py index 3f454a06..8f8f9a74 100644 --- a/spaceKLIP/utils.py +++ b/spaceKLIP/utils.py @@ -654,8 +654,8 @@ def _get_tp_comsubst(instrume, # Compute COM substrate transmission averaged over the respective # filter profile. bandpass_throughput = np.interp(comsubst_wave, bandpass_wave, bandpass_throughput) - int_tp_bandpass = simps(bandpass_throughput, comsubst_wave) - int_tp_bandpass_comsubst = simps(bandpass_throughput * comsubst_throughput, comsubst_wave) + int_tp_bandpass = simpson(bandpass_throughput, comsubst_wave) + int_tp_bandpass_comsubst = simpson(bandpass_throughput * comsubst_throughput, comsubst_wave) tp_comsubst = int_tp_bandpass_comsubst / int_tp_bandpass # Return. From 0ea9999009f1da7004a95fd1ca622862e533b024 Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Tue, 9 Jul 2024 16:29:35 -0700 Subject: [PATCH 09/11] Adding fix_bad_pixels() function back to facilitate backwards compatability --- spaceKLIP/imagetools.py | 167 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) diff --git a/spaceKLIP/imagetools.py b/spaceKLIP/imagetools.py index 93350113..674f0b8f 100644 --- a/spaceKLIP/imagetools.py +++ b/spaceKLIP/imagetools.py @@ -897,6 +897,173 @@ def find_bad_pixels(self, pass def fix_bad_pixels(self, + method='timemed+localmed+medfilt', + sigclip_kwargs={}, + custom_kwargs={}, + timemed_kwargs={}, + localmed_kwargs={}, + medfilt_kwargs={}, + types=['SCI', 'SCI_TA', 'SCI_BG', 'REF', 'REF_TA', 'REF_BG'], + subdir='bpcleaned', + restrict_to=None): + """ + **** DEPRECATED BY FIND_BAD_PIXELS() AND CLEAN_BAD_PIXELS() **** + Identify and fix bad pixels. + + Parameters + ---------- + method : str, optional + Sequence of bad pixel identification and cleaning methods to be run + on the data. Different methods must be joined by a '+' sign without + whitespace. Available methods are: + + - sigclip: use sigma clipping to identify additional bad pixels. + - custom: use a custom bad pixel map. + - timemed: replace pixels which are only bad in some frames with + their median value from the good frames. + - localmed: replace bad pixels with the median value of their + surrounding good pixels. + - medfilt: replace bad pixels with an image plane median filter. + + The default is 'timemed+localmed+medfilt'. + sigclip_kwargs : dict, optional + Keyword arguments for the 'sigclip' method. Available keywords are: + + - sigclip : float, optional + Sigma clipping threshold. The default is 5. + - shift_x : list of int, optional + Pixels in x-direction to which each pixel shall be compared to. + The default is [-1, 0, 1]. + - shift_y : list of int, optional + Pixels in y-direction to which each pixel shall be compared to. + The default is [-1, 0, 1]. + + The default is {}. + custom_kwargs : dict, optional + Keyword arguments for the 'custom' method. The dictionary keys must + match the keys of the observations database and the dictionary + content must be binary bad pixel maps (1 = bad, 0 = good) with the + same shape as the corresponding data. The default is {}. + timemed_kwargs : dict, optional + Keyword arguments for the 'timemed' method. Available keywords are: + + - n/a + + The default is {}. + localmed_kwargs : dict, optional + Keyword arguments for the 'localmed' method. Available keywords are: + + - shift_x : list of int, optional + Pixels in x-direction from which the median shall be computed. + The default is [-1, 0, 1]. + - shift_y : list of int, optional + Pixels in y-direction from which the median shall be computed. + The default is [-1, 0, 1]. + + The default is {}. + medfilt_kwargs : dict, optional + Keyword arguments for the 'medfilt' method. Available keywords are: + + - size : int, optional + Kernel size of the median filter to be used. The default is 4. + + The default is {}. + types : list of str, optional + List of data types for which bad pixels shall be identified and + fixed. The default is ['SCI', 'SCI_TA', 'SCI_BG', 'REF', 'REF_TA', + 'REF_BG']. + subdir : str, optional + Name of the directory where the data products shall be saved. The + default is 'bpcleaned'. + + Returns + ------- + None. + + """ + log.info('--> WARNING! The fix_bad_pixels() routine is deprecated, the ..........') + log.info('--> WARNING! find_bad_pixels() and clean_bad_pixels() are preferred!!!!') + # Set output directory. + output_dir = os.path.join(self.database.output_dir, subdir) + if not os.path.exists(output_dir): + os.makedirs(output_dir) + + # Loop through concatenations. + for i, key in enumerate(self.database.obs.keys()): + # if we limit to only processing some concatenations, check whether this concatenation matches the pattern + if (restrict_to is not None) and (restrict_to not in key): + continue + + log.info('--> Concatenation ' + key) + + # Loop through FITS files. + nfitsfiles = len(self.database.obs[key]) + for j in range(nfitsfiles): + + # Read FITS file and PSF mask. + fitsfile = self.database.obs[key]['FITSFILE'][j] + data, erro, pxdq, head_pri, head_sci, is2d, imshifts, maskoffs = ut.read_obs(fitsfile) + maskfile = self.database.obs[key]['MASKFILE'][j] + mask = ut.read_msk(maskfile) + + # Skip file types that are not in the list of types. + if self.database.obs[key]['TYPE'][j] in types: + + # Call bad pixel cleaning routines. + pxdq_temp = pxdq.copy() + # if self.database.obs[key]['TELESCOP'][j] == 'JWST' and self.database.obs[key]['INSTRUME'][j] == 'NIRCAM': + # pxdq_temp = (pxdq_temp != 0) & np.logical_not(pxdq_temp & 512 == 512) + # else: + pxdq_temp = (np.isnan(data) | (pxdq_temp & 1 == 1)) & np.logical_not(pxdq_temp & 512 == 512) + method_split = method.split('+') + for k in range(len(method_split)): + head, tail = os.path.split(fitsfile) + if method_split[k] == 'sigclip': + log.info(' --> Method ' + method_split[k] + ': ' + tail) + self.find_bad_pixels_sigclip(data, erro, pxdq_temp, pxdq & 512 == 512, sigclip_kwargs) + elif method_split[k] == 'custom': + log.info(' --> Method ' + method_split[k] + ': ' + tail) + if self.database.obs[key]['TYPE'][j] not in ['SCI_TA', 'REF_TA']: + self.find_bad_pixels_custom(data, erro, pxdq_temp, key, custom_kwargs) + else: + log.info(' --> Method ' + method_split[k] + ': skipped because TA file') + elif method_split[k] == 'timemed': + log.info(' --> Method ' + method_split[k] + ': ' + tail) + self.fix_bad_pixels_timemed(data, erro, pxdq_temp, timemed_kwargs) + elif method_split[k] == 'localmed': + log.info(' --> Method ' + method_split[k] + ': ' + tail) + self.fix_bad_pixels_localmed(data, erro, pxdq_temp, localmed_kwargs) + elif method_split[k] == 'medfilt': + log.info(' --> Method ' + method_split[k] + ': ' + tail) + self.fix_bad_pixels_medfilt(data, erro, pxdq_temp, medfilt_kwargs) + else: + log.info(' --> Unknown method ' + method_split[k] + ': skipped') + # if self.database.obs[key]['TELESCOP'][j] == 'JWST' and self.database.obs[key]['INSTRUME'][j] == 'NIRCAM': + # pxdq[(pxdq != 0) & np.logical_not(pxdq & 512 == 512) & (pxdq_temp == 0)] = 0 + # else: + # pxdq[(pxdq & 1 == 1) & np.logical_not(pxdq & 512 == 512) & (pxdq_temp == 0)] = 0 + + # update the pixel DQ bit flags for the output files. + # The pxdq variable here is effectively just the DO_NOT_USE flag, discarding other bits. + # We want to make a new dq which retains the other bits as much as possible. + # first, retain all the other bits (bits greater than 1), then add in the new/cleaned DO_NOT_USE bit + import jwst.datamodels + do_not_use = jwst.datamodels.dqflags.pixel['DO_NOT_USE'] + new_dq = np.bitwise_and(pxdq.copy(), np.invert(do_not_use)) # retain all other bits except the do_not_use bit + new_dq = np.bitwise_or(new_dq, pxdq_temp) # add in the do_not_use bit from the cleaned version + new_dq = new_dq.astype(np.uint32) # ensure correct output type for saving + # (the bitwise steps otherwise return np.int64 which isn't FITS compatible) + + # Write FITS file and PSF mask. + fitsfile = ut.write_obs(fitsfile, output_dir, data, erro, new_dq, head_pri, head_sci, is2d, imshifts, maskoffs) + maskfile = ut.write_msk(maskfile, mask, fitsfile) + + # Update spaceKLIP database. + self.database.update_obs(key, j, fitsfile, maskfile) + + pass + + def clean_bad_pixels(self, method='timemed+medfilt', timemed_kwargs={}, localmed_kwargs={}, From aa46abc4acb40856d5adabc56efe4d12ab6326e2 Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Wed, 7 Aug 2024 13:08:46 -0400 Subject: [PATCH 10/11] Add toggle when finding_bad_pixels to overwrite the existing dq array or start a new clean one. --- spaceKLIP/imagetools.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/spaceKLIP/imagetools.py b/spaceKLIP/imagetools.py index 674f0b8f..9d6116d4 100644 --- a/spaceKLIP/imagetools.py +++ b/spaceKLIP/imagetools.py @@ -769,6 +769,7 @@ def subtract_background(self, def find_bad_pixels(self, method='dqarr', + overwrite_dq=True, dqarr_kwargs={}, sigclip_kwargs={}, custom_kwargs={}, @@ -794,6 +795,10 @@ def find_bad_pixels(self, - custom: use a custom bad pixel map The default is 'dqarr'. + overwrite_dq : bool, optional + Toggle to start a new empty DQ array, or built upon the existing array. + + The default is True dqarr_kwargs : dict, optional Keyword arguments for the 'dqarr' identification method. Available keywords are: @@ -854,8 +859,11 @@ def find_bad_pixels(self, maskfile = self.database.obs[key]['MASKFILE'][j] mask = ut.read_msk(maskfile) - # Make copy of DQ array filled with zeros, i.e. all good pixels - pxdq_temp = np.zeros_like(pxdq) + if overwrite_dq: + # Make copy of DQ array filled with zeros, i.e. all good pixels + pxdq_temp = np.zeros_like(pxdq) + else: + pxdq_temp = pxdq # Skip file types that are not in the list of types. if self.database.obs[key]['TYPE'][j] in types: From 14cabec6eb04d0e8f8e3eec7c313006cf89cfe7a Mon Sep 17 00:00:00 2001 From: AarynnCarter <23636747+AarynnCarter@users.noreply.github.com> Date: Wed, 7 Aug 2024 13:39:29 -0400 Subject: [PATCH 11/11] Fixes for Jens' comments --- spaceKLIP/imagetools.py | 41 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 36 insertions(+), 5 deletions(-) diff --git a/spaceKLIP/imagetools.py b/spaceKLIP/imagetools.py index 9d6116d4..45755ff3 100644 --- a/spaceKLIP/imagetools.py +++ b/spaceKLIP/imagetools.py @@ -863,7 +863,7 @@ def find_bad_pixels(self, # Make copy of DQ array filled with zeros, i.e. all good pixels pxdq_temp = np.zeros_like(pxdq) else: - pxdq_temp = pxdq + pxdq_temp = pxdq.copy() # Skip file types that are not in the list of types. if self.database.obs[key]['TYPE'][j] in types: @@ -1072,7 +1072,7 @@ def fix_bad_pixels(self, pass def clean_bad_pixels(self, - method='timemed+medfilt', + method='timemed+localmed+medfilt', timemed_kwargs={}, localmed_kwargs={}, medfilt_kwargs={}, @@ -1093,14 +1093,30 @@ def clean_bad_pixels(self, - timemed: replace pixels which are only bad in some frames with their median value from the good frames. + - localmed: replace bad pixels with the median value of their + surrounding good pixels. + - medfilt: replace bad pixels with an image plane median filter. - The default is 'timemed+medfilt'. + - interp2d: replace bad pixels with an interpolation of neighbouring pixels. + + The default is 'timemed+localmed+medfilt'. timemed_kwargs : dict, optional Keyword arguments for the 'timemed' method. Available keywords are: - n/a + The default is {}. + localmed_kwargs: dict, optional + Keyword arguments for the 'localmed' method. Available keywords are: + + - shift_x : list of int, optional + Pixels in x-direction from which the median shall be computed. + The default is [-1, 0, 1]. + - shift_y : list of int, optional + Pixels in y-direction from which the median shall be computed. + The default is [-1, 0, 1]. + The default is {}. medfilt_kwargs : dict, optional Keyword arguments for the 'medfilt' method. Available keywords are: @@ -1109,6 +1125,14 @@ def clean_bad_pixels(self, Kernel size of the median filter to be used. The default is 4. The default is {}. + interp2d_kwargs: dict, optional + Keyword arguments for the 'interp2d' method. Available keywords are: + + - size : int, optional + Kernel size of the median filter to be used. The default is 4. + + The default is {}. + types : list of str, optional List of data types for which bad pixels shall be identified and fixed. The default is ['SCI', 'SCI_TA', 'SCI_BG', 'REF', 'REF_TA', @@ -1155,6 +1179,9 @@ def clean_bad_pixels(self, # Make copy of DQ array pxdq_temp = pxdq.copy() + + # Don't want to clean anything that isn't bad or is a non-science pixel + pxdq_temp = (np.isnan(data) | (pxdq_temp & 1 == 1)) & np.logical_not(pxdq_temp & 512 == 512) # Skip file types that are not in the list of types. if self.database.obs[key]['TYPE'][j] in types: @@ -1163,8 +1190,12 @@ def clean_bad_pixels(self, spatial = ['localmed', 'medfilt', 'interp2d'] # If localmed and medfilt in cleaning, can't run both if len(set(method_split) & set(spatial)) > 1: - log.info(' --> The localmed/medfilt/interp2d methods are redundant,') - log.info(' only the first method you listed will be affect the data.') + log.info(' --> WARNING: Multiple spatial cleaning routines detected!') + log.info(' --> The localmed/medfilt/interp2d methods clean data in a similar manner!') + log.info(' --> medfilt and interp2d are redundant') + log.info(' --> only the first method listed will affect the data') + log.info(' --> localmed is partially redundant with other methods') + log.info(' --> if run first, large clusters of bad pixels may not be fully cleaned.') # Loop over methods for k in range(len(method_split)):