diff --git a/jwst/extract_1d/extract.py b/jwst/extract_1d/extract.py index 0bd633ad88..57bca3620c 100644 --- a/jwst/extract_1d/extract.py +++ b/jwst/extract_1d/extract.py @@ -8,6 +8,7 @@ from typing import Union, Tuple, NamedTuple, List from astropy.modeling import polynomial from astropy.io import fits +from scipy.optimize import curve_fit from gwcs import WCS from stdatamodels import DataModel from stdatamodels.properties import ObjectNode @@ -88,6 +89,10 @@ EXACT = "exact match" +def gaussian(x, a, b, c): + return c / a * np.exp(-1. * (x - b)**2 / (2. * a * a)) + + class Aperture(NamedTuple): # When python 3.6 is no longer supported, consider converting to DataClass xstart: Union[int, float] xstop: Union[int, float] @@ -220,7 +225,8 @@ def get_extract_parameters( smoothing_length: Union[int, None], bkg_fit: str, bkg_order: Union[int, None], - use_source_posn: Union[bool, None] + use_source_posn: Union[bool, None], + slit_autocen: Union[bool, None] ) -> dict: """Get extract1d reference file values. @@ -283,6 +289,10 @@ def get_extract_parameters( If None, the value specified in `ref_dict` will be used, or it will be set to True if not found in `ref_dict`. + slit_autocen : bool + Switch to turn on auto-centering for point source spectral extraction + in slit-like modes. Default is True. + Returns ------- extract_params : dict @@ -310,6 +320,7 @@ def get_extract_parameters( extract_params['bkg_fit'] = None # because no background sub. extract_params['bkg_order'] = 0 # because no background sub. extract_params['subtract_background'] = False + extract_params['slit_autocen'] = slit_autocen if use_source_posn is None: extract_params['use_source_posn'] = False @@ -386,7 +397,7 @@ def get_extract_parameters( if use_source_posn is None: # no value set on command line if use_source_posn_aper is None: # no value set in ref file # Use a suitable default - if meta.exposure.type in ['MIR_LRS-FIXEDSLIT', 'MIR_MRS', 'NRS_FIXEDSLIT', 'NRS_IFU', 'NRS_MSASPEC']: + if meta.exposure.type in ['MIR_MRS', 'NRS_IFU']: use_source_posn = True log.info(f"Turning on source position correction for exp_type = {meta.exposure.type}") else: @@ -394,6 +405,7 @@ def get_extract_parameters( else: # use the value from the ref file use_source_posn = use_source_posn_aper extract_params['use_source_posn'] = use_source_posn + extract_params['slit_autocen'] = slit_autocen extract_params['extract_width'] = aper.get('extract_width') extract_params['position_correction'] = 0 # default value @@ -1056,7 +1068,8 @@ def __init__( subtract_background: Union[bool, None] = None, use_source_posn: Union[bool, None] = None, match: Union[str, None] = None, - ref_file_type: Union[str, None] = None + ref_file_type: Union[str, None] = None, + slit_autocen: Union[bool, None] = None ): """ Parameters @@ -1198,6 +1211,7 @@ def __init__( self.use_source_posn = use_source_posn self.position_correction = position_correction self.subtract_background = subtract_background + self.slit_autocen = slit_autocen self.wcs = None # initial value @@ -1417,7 +1431,7 @@ def locn_from_wcs( # Width (height) in the cross-dispersion direction, from the start of the 2-D cutout (or of the full image) # to the upper limit of the bounding box. # This may be smaller than the full width of the image, but it's all we need to consider. - xd_width = int(round(bb[1][1])) # must be an int + xd_width = int(round(bb[1][1])+0.0001) # must be an int middle = int((bb[0][0] + bb[0][1]) / 2.) # Middle of the bounding_box in the dispersion direction. x = np.empty(xd_width, dtype=np.float64) x[:] = float(middle) @@ -2524,6 +2538,7 @@ def run_extract1d( subtract_background: Union[bool, None], use_source_posn: Union[bool, None], center_xy: Union[float, None], + slit_autocen: Union[bool, None], ifu_autocen: Union[bool, None], ifu_rfcorr: Union[bool, None], was_source_model: bool = False, @@ -2578,6 +2593,10 @@ def run_extract1d( done by the step. Two values, in x,y order, are used for extraction from IFU cubes. Default is None. + slit_autocen : bool + Switch to turn on auto-centering for point source spectral extraction + in slit-like modes. Default is True. + ifu_autocen : bool Switch to turn on auto-centering for point source spectral extraction in IFU mode. Default is False. @@ -2632,6 +2651,7 @@ def run_extract1d( subtract_background, use_source_posn, center_xy, + slit_autocen, ifu_autocen, ifu_rfcorr, was_source_model, @@ -2694,6 +2714,7 @@ def do_extract1d( subtract_background: Union[int, None] = None, use_source_posn: Union[bool, None] = None, center_xy: Union[int, None] = None, + slit_autocen: Union[bool, None] = None, ifu_autocen: Union[bool, None] = None, ifu_rfcorr: Union[bool, None] = None, was_source_model: bool = False @@ -2758,6 +2779,10 @@ def do_extract1d( Two values, in x,y order, are used for extraction from IFU cubes. Default is None. + slit_autocen : bool + Switch to turn on auto-centering for point source spectral extraction + in slit-like modes. Default is True. + ifu_autocen : bool Switch to turn on auto-centering for point source spectral extraction in IFU mode. Default is False. @@ -2869,7 +2894,7 @@ def do_extract1d( smoothing_length, bkg_fit, bkg_order, use_source_posn, prev_offset, exp_type, subtract_background, input_model, output_model, apcorr_ref_model, log_increment, - is_multiple_slits + is_multiple_slits, slit_autocen ) except ContinueError: continue @@ -2918,7 +2943,7 @@ def do_extract1d( smoothing_length, bkg_fit, bkg_order, use_source_posn, prev_offset, exp_type, subtract_background, input_model, output_model, apcorr_ref_model, log_increment, - is_multiple_slits + is_multiple_slits, slit_autocen ) except ContinueError: continue @@ -2959,7 +2984,7 @@ def do_extract1d( smoothing_length, bkg_fit, bkg_order, use_source_posn, prev_offset, exp_type, subtract_background, input_model, output_model, apcorr_ref_model, log_increment, - is_multiple_slits + is_multiple_slits, slit_autocen ) except ContinueError: continue @@ -3451,6 +3476,54 @@ def extract_one_slit( extract_model.log_extraction_parameters() extract_model.assign_polynomial_limits() + if extract_model.slit_autocen is True: + good_pixels = np.where(data == data) + y_range = [np.min(good_pixels[0]), np.max(good_pixels[0]) + 1] + x_range = [np.min(good_pixels[1]), np.max(good_pixels[1]) + 1] + extract_model.xstart = x_range[0] + extract_model.xstop = x_range[1] + extract_model.ystart = y_range[0] + extract_model.ystop = y_range[1] + + absdata = data.copy() + absdata = absdata[y_range[0]: y_range[1], x_range[0]: x_range[1]] + + absdata[absdata < 0] = 0 + + normed = absdata / np.nanmax(absdata, axis=extract_params['dispaxis'] - 1, keepdims=True) + + y_fit = np.nanmedian(normed, axis=int(-1 * extract_params['dispaxis'])) + x_fit = np.arange(len(y_fit)) + y_peak = np.argmax(y_fit) + # Provide bounds to fit to avoid wasting time in region where a fit is undesirable + try: + fit_params, _ = curve_fit(gaussian, x_fit, y_fit, bounds=(np.array([0, y_peak - 3, 0.5]), + np.array([8, y_peak + 3, 10])) + ) + fit_width, fit_center, _ = fit_params + n_sig = 3. + if extract_model.dispaxis == HORIZONTAL: + extract_model.ystart = max(0, fit_center - (n_sig * fit_width)) + y_range[0] + extract_model.ystop = min(data.shape[-2] - 1, fit_center + (n_sig * fit_width)) + y_range[0] + extract_model.p_src = [[create_poly([extract_model.ystart]), + create_poly([extract_model.ystop])]] + log.debug(f"Slit auto-centering fit_center: {fit_center + y_range[0]}" + f"and fit_width: {fit_width}") + + else: + extract_model.xstart = max(0, fit_center - (n_sig * fit_width)) + x_range[0] + extract_model.xstop = min(data.shape[-1] - 1, fit_center + (n_sig * fit_width)) + x_range[0] + extract_model.p_src = [[create_poly([extract_model.xstart]), + create_poly([extract_model.xstop])]] + + log.debug(f"Slit auto-centering fit_center: {fit_center + x_range[0]}" + f"and fit_width: {fit_width}") + + log.info("Overriding cross-dispersion bounds with slit auto-centering.") + except (ValueError, RuntimeError) as e: + log.warning("Auto-centering failed; reverting to default behavior.") + log.debug(f"Auto-centering error: {e}") + # Log the extraction limits being used if integ < 1: if extract_model.src_coeff is not None: @@ -3599,7 +3672,8 @@ def create_extraction(extract_ref_dict, output_model, apcorr_ref_model, log_increment, - is_multiple_slits): + is_multiple_slits, + slit_autocen): if slit is None: meta_source = input_model else: @@ -3648,7 +3722,8 @@ def create_extraction(extract_ref_dict, # Turn off use_source_posn if the source is not POINT if source_type != 'POINT': use_source_posn = False - log.info(f"Setting use_source_posn to False for source type {source_type}") + slit_autocen = False + log.info(f"Setting use_source_posn and slit auto-centering to False for source type {source_type}") # Turn off use_source_posn if working on non-primary NRS fixed slits if is_multiple_slits: @@ -3674,7 +3749,8 @@ def create_extraction(extract_ref_dict, smoothing_length, bkg_fit, bkg_order, - use_source_posn + use_source_posn, + slit_autocen ) if subtract_background is not None: diff --git a/jwst/extract_1d/extract_1d_step.py b/jwst/extract_1d/extract_1d_step.py index 22f9b0e980..46e0231368 100644 --- a/jwst/extract_1d/extract_1d_step.py +++ b/jwst/extract_1d/extract_1d_step.py @@ -147,6 +147,7 @@ class Extract1dStep(Step): use_source_posn = boolean(default=None) # use source coords to center extractions? center_xy = float_list(min=2, max=2, default=None) # IFU extraction x/y center apply_apcorr = boolean(default=True) # apply aperture corrections? + slit_autocen = boolean(default=True) # Auto source centering for slit-like point source data. ifu_autocen = boolean(default=False) # Auto source centering for IFU point source data. ifu_rfcorr = boolean(default=False) # Apply 1d residual fringe correction soss_atoca = boolean(default=True) # use ATOCA algorithm @@ -258,6 +259,7 @@ def process(self, input): self.subtract_background, self.use_source_posn, self.center_xy, + self.slit_autocen, self.ifu_autocen, self.ifu_rfcorr, was_source_model=was_source_model @@ -294,6 +296,7 @@ def process(self, input): self.subtract_background, self.use_source_posn, self.center_xy, + self.slit_autocen, self.ifu_autocen, self.ifu_rfcorr, was_source_model=was_source_model, @@ -333,6 +336,7 @@ def process(self, input): self.subtract_background, self.use_source_posn, self.center_xy, + self.slit_autocen, self.ifu_autocen, self.ifu_rfcorr, was_source_model=was_source_model, @@ -474,6 +478,7 @@ def process(self, input): self.subtract_background, self.use_source_posn, self.center_xy, + self.slit_autocen, self.ifu_autocen, self.ifu_rfcorr, was_source_model=False,