From 7d2684425417dc89c2d289944cfe2678d8ce0999 Mon Sep 17 00:00:00 2001 From: Gabriel Brammer Date: Thu, 14 Sep 2023 16:25:23 +0200 Subject: [PATCH] options to precompute a contamination model before median filtering --- grizli/multifit.py | 55 ++++++++++++++++++++++++++-------- grizli/pipeline/auto_script.py | 18 +++++++++-- 2 files changed, 58 insertions(+), 15 deletions(-) diff --git a/grizli/multifit.py b/grizli/multifit.py index 770fbc8c..5344f3cd 100644 --- a/grizli/multifit.py +++ b/grizli/multifit.py @@ -1186,7 +1186,7 @@ def subtract_sep_background(self, mask_threshold=3, mask_iter=3, bw=256, bh=256, flt.grism.header['BKGBH'] = bh, 'sep background bh' - def subtract_median_filter(self, filter_size=71, filter_central=10, revert=True, filter_footprint=None, subtract_model=False, second_pass_filtering=False, box_filter_sn=3, box_filter_width=3): + def subtract_median_filter(self, filter_size=71, filter_central=10, revert=True, filter_footprint=None, subtract_model=False, second_pass_filtering=False, box_filter_sn=3, box_filter_width=3, put_model_in_median=False, verbose=True, mask_sn_threshold=None, mask_sn_dilate_iters=5): """ Remove a median filter calculated along the dispersion axis """ @@ -1210,16 +1210,13 @@ def subtract_median_filter(self, filter_size=71, filter_central=10, revert=True, msg += f' filter_size={filter_size} filter_central={filter_central}' msg += f' [{_filter_name}]' - utils.log_comment(utils.LOGFILE, msg, verbose=True) + utils.log_comment(utils.LOGFILE, msg, verbose=verbose) sci_i = flt.grism.data['SCI']*1 if revert & ('MED' in flt.grism.data): sci_i += flt.grism.data['MED'] - - if subtract_model: - sci_i -= flt.model - + err_i = flt.grism.data['ERR'] dq_i = flt.grism.data['DQ'] @@ -1229,6 +1226,12 @@ def subtract_median_filter(self, filter_size=71, filter_central=10, revert=True, ivar = 1/err_i**2 ivar[~ok] = 0 + if subtract_model: + msg = 'subtract_median_filter: subtract model before median' + utils.log_comment(utils.LOGFILE, msg, verbose=verbose) + + sci_i -= flt.model + #filter_sci = np.zeros_like(sci_i) #sh = sci_i.shape @@ -1245,16 +1248,18 @@ def subtract_median_filter(self, filter_size=71, filter_central=10, revert=True, if second_pass_filtering: if nbutils is None: # need numba installed - msg = 'subtract_median_filter: `numba` not found, skip second filter pass.' - utils.log_comment(utils.LOGFILE, msg, verbose=True) + msg = 'subtract_median_filter: `numba` not found, ' + msg += 'skip second filter pass.' + utils.log_comment(utils.LOGFILE, msg, verbose=verbose) else: # run filter again, but mask pixels that show significant residuals (e.g. strong emission lines) msg = f'subtract_median_filter: rerun filtering masking ' msg += f' S/N>{box_filter_sn} pixels in residual' - utils.log_comment(utils.LOGFILE, msg, verbose=True) + utils.log_comment(utils.LOGFILE, msg, verbose=verbose) # first do some binning/box filtering to identify significantly detected lines in individual exposures - box_filter_footprint = np.ones((box_filter_width, box_filter_width), + box_filter_footprint = np.ones((box_filter_width, + box_filter_width), dtype=int) box_filter_clean = nd.generic_filter(sci_i-filter_sci, nbutils.nansum, @@ -1272,9 +1277,35 @@ def subtract_median_filter(self, filter_size=71, filter_central=10, revert=True, footprint=filter_footprint) filter_sci[~np.isfinite(filter_sci)] = 0 + + if mask_sn_threshold: + flt.grism.header['MEDSNTH'] = (mask_sn_threshold, + 'Median mask threshold') + flt.grism.header['MEDSNIT'] = (mask_sn_dilate_iters, + 'Median mask threshold dilations') + + _msk = np.abs(filter_sci) > mask_sn_threshold*err_i + _msk = nd.binary_dilation(_msk, iterations=mask_sn_dilate_iters) + filter_sci[~_msk] = 0 + + if put_model_in_median & subtract_model: + msg = 'subtract_median_filter: put model in median and reset ' + msg += '`object_dispersers`' + utils.log_comment(utils.LOGFILE, msg, verbose=verbose) + + flt.grism.data['MED'] = (filter_sci + flt.model)*ok - flt.grism.data['SCI'] -= filter_sci*ok - flt.grism.data['MED'] = filter_sci*ok + # Reset model and object_dispersers + flt.object_dispersers = OrderedDict() + flt.model *= 0 + else: + flt.grism.data['MED'] = filter_sci*ok + + flt.grism.data['SCI'] -= flt.grism.data['MED'] + + flt.grism.header['MEDSMOD'] = subtract_model, 'Model subtracted first' + flt.grism.header['MEDWMOD'] = (put_model_in_median, + 'Model included in MED') flt.grism.header['MEDSIZE'] = filter_size, 'Median filter size' flt.grism.header['MEDCLIP'] = (filter_central, 'Masked central pixels of the filter') diff --git a/grizli/pipeline/auto_script.py b/grizli/pipeline/auto_script.py index 988a3490..2acf3180 100644 --- a/grizli/pipeline/auto_script.py +++ b/grizli/pipeline/auto_script.py @@ -3201,7 +3201,7 @@ def load_GroupFLT(field_root='j142724+334246', PREP_PATH='../Prep', force_ref=No return [grp] -def grism_prep(field_root='j142724+334246', PREP_PATH='../Prep', EXTRACT_PATH='../Extractions', ds9=None, refine_niter=3, gris_ref_filters=GRIS_REF_FILTERS, force_ref=None, files=None, split_by_grism=True, refine_poly_order=1, refine_fcontam=0.5, cpu_count=0, mask_mosaic_edges=False, prelim_mag_limit=25, refine_mag_limits=[18, 24], init_coeffs=[1.1, -0.5], grisms_to_process=None, pad=(64, 256), model_kwargs={'compute_size': True}, sep_background_kwargs=None, subtract_median_filter=False, median_filter_size=71, median_filter_central=10, second_pass_filtering=False, box_filter_sn=3, box_filter_width=3, use_jwst_crds=False): +def grism_prep(field_root='j142724+334246', PREP_PATH='../Prep', EXTRACT_PATH='../Extractions', ds9=None, refine_niter=3, gris_ref_filters=GRIS_REF_FILTERS, force_ref=None, files=None, split_by_grism=True, refine_poly_order=1, refine_fcontam=0.5, cpu_count=0, mask_mosaic_edges=False, prelim_mag_limit=25, refine_mag_limits=[18, 24], init_coeffs=[1.1, -0.5], grisms_to_process=None, pad=(64, 256), model_kwargs={'compute_size': True}, sep_background_kwargs=None, subtract_median_filter=False, median_filter_size=71, median_filter_central=10, second_pass_filtering=False, box_filter_sn=3, box_filter_width=3, median_mask_sn_threshold=None, median_mask_dilate=8, prelim_model_for_median=False, use_jwst_crds=False): """ Contamination model for grism exposures """ @@ -3236,14 +3236,26 @@ def grism_prep(field_root='j142724+334246', PREP_PATH='../Prep', EXTRACT_PATH='. if subtract_median_filter: + ################ + # Compute preliminary model before median? + if prelim_model_for_median: + grp.compute_full_model(fit_info=None, verbose=True, store=False, + mag_limit=prelim_mag_limit, + coeffs=init_coeffs, + cpu_count=cpu_count, + model_kwargs=model_kwargs) + # Subtract a median along the dispersion direction - # and don't do any contamination modeling refine_niter = 0 grp.subtract_median_filter(filter_size=median_filter_size, filter_central=median_filter_central, second_pass_filtering=second_pass_filtering, box_filter_sn=box_filter_sn, - box_filter_width=box_filter_width) + box_filter_width=box_filter_width, + mask_sn_threshold=median_mask_sn_threshold, + mask_sn_dilate_iters=median_mask_dilate, + subtract_model=prelim_model_for_median, + put_model_in_median=prelim_model_for_median) else: ################