diff --git a/doc/calibrations/slit_tracing.rst b/doc/calibrations/slit_tracing.rst index 84bbdbd574..af2b1ba969 100644 --- a/doc/calibrations/slit_tracing.rst +++ b/doc/calibrations/slit_tracing.rst @@ -369,7 +369,7 @@ case for low-dispersion data, e.g. LRISb 300 grism spectra .. code-block:: ini [calibrations] - [[slits]] + [[slitedges]] smash_range = 0.5,1. Algorithm diff --git a/pypeit/alignframe.py b/pypeit/alignframe.py index 9a28491d21..2cfc6b28d8 100644 --- a/pypeit/alignframe.py +++ b/pypeit/alignframe.py @@ -181,8 +181,8 @@ def build_traces(self, show_peaks=False, debug=False): dict: self.align_dict """ # Generate slits - slitid_img_init = self.slits.slit_img(initial=True) - left, right, _ = self.slits.select_edges(initial=True) + slitid_img_init = self.slits.slit_img() + left, right, _ = self.slits.select_edges() align_prof = dict({}) # Go through the slits diff --git a/pypeit/calibrations.py b/pypeit/calibrations.py index e8875595df..0d55d99674 100644 --- a/pypeit/calibrations.py +++ b/pypeit/calibrations.py @@ -607,7 +607,7 @@ def get_scattlight(self): spatbin = parse.parse_binning(binning)[1] pad = self.par['scattlight_pad'] // spatbin - offslitmask = self.slits.slit_img(pad=pad, initial=True, flexure=None) == -1 + offslitmask = self.slits.slit_img(pad=pad, flexure=None) == -1 # Get starting parameters for the scattered light model x0, bounds = self.spectrograph.scattered_light_archive(binning, dispname) @@ -766,6 +766,7 @@ def get_flats(self): msgs.info('Creating slit-illumination flat calibration frame using files: ') for f in raw_illum_files: msgs.prindent(f'{Path(f).name}') + illum_flat = buildimage.buildimage_fromlist(self.spectrograph, self.det, self.par['illumflatframe'], raw_illum_files, dark=self.msdark, bias=self.msbias, scattlight=self.msscattlight, diff --git a/pypeit/coadd3d.py b/pypeit/coadd3d.py index 612d91562b..225747d76d 100644 --- a/pypeit/coadd3d.py +++ b/pypeit/coadd3d.py @@ -838,8 +838,8 @@ def add_grating_corr(self, flatfile, waveimg, slits, spat_flexure=None): msgs.info("Calculating relative sensitivity for grating correction") # Load the Flat file flatimages = flatfield.FlatImages.from_file(flatfile, chk_version=self.chk_version) - total_illum = flatimages.fit2illumflat(slits, finecorr=False, frametype='illum', initial=True, spat_flexure=spat_flexure) * \ - flatimages.fit2illumflat(slits, finecorr=True, frametype='illum', initial=True, spat_flexure=spat_flexure) + total_illum = flatimages.fit2illumflat(slits, finecorr=False, frametype='illum', spat_flexure=spat_flexure) * \ + flatimages.fit2illumflat(slits, finecorr=True, frametype='illum', spat_flexure=spat_flexure) flatframe = flatimages.pixelflat_raw / total_illum if flatimages.pixelflat_spec_illum is None: # Calculate the relative scale @@ -950,7 +950,7 @@ def get_alignments(self, spec2DObj, slits, spat_flexure=None): msgs.info("Using slit edges for astrometric transform") # If nothing better was provided, use the slit edges if alignments is None: - left, right, _ = slits.select_edges(initial=True, flexure=spat_flexure) + left, right, _ = slits.select_edges(flexure=spat_flexure) locations = [0.0, 1.0] traces = np.append(left[:, None, :], right[:, None, :], axis=1) else: @@ -1015,8 +1015,8 @@ def load(self): # Initialise the slit edges msgs.info("Constructing slit image") slits = spec2DObj.slits - slitid_img_init = slits.slit_img(pad=0, initial=True, flexure=spat_flexure) - slits_left, slits_right, _ = slits.select_edges(initial=True, flexure=spat_flexure) + slitid_img = slits.slit_img(pad=0, flexure=spat_flexure) + slits_left, slits_right, _ = slits.select_edges(flexure=spat_flexure) # The order of operations below proceeds as follows: # (1) Get science image @@ -1041,7 +1041,7 @@ def load(self): bpmmask = spec2DObj.bpmmask # Mask the edges of the spectrum where the sky model is bad - sky_is_good = datacube.make_good_skymask(slitid_img_init, spec2DObj.tilts) + sky_is_good = datacube.make_good_skymask(slitid_img, spec2DObj.tilts) # TODO :: Really need to write some detailed information in the docs about all of the various corrections that can optionally be applied @@ -1074,7 +1074,7 @@ def load(self): if self.mnmx_wv is None: self.mnmx_wv = np.zeros((len(self.spec2d), slits.nslits, 2)) for slit_idx, slit_spat in enumerate(slits.spat_id): - onslit_init = (slitid_img_init == slit_spat) + onslit_init = (slitid_img == slit_spat) self.mnmx_wv[ff, slit_idx, 0] = np.min(waveimg[onslit_init]) self.mnmx_wv[ff, slit_idx, 1] = np.max(waveimg[onslit_init]) @@ -1097,7 +1097,7 @@ def load(self): # Construct a good pixel mask # TODO: This should use the mask function to figure out which elements are masked. - onslit_gpm = (slitid_img_init > 0) & (bpmmask.mask == 0) & sky_is_good + onslit_gpm = (slitid_img > 0) & (bpmmask.mask == 0) & sky_is_good # Generate the alignment splines, and then retrieve images of the RA and Dec of every pixel, # and the number of spatial pixels in each slit @@ -1109,7 +1109,7 @@ def load(self): crval_wv = self.cubepar['wave_min'] if self.cubepar['wave_min'] is not None else wave0 cd_wv = self.cubepar['wave_delta'] if self.cubepar['wave_delta'] is not None else dwv self.all_wcs.append(self.spec.get_wcs(spec2DObj.head0, slits, detector.platescale, crval_wv, cd_wv)) - ra_img, dec_img, minmax = slits.get_radec_image(self.all_wcs[ff], alignSplines, spec2DObj.tilts, initial=True, flexure=spat_flexure) + ra_img, dec_img, minmax = slits.get_radec_image(self.all_wcs[ff], alignSplines, spec2DObj.tilts, flexure=spat_flexure) # Extract wavelength and delta wavelength arrays from the images wave_ext = waveimg[onslit_gpm] @@ -1198,7 +1198,7 @@ def load(self): wghts = self.weights[ff] * np.ones(sciImg.shape) # Get the slit image and then unset pixels in the slit image that are bad - slitid_img_gpm = slitid_img_init * onslit_gpm.astype(int) + slitid_img_gpm = slitid_img * onslit_gpm.astype(int) # If individual frames are to be output without aligning them, # there's no need to store information, just make the cubes now @@ -1209,7 +1209,7 @@ def load(self): else: outfile = datacube.get_output_filename(fil, self.cubepar['output_filename'], self.combine, ff + 1) # Get the coordinate bounds - slitlength = int(np.round(np.median(slits.get_slitlengths(initial=True, median=True)))) + slitlength = int(np.round(np.median(slits.get_slitlengths(median=True)))) numwav = int((np.max(waveimg) - wave0) / dwv) bins = self.spec.get_datacube_bins(slitlength, minmax, numwav) # Set the wavelength range of the white light image. diff --git a/pypeit/core/datacube.py b/pypeit/core/datacube.py index 66b3c419df..7688753ade 100644 --- a/pypeit/core/datacube.py +++ b/pypeit/core/datacube.py @@ -596,13 +596,16 @@ def get_whitelight_range(wavemin, wavemax, wl_range): return wlrng -def make_whitelight_fromcube(cube, wave=None, wavemin=None, wavemax=None): +def make_whitelight_fromcube(cube, bpmcube, wave=None, wavemin=None, wavemax=None): """ Generate a white light image using an input cube. Args: cube (`numpy.ndarray`_): 3D datacube (the final element contains the wavelength dimension) + bpmcube (`numpy.ndarray`_): + 3D bad pixel mask cube (the final element contains the wavelength dimension). + A value of 1 indicates a bad pixel. wave (`numpy.ndarray`_, optional): 1D wavelength array. Only required if wavemin or wavemax are not None. @@ -619,7 +622,6 @@ def make_whitelight_fromcube(cube, wave=None, wavemin=None, wavemax=None): A whitelight image of the input cube (of type `numpy.ndarray`_). """ # Make a wavelength cut, if requested - cutcube = cube.copy() if wavemin is not None or wavemax is not None: # Make some checks on the input if wave is None: @@ -635,10 +637,15 @@ def make_whitelight_fromcube(cube, wave=None, wavemin=None, wavemax=None): ww = np.where((wave >= wavemin) & (wave <= wavemax))[0] wmin, wmax = ww[0], ww[-1]+1 cutcube = cube[:, :, wmin:wmax] + # Cut the bad pixel mask and convert it to a good pixel mask + cutgpmcube = np.logical_not(bpmcube[:, :, wmin:wmax]) + else: + cutcube = cube.copy() + cutgpmcube = np.logical_not(bpmcube) # Now sum along the wavelength axis - nrmval = np.sum(cutcube != 0.0, axis=2) - nrmval[nrmval == 0.0] = 1.0 - wl_img = np.sum(cutcube, axis=2) / nrmval + nrmval = np.sum(cutgpmcube, axis=2) + nrmval[nrmval == 0] = 1.0 + wl_img = np.sum(cutcube*cutgpmcube, axis=2) / nrmval return wl_img @@ -1594,7 +1601,7 @@ def generate_cube_subpixel(output_wcs, bins, sciImg, ivarImg, waveImg, slitid_im whitelight_range[0], whitelight_range[1])) # Get the output filename for the white light image out_whitelight = get_output_whitelight_filename(outfile) - whitelight_img = make_whitelight_fromcube(flxcube, wave=wave, wavemin=whitelight_range[0], wavemax=whitelight_range[1]) + whitelight_img = make_whitelight_fromcube(flxcube, bpmcube, wave=wave, wavemin=whitelight_range[0], wavemax=whitelight_range[1]) msgs.info("Saving white light image as: {0:s}".format(out_whitelight)) img_hdu = fits.PrimaryHDU(whitelight_img.T, header=whitelight_wcs.to_header()) img_hdu.writeto(out_whitelight, overwrite=overwrite) @@ -1698,8 +1705,7 @@ def subpixellate(output_wcs, bins, sciImg, ivarImg, waveImg, slitid_img_gpm, wgh Returns: :obj:`tuple`: Three or four `numpy.ndarray`_ objects containing (1) the datacube generated from the subpixellated inputs, (2) the corresponding - variance cube, (3) the corresponding bad pixel mask cube, and (4) the - residual cube. The latter is only returned if debug is True. + variance cube, and (3) the corresponding bad pixel mask cube. """ # Check the inputs for combinations of lists or not _sciImg, _ivarImg, _waveImg, _gpmImg, _wghtImg, _all_wcs, _tilts, _slits, _astrom_trans, _all_dar, _ra_offset, _dec_offset = \ diff --git a/pypeit/core/flat.py b/pypeit/core/flat.py index 1c104d1d96..46495f56f2 100644 --- a/pypeit/core/flat.py +++ b/pypeit/core/flat.py @@ -16,44 +16,8 @@ from IPython import embed from pypeit import msgs -from pypeit.core import parse -from pypeit.core import pixels -from pypeit.core import tracewave from pypeit.core import coadd from pypeit import utils -from pypeit.core import pydl - -# TODO: Put this in utils -def linear_interpolate(x1, y1, x2, y2, x): - r""" - Interplate or extrapolate between two points. - - Given a line defined two points, :math:`(x_1,y_1)` and - :math:`(x_2,y_2)`, return the :math:`y` value of a new point on - the line at coordinate :math:`x`. - - This function is meant for speed. No type checking is performed and - the only check is that the two provided ordinate coordinates are not - numerically identical. By definition, the function will extrapolate - without any warning. - - Args: - x1 (:obj:`float`): - First abscissa position - y1 (:obj:`float`): - First ordinate position - x2 (:obj:`float`): - Second abscissa position - y3 (:obj:`float`): - Second ordinate position - x (:obj:`float`): - Abcissa for new value - - Returns: - :obj:`float`: Interpolated/extrapolated value of ordinate at - :math:`x`. - """ - return y1 if np.isclose(x1,x2) else y1 + (x-x1)*(y2-y1)/(x2-x1) # TODO: Make this function more general and put it in utils. @@ -505,10 +469,120 @@ def poly_map(rawimg, rawivar, waveimg, slitmask, slitmask_trim, modelimg, deg=3, return modelmap, relscale +def tweak_slit_edges_gradient(left, right, spat_coo, norm_flat, maxfrac=0.1, debug=False): + r""" Adjust slit edges based on the gradient of the normalized + flat-field illumination profile. + + Args: + left (`numpy.ndarray`_): + Array with the left slit edge for a single slit. Shape is + :math:`(N_{\rm spec},)`. + right (`numpy.ndarray`_): + Array with the right slit edge for a single slit. Shape + is :math:`(N_{\rm spec},)`. + spat_coo (`numpy.ndarray`_): + Spatial pixel coordinates in fractions of the slit width + at each spectral row for the provided normalized flat + data. Coordinates are relative to the left edge (with the + left edge at 0.). Shape is :math:`(N_{\rm flat},)`. + Function assumes the coordinate array is sorted. + norm_flat (`numpy.ndarray`_) + Normalized flat data that provide the slit illumination + profile. Shape is :math:`(N_{\rm flat},)`. + maxfrac (:obj:`float`, optional): + The maximum fraction of the slit width that the slit edge + can be adjusted by this algorithm. If ``maxfrac = 0.1``, + this means the maximum change in the slit width (either + narrowing or broadening) is 20% (i.e., 10% for either + edge). + debug (:obj:`bool`, optional): + If True, the function will output plots to test if the + fitting is working correctly. + + Returns: + tuple: Returns six objects: + + - The threshold used to set the left edge + - The fraction of the slit that the left edge is shifted to + the right + - The adjusted left edge + - The threshold used to set the right edge + - The fraction of the slit that the right edge is shifted to + the left + - The adjusted right edge + """ + # Check input + nspec = len(left) + if len(right) != nspec: + msgs.error('Input left and right traces must have the same length!') + + # Median slit width + slitwidth = np.median(right - left) + + # Calculate the gradient of the normalized flat profile + grad_norm_flat = np.gradient(norm_flat) + # Smooth with a Gaussian kernel + # The standard deviation of the kernel is set to be one detector pixel. Since the norm_flat array is oversampled, + # we need to set the kernel width (sig_res) to be the oversampling factor. + sig_res = norm_flat.size / slitwidth + # The scipy.ndimage module is faster than the astropy convolution module + grad_norm_flat_smooth = scipy.ndimage.gaussian_filter1d(grad_norm_flat, sig_res, mode='nearest') + + # Find the location of the minimum/maximum gradient - this is the amount of shift required + left_shift = spat_coo[np.argmax(grad_norm_flat_smooth)] + right_shift = spat_coo[np.argmin(grad_norm_flat_smooth)]-1.0 + + # Check if the shift is within the allowed range + if np.abs(left_shift) > maxfrac: + msgs.warn('Left slit edge shift of {0:.1f}% exceeds the maximum allowed of {1:.1f}%'.format( + 100*left_shift, 100*maxfrac) + msgs.newline() + + 'The left edge will not be tweaked.') + left_shift = 0.0 + else: + msgs.info('Tweaking left slit boundary by {0:.1f}%'.format(100 * left_shift) + + ' ({0:.2f} pixels)'.format(left_shift * slitwidth)) + if np.abs(right_shift) > maxfrac: + msgs.warn('Right slit edge shift of {0:.1f}% exceeds the maximum allowed of {1:.1f}%'.format( + 100*right_shift, 100*maxfrac) + msgs.newline() + + 'The right edge will not be tweaked.') + right_shift = 0.0 + else: + msgs.info('Tweaking right slit boundary by {0:.1f}%'.format(100 * right_shift) + + ' ({0:.2f} pixels)'.format(right_shift * slitwidth)) + + # Calculate the tweak for the left edge + new_left = left + left_shift * slitwidth + new_right = right + right_shift * slitwidth + + # Calculate the value of the threshold at the new slit edges + left_thresh = np.interp(left_shift, spat_coo, norm_flat) + right_thresh = np.interp(1+right_shift, spat_coo, norm_flat) + + if debug: + plt.subplot(211) + plt.plot(spat_coo, norm_flat, 'k-') + plt.axvline(0.0, color='b', linestyle='-', label='initial') + plt.axvline(1.0, color='b', linestyle='-') + plt.axvline(left_shift, color='g', linestyle='-', label='tweak (gradient)') + plt.axvline(1+right_shift, color='g', linestyle='-') + plt.axhline(left_thresh, xmax=0.5, color='lightgreen', linewidth=3.0, zorder=10) + plt.axhline(right_thresh, xmin=0.5, color='lightgreen', linewidth=3.0, zorder=10) + plt.legend() + plt.subplot(212) + plt.plot(spat_coo, grad_norm_flat, 'k-') + plt.plot(spat_coo, grad_norm_flat_smooth, 'm-') + plt.axvline(0.0, color='b', linestyle='-') + plt.axvline(1.0, color='b', linestyle='-') + plt.axvline(left_shift, color='g', linestyle='-') + plt.axvline(1+right_shift, color='g', linestyle='-') + plt.show() + return left_thresh, left_shift, new_left, right_thresh, right_shift, new_right + + # TODO: See pypeit/deprecated/flat.py for the previous version. We need # to continue to vet this algorithm to make sure there are no # unforeseen corner cases that cause errors. -def tweak_slit_edges(left, right, spat_coo, norm_flat, thresh=0.93, maxfrac=0.1, debug=False): +def tweak_slit_edges_threshold(left, right, spat_coo, norm_flat, thresh=0.93, maxfrac=0.1, debug=False): r""" Adjust slit edges based on the normalized slit illumination profile. @@ -614,10 +688,10 @@ def tweak_slit_edges(left, right, spat_coo, norm_flat, thresh=0.93, maxfrac=0.1, 100*maxfrac)) left_shift = maxfrac else: - left_shift = linear_interpolate(norm_flat[i], spat_coo[i], norm_flat[i+1], - spat_coo[i+1], left_thresh) + left_shift = utils.linear_interpolate(norm_flat[i], spat_coo[i], norm_flat[i+1], + spat_coo[i+1], left_thresh) msgs.info('Tweaking left slit boundary by {0:.1f}%'.format(100*left_shift) + - ' % ({0:.2f} pixels)'.format(left_shift*slitwidth)) + ' ({0:.2f} pixels)'.format(left_shift*slitwidth)) new_left += left_shift * slitwidth # ------------------------------------------------------------------ @@ -668,15 +742,15 @@ def tweak_slit_edges(left, right, spat_coo, norm_flat, thresh=0.93, maxfrac=0.1, 100*maxfrac)) right_shift = maxfrac else: - right_shift = 1-linear_interpolate(norm_flat[i-1], spat_coo[i-1], norm_flat[i], - spat_coo[i], right_thresh) + right_shift = 1-utils.linear_interpolate(norm_flat[i-1], spat_coo[i-1], norm_flat[i], + spat_coo[i], right_thresh) msgs.info('Tweaking right slit boundary by {0:.1f}%'.format(100*right_shift) + - ' % ({0:.2f} pixels)'.format(right_shift*slitwidth)) + ' ({0:.2f} pixels)'.format(right_shift*slitwidth)) new_right -= right_shift * slitwidth return left_thresh, left_shift, new_left, right_thresh, right_shift, new_right -#def flatfield(sciframe, flatframe, bpm=None, illum_flat=None, snframe=None, varframe=None): + def flatfield(sciframe, flatframe, varframe=None): r""" Field flatten the input image. diff --git a/pypeit/core/procimg.py b/pypeit/core/procimg.py index 6661c3d257..c68be8d379 100644 --- a/pypeit/core/procimg.py +++ b/pypeit/core/procimg.py @@ -813,6 +813,7 @@ def subtract_pattern(rawframe, datasec_img, oscansec_img, frequency=None, axis=1 frequency = np.mean(frq) # Perform the overscan subtraction for each amplifier + full_model = np.zeros_like(frame_orig) # Store the model pattern for all amplifiers in this array for aa, amp in enumerate(amps): # Get the frequency to use for this amplifier if isinstance(frequency, list): @@ -823,9 +824,9 @@ def subtract_pattern(rawframe, datasec_img, oscansec_img, frequency=None, axis=1 use_fr = frequency # Extract overscan - overscan, os_slice = rect_slice_with_mask(frame_orig, tmp_oscan, amp) + overscan, os_slice = rect_slice_with_mask(frame_orig.copy(), tmp_oscan, amp) # Extract overscan+data - oscandata, osd_slice = rect_slice_with_mask(frame_orig, tmp_oscan+tmp_data, amp) + oscandata, osd_slice = rect_slice_with_mask(frame_orig.copy(), tmp_oscan+tmp_data, amp) # Subtract the DC offset overscan -= np.median(overscan, axis=1)[:, np.newaxis] @@ -854,7 +855,7 @@ def subtract_pattern(rawframe, datasec_img, oscansec_img, frequency=None, axis=1 tmpamp = np.fft.rfft(overscan, axis=1) idx = (np.arange(overscan.shape[0]), np.argmax(np.abs(tmpamp), axis=1)) # Convert result to amplitude and phase - amps = (np.abs(tmpamp))[idx] * (2.0 / overscan.shape[1]) + ampls = (np.abs(tmpamp))[idx] * (2.0 / overscan.shape[1]) # STEP 2 - Using the model frequency, calculate how amplitude depends on pixel row (usually constant) # Use the above to as initial guess parameters for a chi-squared minimisation of the amplitudes @@ -877,7 +878,7 @@ def subtract_pattern(rawframe, datasec_img, oscansec_img, frequency=None, axis=1 try: # Now fit it popt, pcov = scipy.optimize.curve_fit( - cosfunc, cent[wgd], hist[wgd], p0=[amps[ii], 0.0], + cosfunc, cent[wgd], hist[wgd], p0=[ampls[ii], 0.0], bounds=([0, -np.inf],[np.inf, np.inf]) ) except ValueError: @@ -920,21 +921,15 @@ def subtract_pattern(rawframe, datasec_img, oscansec_img, frequency=None, axis=1 model_pattern[ii, :] = cosfunc_full(xdata_all, amp_mod[ii], frq_mod[ii], popt[0]) # Estimate the improvement of the effective read noise - tmp = outframe.copy() - tmp[osd_slice] -= model_pattern - mod_oscan, _ = rect_slice_with_mask(tmp, tmp_oscan, amp) - old_ron = astropy.stats.sigma_clipped_stats(overscan, sigma=5)[-1] - new_ron = astropy.stats.sigma_clipped_stats(overscan-mod_oscan, sigma=5)[-1] + full_model[osd_slice] = model_pattern + old_ron = astropy.stats.sigma_clipped_stats(overscan, sigma=5, stdfunc='mad_std')[-1] + new_ron = astropy.stats.sigma_clipped_stats(overscan-full_model[os_slice], sigma=5, stdfunc='mad_std')[-1] msgs.info(f'Effective read noise of amplifier {amp} reduced by a factor of {old_ron/new_ron:.2f}x') - # Subtract the model pattern from the full datasec - outframe[osd_slice] -= model_pattern - # Transpose if the input frame if applied along a different axis if axis == 0: - outframe = outframe.T - # Return the result - return outframe + return (outframe - full_model).T + return outframe - full_model def pattern_frequency(frame, axis=1): diff --git a/pypeit/find_objects.py b/pypeit/find_objects.py index 120a800bcc..5469a5ceba 100644 --- a/pypeit/find_objects.py +++ b/pypeit/find_objects.py @@ -972,22 +972,6 @@ class SlicerIFUFindObjects(MultiSlitFindObjects): def __init__(self, sciImg, slits, spectrograph, par, objtype, **kwargs): super().__init__(sciImg, slits, spectrograph, par, objtype, **kwargs) - def initialize_slits(self, slits, initial=True): - """ - Gather all the :class:`~pypeit.slittrace.SlitTraceSet` attributes that - we'll use here in :class:`FindObjects`. Identical to the parent but the - slits are not trimmed. - - Args: - slits (:class:`~pypeit.slittrace.SlitTraceSet`): - SlitTraceSet object containing the slit boundaries that will be - initialized. - initial (:obj:`bool`, optional): - Use the initial definition of the slits. If False, - tweaked slits are used. - """ - super().initialize_slits(slits, initial=True) - def global_skysub(self, skymask=None, bkg_redux_sciimg=None, update_crmask=True, previous_sky=None, show_fit=False, show=False, show_objs=False, objs_not_masked=False, reinit_bpm: bool = True): @@ -1074,7 +1058,7 @@ def joint_skysub(self, skymask=None, update_crmask=True, trim_edg=(0,0), # Use the FWHM map determined from the arc lines to convert the science frame # to have the same effective spectral resolution. - fwhm_map = self.wv_calib.build_fwhmimg(self.tilts, self.slits, initial=True, spat_flexure=self.spat_flexure_shift) + fwhm_map = self.wv_calib.build_fwhmimg(self.tilts, self.slits, spat_flexure=self.spat_flexure_shift) thismask = thismask & (fwhm_map != 0.0) # Need to include S/N for deconvolution sciimg = skysub.convolve_skymodel(self.sciImg.image, fwhm_map, thismask) @@ -1083,8 +1067,8 @@ def joint_skysub(self, skymask=None, update_crmask=True, trim_edg=(0,0), model_ivar = self.sciImg.ivar sl_ref = self.par['calibrations']['flatfield']['slit_illum_ref_idx'] # Prepare the slitmasks for the relative spectral illumination - slitmask = self.slits.slit_img(pad=0, initial=True, flexure=self.spat_flexure_shift) - slitmask_trim = self.slits.slit_img(pad=-3, initial=True, flexure=self.spat_flexure_shift) + slitmask = self.slits.slit_img(pad=0, flexure=self.spat_flexure_shift) + slitmask_trim = self.slits.slit_img(pad=-3, flexure=self.spat_flexure_shift) for nn in range(numiter): msgs.info("Performing iterative joint sky subtraction - ITERATION {0:d}/{1:d}".format(nn+1, numiter)) # TODO trim_edg is in the parset so it should be passed in here via trim_edg=tuple(self.par['reduce']['trim_edge']), diff --git a/pypeit/flatfield.py b/pypeit/flatfield.py index e73f8462ae..432ed4d469 100644 --- a/pypeit/flatfield.py +++ b/pypeit/flatfield.py @@ -803,6 +803,7 @@ def fit(self, spat_illum_only=False, doqa=True, debug=False): # Set parameters (for convenience; spec_samp_fine = self.flatpar['spec_samp_fine'] spec_samp_coarse = self.flatpar['spec_samp_coarse'] + tweak_method = self.flatpar['tweak_method'] tweak_slits = self.flatpar['tweak_slits'] tweak_slits_thresh = self.flatpar['tweak_slits_thresh'] tweak_slits_maxfrac = self.flatpar['tweak_slits_maxfrac'] @@ -1089,9 +1090,10 @@ def fit(self, spat_illum_only=False, doqa=True, debug=False): # TODO: Will this break if left_thresh, left_shift, self.slits.left_tweak[:,slit_idx], right_thresh, \ right_shift, self.slits.right_tweak[:,slit_idx] \ - = flat.tweak_slit_edges(self.slits.left_init[:,slit_idx], + = self.tweak_slit_edges(self.slits.left_init[:,slit_idx], self.slits.right_init[:,slit_idx], spat_coo_data, spat_flat_data, + method=tweak_method, thresh=tweak_slits_thresh, maxfrac=tweak_slits_maxfrac, debug=debug) # TODO: Because the padding doesn't consider adjacent @@ -1102,14 +1104,15 @@ def fit(self, spat_illum_only=False, doqa=True, debug=False): # Update the onslit mask _slitid_img = self.slits.slit_img(slitidx=slit_idx, initial=False) onslit_tweak = _slitid_img == slit_spat - spat_coo_tweak = self.slits.spatial_coordinate_image(slitidx=slit_idx, - slitid_img=_slitid_img) + # Note, we need to get the full image with the coordinates similar to spat_coo_init, otherwise, the + # tweaked locations are biased. + spat_coo_tweak = self.slits.spatial_coordinate_image(slitidx=slit_idx, full=True, slitid_img=_slitid_img) # Construct the empirical illumination profile # TODO This is extremely inefficient, because we only need to re-fit the illumflat, but # spatial_fit does both the reconstruction of the illumination function and the bspline fitting. # Only the b-spline fitting needs be reddone with the new tweaked spatial coordinates, so that would - # save a ton of runtime. It is not a trivial change becauase the coords are sorted, etc. + # save a ton of runtime. It is not a trivial change because the coords are sorted, etc. exit_status, spat_coo_data, spat_flat_data, spat_bspl, spat_gpm_fit, \ spat_flat_fit, spat_flat_data_raw = self.spatial_fit( norm_spec, spat_coo_tweak, median_slit_widths[slit_idx], spat_gpm, gpm, debug=False) @@ -1160,7 +1163,7 @@ def fit(self, spat_illum_only=False, doqa=True, debug=False): spat_model[onslit_padded] = spat_bspl.value(spat_coo_final[onslit_padded])[0] specspat_illum = np.fmax(spec_model, 1.0) * spat_model norm_spatspec = rawflat / specspat_illum - self.spatial_fit_finecorr(norm_spatspec, onslit_tweak, slit_idx, slit_spat, gpm, doqa=doqa) + spat_illum_fine = self.spatial_fit_finecorr(norm_spatspec, onslit_tweak, slit_idx, slit_spat, gpm, doqa=doqa)[onslit_tweak] # ---------------------------------------------------------- # Construct the illumination profile with the tweaked edges @@ -1401,6 +1404,11 @@ def spatial_fit_finecorr(self, normed, onslit_tweak, slit_idx, slit_spat, gpm, s A positive number trims the slit edges, a negative number pads the slit edges. doqa : :obj:`bool`, optional: Save the QA? + + Returns + ------- + illumflat_finecorr: `numpy.ndarray`_ + An image (same shape as normed) containing the fine correction to the spatial illumination profile """ # TODO :: Include fit_order in the parset?? fit_order = np.array([3, 6]) @@ -1413,7 +1421,7 @@ def spatial_fit_finecorr(self, normed, onslit_tweak, slit_idx, slit_spat, gpm, s onslit_tweak_trim = self.slits.slit_img(pad=-slit_trim, slitidx=slit_idx, initial=False) == slit_spat # Setup slitimg = (slit_spat + 1) * onslit_tweak.astype(int) - 1 # Need to +1 and -1 so that slitimg=-1 when off the slit - left, right, msk = self.slits.select_edges(initial=True, flexure=self.wavetilts.spat_flexure) + left, right, msk = self.slits.select_edges(flexure=self.wavetilts.spat_flexure) this_left = left[:, slit_idx] this_right = right[:, slit_idx] slitlen = int(np.median(this_right - this_left)) @@ -1422,7 +1430,6 @@ def spatial_fit_finecorr(self, normed, onslit_tweak, slit_idx, slit_spat, gpm, s this_slit = np.where(onslit_tweak & self.rawflatimg.select_flag(invert=True) & (self.waveimg!=0.0)) this_wave = self.waveimg[this_slit] xpos_img = self.slits.spatial_coordinate_image(slitidx=slit_idx, - initial=True, slitid_img=slitimg, flexure_shift=self.wavetilts.spat_flexure) # Generate the trimmed versions for fitting @@ -1465,7 +1472,7 @@ def spatial_fit_finecorr(self, normed, onslit_tweak, slit_idx, slit_spat, gpm, s normed[np.logical_not(onslit_tweak)] = 1 # For the QA, make everything off the slit equal to 1 spatillum_finecorr_qa(normed, illumflat_finecorr, this_left, this_right, ypos_fit, this_slit_trim, outfile=outfile, title=title, half_slen=slitlen//2) - return + return illumflat_finecorr def extract_structure(self, rawflat_orig, slit_trim=3): """ @@ -1574,7 +1581,72 @@ def spectral_illumination(self, gpm=None, debug=False): slit_illum_ref_idx=self.flatpar['slit_illum_ref_idx'], model=None, gpmask=gpm, skymask=None, trim=trim, flexure=self.wavetilts.spat_flexure, - smooth_npix=self.flatpar['slit_illum_smooth_npix']) + smooth_npix=self.flatpar['slit_illum_smooth_npix'], + debug=debug) + + def tweak_slit_edges(self, left, right, spat_coo, norm_flat, method='threshold', thresh=0.93, maxfrac=0.1, debug=False): + """ + Tweak the slit edges based on the normalized slit illumination profile. + + Args: + left (`numpy.ndarray`_): + Array with the left slit edge for a single slit. Shape is + :math:`(N_{\rm spec},)`. + right (`numpy.ndarray`_): + Array with the right slit edge for a single slit. Shape + is :math:`(N_{\rm spec},)`. + spat_coo (`numpy.ndarray`_): + Spatial pixel coordinates in fractions of the slit width + at each spectral row for the provided normalized flat + data. Coordinates are relative to the left edge (with the + left edge at 0.). Shape is :math:`(N_{\rm flat},)`. + Function assumes the coordinate array is sorted. + norm_flat (`numpy.ndarray`_) + Normalized flat data that provide the slit illumination + profile. Shape is :math:`(N_{\rm flat},)`. + method (:obj:`str`, optional): + Method to use for tweaking the slit edges. Options are: + - 'threshold': Use the threshold to set the slit edge + and then shift it to the left or right based on the + illumination profile. + - 'gradient': Use the gradient of the illumination + profile to set the slit edge and then shift it to + the left or right based on the illumination profile. + thresh (:obj:`float`, optional): + Threshold of the normalized flat profile at which to + place the two slit edges. + maxfrac (:obj:`float`, optional): + The maximum fraction of the slit width that the slit edge + can be adjusted by this algorithm. If ``maxfrac = 0.1``, + this means the maximum change in the slit width (either + narrowing or broadening) is 20% (i.e., 10% for either + edge). + debug (:obj:`bool`, optional): + Show flow interrupting plots that show illumination + profile in the case of a failure and the placement of the + tweaked edge for each side of the slit regardless. + + Returns: + tuple: Returns six objects: + + - The threshold used to set the left edge + - The fraction of the slit that the left edge is shifted to + the right + - The adjusted left edge + - The threshold used to set the right edge + - The fraction of the slit that the right edge is shifted to + the left + - The adjusted right edge + """ + # TODO :: Since this is just a wrapper, and not really "core", maybe it should be moved to pypeit.flatfield? + # Tweak the edges via the specified method + if method == "threshold": + return flat.tweak_slit_edges_threshold(left, right, spat_coo, norm_flat, + thresh=thresh, maxfrac=maxfrac, debug=debug) + elif method == "gradient": + return flat.tweak_slit_edges_gradient(left, right, spat_coo, norm_flat, maxfrac=maxfrac, debug=debug) + else: + msgs.error("Method for tweaking slit edges not recognized: {0}".format(method)) def spatillum_finecorr_qa(normed, finecorr, left, right, ypos, cut, outfile=None, title=None, half_slen=50): @@ -1801,7 +1873,7 @@ def show_flats(image_list, wcs_match=True, slits=None, waveimg=None): # TODO :: This could possibly be moved to core.flat def illum_profile_spectral(rawimg, waveimg, slits, slit_illum_ref_idx=0, smooth_npix=None, polydeg=None, - model=None, gpmask=None, skymask=None, trim=3, flexure=None, maxiter=5): + model=None, gpmask=None, skymask=None, trim=3, flexure=None, maxiter=5, debug=False): """ Determine the relative spectral illumination of all slits. Currently only used for image slicer IFUs. @@ -1834,6 +1906,8 @@ def illum_profile_spectral(rawimg, waveimg, slits, slit_illum_ref_idx=0, smooth_ Spatial flexure maxiter : :obj:`int` Maximum number of iterations to perform + debug : :obj:`bool` + Show the results of the relative spectral illumination correction Returns ------- @@ -1850,14 +1924,14 @@ def illum_profile_spectral(rawimg, waveimg, slits, slit_illum_ref_idx=0, smooth_ gpm = gpmask if (gpmask is not None) else np.ones_like(rawimg, dtype=bool) modelimg = model if (model is not None) else rawimg.copy() # Setup the slits - slitid_img_init = slits.slit_img(pad=0, initial=True, flexure=flexure) - slitid_img_trim = slits.slit_img(pad=-trim, initial=True, flexure=flexure) + slitid_img = slits.slit_img(pad=0, flexure=flexure) + slitid_img_trim = slits.slit_img(pad=-trim, flexure=flexure) scaleImg = np.ones_like(rawimg) modelimg_copy = modelimg.copy() # Obtain the minimum and maximum wavelength of all slits mnmx_wv = np.zeros((slits.nslits, 2)) for slit_idx, slit_spat in enumerate(slits.spat_id): - onslit_init = (slitid_img_init == slit_spat) + onslit_init = (slitid_img == slit_spat) mnmx_wv[slit_idx, 0] = np.min(waveimg[onslit_init]) mnmx_wv[slit_idx, 1] = np.max(waveimg[onslit_init]) wavecen = np.mean(mnmx_wv, axis=1) @@ -1892,7 +1966,7 @@ def illum_profile_spectral(rawimg, waveimg, slits, slit_illum_ref_idx=0, smooth_ for ss in range(1, slits.spat_id.size): # Calculate the region of overlap onslit_b = (slitid_img_trim == slits.spat_id[wvsrt[ss]]) - onslit_b_init = (slitid_img_init == slits.spat_id[wvsrt[ss]]) + onslit_b_init = (slitid_img == slits.spat_id[wvsrt[ss]]) onslit_b_olap = onslit_b & gpm & (waveimg >= mnmx_wv[wvsrt[ss], 0]) & (waveimg <= mnmx_wv[wvsrt[ss], 1]) & skymask_now hist, edge = np.histogram(waveimg[onslit_b_olap], bins=wavebins, weights=modelimg_copy[onslit_b_olap]) cntr, edge = np.histogram(waveimg[onslit_b_olap], bins=wavebins) @@ -1943,7 +2017,6 @@ def illum_profile_spectral(rawimg, waveimg, slits, slit_illum_ref_idx=0, smooth_ modelimg_copy /= relscl_model if max(abs(1/minv), abs(maxv)) < 1.005: # Relative accuracy of 0.5% is sufficient break - debug = False if debug: embed() ricp = rawimg.copy() @@ -1958,12 +2031,24 @@ def illum_profile_spectral(rawimg, waveimg, slits, slit_illum_ref_idx=0, smooth_ scale_ref = histScl * norm plt.subplot(211) plt.plot(wave_ref, this_spec) + plt.xlim([3600, 4500]) plt.subplot(212) plt.plot(wave_ref, scale_ref) + plt.xlim([3600, 4500]) plt.subplot(211) plt.plot(wave_ref, spec_ref, 'k--') + plt.xlim([3600, 4500]) + plt.show() + # Plot the relative scales of each slit + scales_med, scales_avg = np.zeros(slits.spat_id.size), np.zeros(slits.spat_id.size) + for ss in range(slits.spat_id.size): + onslit_ref_trim = (slitid_img_trim == slits.spat_id[ss]) & gpm & skymask_now & (waveimg>3628) & (waveimg<4510) + scales_med[ss] = np.median(ricp[onslit_ref_trim]/scaleImg[onslit_ref_trim]) + scales_avg[ss] = np.mean(ricp[onslit_ref_trim]/scaleImg[onslit_ref_trim]) + plt.plot(slits.spat_id, scales_med, 'bo-', label='Median') + plt.plot(slits.spat_id, scales_avg, 'ro-', label='Mean') + plt.legend() plt.show() - return scaleImg diff --git a/pypeit/images/rawimage.py b/pypeit/images/rawimage.py index 809187e5bb..4cb750c0b1 100644 --- a/pypeit/images/rawimage.py +++ b/pypeit/images/rawimage.py @@ -1218,7 +1218,7 @@ def subtract_scattlight(self, msscattlight, slits, debug=False): f" {tmp[13]}, {tmp[14]}, {tmp[15]}]) # Polynomial terms (coefficients of spec**index)\n" print(strprint) pad = msscattlight.pad // spatbin - offslitmask = slits.slit_img(pad=pad, initial=True, flexure=None) == -1 + offslitmask = slits.slit_img(pad=pad, flexure=None) == -1 from matplotlib import pyplot as plt _frame = self.image[ii, ...] vmin, vmax = 0, np.max(scatt_img) @@ -1243,7 +1243,7 @@ def subtract_scattlight(self, msscattlight, slits, debug=False): elif self.par["scattlight"]["method"] == "frame": # Calculate a model specific for this frame pad = msscattlight.pad // spatbin - offslitmask = slits.slit_img(pad=pad, initial=True, flexure=None) == -1 + offslitmask = slits.slit_img(pad=pad, flexure=None) == -1 # Get starting parameters for the scattered light model x0, bounds = self.spectrograph.scattered_light_archive(binning, dispname) # Perform a fit to the scattered light @@ -1267,11 +1267,11 @@ def subtract_scattlight(self, msscattlight, slits, debug=False): # Check if a fine correction to the scattered light should be applied if do_finecorr: pad = self.par['scattlight']['finecorr_pad'] // spatbin - offslitmask = slits.slit_img(pad=pad, initial=True, flexure=None) == -1 + offslitmask = slits.slit_img(pad=pad, flexure=None) == -1 # Check if the user wishes to mask some inter-slit regions if self.par['scattlight']['finecorr_mask'] is not None: # Get the central trace of each slit - left, right, _ = slits.select_edges(initial=True, flexure=None) + left, right, _ = slits.select_edges(flexure=None) centrace = 0.5*(left+right) # Now mask user-defined inter-slit regions offslitmask = scattlight.mask_slit_regions(offslitmask, centrace, diff --git a/pypeit/par/pypeitpar.py b/pypeit/par/pypeitpar.py index 4c1a21a316..c9de2f60b3 100644 --- a/pypeit/par/pypeitpar.py +++ b/pypeit/par/pypeitpar.py @@ -586,7 +586,7 @@ class FlatFieldPar(ParSet): see :ref:`parameters`. """ def __init__(self, method=None, pixelflat_file=None, spec_samp_fine=None, - spec_samp_coarse=None, spat_samp=None, tweak_slits=None, tweak_slits_thresh=None, + spec_samp_coarse=None, spat_samp=None, tweak_slits=None, tweak_method=None, tweak_slits_thresh=None, tweak_slits_maxfrac=None, rej_sticky=None, slit_trim=None, slit_illum_pad=None, illum_iter=None, illum_rej=None, twod_fit_npoly=None, saturated_slits=None, slit_illum_relative=None, slit_illum_ref_idx=None, slit_illum_smooth_npix=None, @@ -652,6 +652,17 @@ def __init__(self, method=None, pixelflat_file=None, spec_samp_fine=None, descr['tweak_slits'] = 'Use the illumination flat field to tweak the slit edges. ' \ 'This will work even if illumflatten is set to False ' + defaults['tweak_method'] = 'threshold' + options['tweak_method'] = FlatFieldPar.valid_tweak_methods() + dtypes['tweak_method'] = str + descr['tweak_method'] = 'Method used to tweak the slit edges (when "tweak_slits" is set to True). ' \ + 'Options include: {0:s}. '.format(', '.join(options['tweak_method'])) + \ + 'The "threshold" method determines when the left and right slit edges ' \ + 'fall below a threshold relative to the peak illumination. ' \ + 'The "gradient" method determines where the gradient is the highest at ' \ + 'the left and right slit edges. This method performs better when there is ' \ + 'systematic vignetting in the spatial direction. ' \ + defaults['tweak_slits_thresh'] = 0.93 dtypes['tweak_slits_thresh'] = float descr['tweak_slits_thresh'] = 'If tweak_slits is True, this sets the illumination function threshold used to ' \ @@ -769,7 +780,7 @@ def from_dict(cls, cfg): k = np.array([*cfg.keys()]) parkeys = ['method', 'pixelflat_file', 'spec_samp_fine', 'spec_samp_coarse', 'spat_samp', 'pixelflat_min_wave', 'pixelflat_max_wave', - 'tweak_slits', 'tweak_slits_thresh', 'tweak_slits_maxfrac', + 'tweak_slits', 'tweak_method', 'tweak_slits_thresh', 'tweak_slits_maxfrac', 'rej_sticky', 'slit_trim', 'slit_illum_pad', 'slit_illum_relative', 'illum_iter', 'illum_rej', 'twod_fit_npoly', 'saturated_slits', 'slit_illum_ref_idx', 'slit_illum_smooth_npix', 'slit_illum_finecorr', 'fit_2d_det_response'] @@ -790,6 +801,13 @@ def valid_methods(): """ return ['bspline', 'skip'] # [ 'PolyScan', 'bspline' ]. Same here. Not sure what PolyScan is + @staticmethod + def valid_tweak_methods(): + """ + Return the valid options for tweaking slits. + """ + return ['threshold', 'gradient'] + @staticmethod def valid_saturated_slits_methods(): """ diff --git a/pypeit/scattlight.py b/pypeit/scattlight.py index a6f149a6d4..675da0d07b 100644 --- a/pypeit/scattlight.py +++ b/pypeit/scattlight.py @@ -122,7 +122,7 @@ def show(self, image=None, slits=None, mask=False, wcs_match=True): wcs_match : :obj:`bool`, optional Use a reference image for the WCS and match all image in other channels to it. """ - offslitmask = slits.slit_img(pad=0, initial=True, flexure=None) == -1 if mask else 1 + offslitmask = slits.slit_img(pad=0, flexure=None) == -1 if mask else 1 # Prepare the frames _data = self.scattlight_raw if image is None else image diff --git a/pypeit/slittrace.py b/pypeit/slittrace.py index 93b4b7df32..9405b51ad2 100644 --- a/pypeit/slittrace.py +++ b/pypeit/slittrace.py @@ -479,7 +479,7 @@ def get_slitlengths(self, initial=False, median=False): slitlen = right - left return np.median(slitlen, axis=1) if median else slitlen - def get_radec_image(self, wcs, alignSplines, tilts, slit_compute=None, slice_offset=None, initial=True, flexure=None): + def get_radec_image(self, wcs, alignSplines, tilts, slit_compute=None, slice_offset=None, initial=False, flexure=None): """Generate an RA and DEC image for every pixel in the frame NOTE: This function is currently only used for SlicerIFU reductions. diff --git a/pypeit/spectrographs/gemini_gnirs.py b/pypeit/spectrographs/gemini_gnirs.py index dd56f426c3..5a88298847 100644 --- a/pypeit/spectrographs/gemini_gnirs.py +++ b/pypeit/spectrographs/gemini_gnirs.py @@ -612,11 +612,13 @@ def default_pypeit_par(cls): # Don't do 1D extraction for 3D data - it's meaningless because the DAR correction must be performed on the 3D data. par['reduce']['extraction']['skip_extraction'] = True # Because extraction occurs before the DAR correction, don't extract - #par['calibrations']['flatfield']['tweak_slits'] = False # Do not tweak the slit edges (we want to use the full slit) + # Tweak the slit edges using the gradient method for SlicerIFU + par['calibrations']['flatfield']['tweak_slits'] = True # Tweak the slit edges + par['calibrations']['flatfield']['tweak_method'] = 'gradient' # The gradient method is better for SlicerIFU. par['calibrations']['flatfield']['tweak_slits_thresh'] = 0.0 # Make sure the full slit is used (i.e. when the illumination fraction is > 0.5) par['calibrations']['flatfield']['tweak_slits_maxfrac'] = 0.0 # Make sure the full slit is used (i.e. no padding) par['calibrations']['flatfield']['slit_trim'] = 2 # Trim the slit edges - par['calibrations']['slitedges']['pad'] = 2 # Need to pad out the tilts for the astrometric transform when creating a datacube. + par['calibrations']['slitedges']['pad'] = 0 # Do not pad the slits - this ensures that the tweak_edges method=gradient guarantees that the edges are defined at the maximum gradient. # Decrease the wave tilts order, given the shorter slits of the IFU par['calibrations']['tilts']['spat_order'] = 1 @@ -716,7 +718,7 @@ def get_wcs(self, hdr, slits, platescale, wave0, dwv, spatial_scale=None): pxscl = spatial_scale / 3600.0 # 3600 is to convert arcsec to degrees # Get the typical slit length (this changes by ~0.3% over all slits, so a constant is fine for now) - slitlength = int(np.round(np.median(slits.get_slitlengths(initial=True, median=True)))) + slitlength = int(np.round(np.median(slits.get_slitlengths(median=True)))) # Get RA/DEC raval = self.get_meta_value([hdr], 'ra') diff --git a/pypeit/spectrographs/gtc_osiris.py b/pypeit/spectrographs/gtc_osiris.py index f622ad99de..78bc81fa21 100644 --- a/pypeit/spectrographs/gtc_osiris.py +++ b/pypeit/spectrographs/gtc_osiris.py @@ -471,6 +471,13 @@ def default_pypeit_par(cls): par['calibrations']['tilts']['spat_order'] = 1 par['calibrations']['tilts']['spec_order'] = 1 + # Tweak the slit edges using the gradient method for SlicerIFU + par['calibrations']['slitedges']['pad'] = 0 # Do not pad the slits - this ensures that the tweak_edges method=gradient guarantees that the edges are defined at the maximum gradient. + par['calibrations']['flatfield']['tweak_slits'] = True # Tweak the slit edges + par['calibrations']['flatfield']['tweak_method'] = 'gradient' # The gradient method is better for SlicerIFU. + par['calibrations']['flatfield']['tweak_slits_thresh'] = 0.0 # Make sure the full slit is used (i.e. when the illumination fraction is > 0.5) + par['calibrations']['flatfield']['tweak_slits_maxfrac'] = 0.0 # Make sure the full slit is used (i.e. no padding) + # Make sure that this is reduced as a slit (as opposed to fiber) spectrograph par['reduce']['cube']['slit_spec'] = True par['reduce']['cube']['combine'] = False # Make separate spec3d files from the input spec2d files @@ -527,7 +534,7 @@ def get_wcs(self, hdr, slits, platescale, wave0, dwv, spatial_scale=None): pxscl = spatial_scale / 3600.0 # 3600 is to convert arcsec to degrees # Get the typical slit length (this changes by ~0.3% over all slits, so a constant is fine for now) - slitlength = int(np.round(np.median(slits.get_slitlengths(initial=True, median=True)))) + slitlength = int(np.round(np.median(slits.get_slitlengths(median=True)))) # Get RA/DEC raval = self.get_meta_value([hdr], 'ra') diff --git a/pypeit/spectrographs/keck_kcwi.py b/pypeit/spectrographs/keck_kcwi.py index 74f2a7cff9..cdf4a20b10 100644 --- a/pypeit/spectrographs/keck_kcwi.py +++ b/pypeit/spectrographs/keck_kcwi.py @@ -283,6 +283,34 @@ def default_pypeit_par(cls): # Set the number of alignments in the align frames par['calibrations']['alignment']['locations'] = [0.1, 0.3, 0.5, 0.7, 0.9] # TODO:: Check this - is this accurate enough? + # Correct the illumflat for pixel-to-pixel sensitivity variations + par['calibrations']['illumflatframe']['process']['use_pixelflat'] = True + + # Make sure the overscan is subtracted from the dark + par['calibrations']['darkframe']['process']['use_overscan'] = True + + # Set the slit edge parameters + par['calibrations']['slitedges']['fit_order'] = 4 + par['calibrations']['slitedges']['pad'] = 0 # Do not pad the slits - this ensures that the tweak_edges method=gradient guarantees that the edges are defined at the maximum gradient. + par['calibrations']['slitedges']['edge_thresh'] = 5 # 5 works well with a range of setups tested by RJC (mostly 1x1 binning) + + # KCWI has non-uniform spectral resolution across the field-of-view + par['calibrations']['wavelengths']['fwhm_spec_order'] = 1 + par['calibrations']['wavelengths']['fwhm_spat_order'] = 2 + + # Alter the method used to combine pixel flats + par['calibrations']['pixelflatframe']['process']['combine'] = 'median' + par['calibrations']['flatfield']['spec_samp_coarse'] = 20.0 + par['calibrations']['flatfield']['tweak_slits'] = True # Tweak the slit edges + par['calibrations']['flatfield']['tweak_method'] = 'gradient' # The gradient method is better for SlicerIFU. + par['calibrations']['flatfield']['tweak_slits_thresh'] = 0.0 # Make sure the full slit is used (i.e. when the illumination fraction is > 0.5) + par['calibrations']['flatfield']['tweak_slits_maxfrac'] = 0.0 # Make sure the full slit is used (i.e. no padding) + par['calibrations']['flatfield']['slit_trim'] = 3 # Trim the slit edges + # Relative illumination correction + par['calibrations']['flatfield']['slit_illum_relative'] = True # Calculate the relative slit illumination + par['calibrations']['flatfield']['slit_illum_ref_idx'] = 14 # The reference index - this should probably be the same for the science frame + par['calibrations']['flatfield']['slit_illum_smooth_npix'] = 10 # Sufficiently small value so less structure in relative weights + # LACosmics parameters par['scienceframe']['process']['sigclip'] = 4.0 par['scienceframe']['process']['objlim'] = 1.5 @@ -607,7 +635,7 @@ def get_wcs(self, hdr, slits, platescale, wave0, dwv, spatial_scale=None): pxscl = spatial_scale / 3600.0 # 3600 is to convert arcsec to degrees # Get the typical slit length (this changes by ~0.3% over all slits, so a constant is fine for now) - slitlength = int(np.round(np.median(slits.get_slitlengths(initial=True, median=True)))) + slitlength = int(np.round(np.median(slits.get_slitlengths(median=True)))) # Get RA/DEC ra = self.compound_meta([hdr], 'ra') @@ -920,7 +948,7 @@ def default_pypeit_par(cls): par['calibrations']['scattlight_pad'] = 6 # This is the unbinned number of pixels to pad par['calibrations']['pixelflatframe']['process']['subtract_scattlight'] = True par['calibrations']['illumflatframe']['process']['subtract_scattlight'] = True - par['scienceframe']['process']['subtract_scattlight'] = True + par['scienceframe']['process']['subtract_scattlight'] = False par['scienceframe']['process']['scattlight']['finecorr_method'] = 'median' par['scienceframe']['process']['scattlight']['finecorr_pad'] = 4 # This is the unbinned number of pixels to pad par['scienceframe']['process']['scattlight']['finecorr_order'] = 2 @@ -935,33 +963,11 @@ def default_pypeit_par(cls): par['calibrations']['standardframe']['process']['correct_nonlinear'] = nonlin_array par['scienceframe']['process']['correct_nonlinear'] = nonlin_array - # Correct the illumflat for pixel-to-pixel sensitivity variations - par['calibrations']['illumflatframe']['process']['use_pixelflat'] = True - - # Make sure the overscan is subtracted from the dark - par['calibrations']['darkframe']['process']['use_overscan'] = True - - # Set the slit edge parameters - par['calibrations']['slitedges']['fit_order'] = 4 - par['calibrations']['slitedges']['pad'] = 2 # Need to pad out the tilts for the astrometric transform when creating a datacube. - par['calibrations']['slitedges']['edge_thresh'] = 5 # 5 works well with a range of setups tested by RJC (mostly 1x1 binning) - - # KCWI has non-uniform spectral resolution across the field-of-view - par['calibrations']['wavelengths']['fwhm_spec_order'] = 1 - par['calibrations']['wavelengths']['fwhm_spat_order'] = 2 - # Alter the method used to combine pixel flats - par['calibrations']['pixelflatframe']['process']['combine'] = 'median' - par['calibrations']['flatfield']['spec_samp_coarse'] = 20.0 + par['calibrations']['pixelflatframe']['process']['combine'] = 'mean' par['calibrations']['flatfield']['spat_samp'] = 1.0 # This should give 1% accuracy in the spatial illumination correction for 2x2 binning, and <0.5% accuracy for 1x1 binning - #par['calibrations']['flatfield']['tweak_slits'] = False # Do not tweak the slit edges (we want to use the full slit) - par['calibrations']['flatfield']['tweak_slits_thresh'] = 0.0 # Make sure the full slit is used (i.e. when the illumination fraction is > 0.5) - par['calibrations']['flatfield']['tweak_slits_maxfrac'] = 0.0 # Make sure the full slit is used (i.e. no padding) - par['calibrations']['flatfield']['slit_trim'] = 3 # Trim the slit edges - # Relative illumination correction - par['calibrations']['flatfield']['slit_illum_relative'] = True # Calculate the relative slit illumination - par['calibrations']['flatfield']['slit_illum_ref_idx'] = 14 # The reference index - this should probably be the same for the science frame - par['calibrations']['flatfield']['slit_illum_smooth_npix'] = 5 # Sufficiently small value so less structure in relative weights + + # Need to fit sinusoidal sensitivity pattern, and include in the relative pixel response par['calibrations']['flatfield']['fit_2d_det_response'] = True # Include the 2D detector response in the pixelflat. return par @@ -1351,33 +1357,8 @@ def default_pypeit_par(cls): par['calibrations']['standardframe']['process']['use_pattern'] = False par['scienceframe']['process']['use_pattern'] = False - # Correct the illumflat for pixel-to-pixel sensitivity variations - par['calibrations']['illumflatframe']['process']['use_pixelflat'] = True - - # Make sure the overscan is subtracted from the dark - par['calibrations']['darkframe']['process']['use_overscan'] = True - - # Set the slit edge parameters - par['calibrations']['slitedges']['fit_order'] = 4 - par['calibrations']['slitedges']['pad'] = 2 # Need to pad out the tilts for the astrometric transform when creating a datacube. - par['calibrations']['slitedges']['edge_thresh'] = 5 # 5 works well with a range of setups tested by RJC (mostly 1x1 binning) - - # KCWI has non-uniform spectral resolution across the field-of-view - par['calibrations']['wavelengths']['fwhm_spec_order'] = 1 - par['calibrations']['wavelengths']['fwhm_spat_order'] = 2 - - # Alter the method used to combine pixel flats - par['calibrations']['pixelflatframe']['process']['combine'] = 'median' - par['calibrations']['flatfield']['spec_samp_coarse'] = 20.0 - #par['calibrations']['flatfield']['tweak_slits'] = False # Do not tweak the slit edges (we want to use the full slit) - par['calibrations']['flatfield']['tweak_slits_thresh'] = 0.0 # Make sure the full slit is used (i.e. when the illumination fraction is > 0.5) - par['calibrations']['flatfield']['tweak_slits_maxfrac'] = 0.0 # Make sure the full slit is used (i.e. no padding) - par['calibrations']['flatfield']['slit_trim'] = 3 # Trim the slit edges - # Relative illumination correction - par['calibrations']['flatfield']['slit_illum_relative'] = True # Calculate the relative slit illumination - par['calibrations']['flatfield']['slit_illum_ref_idx'] = 14 # The reference index - this should probably be the same for the science frame - par['calibrations']['flatfield']['slit_illum_smooth_npix'] = 5 # Sufficiently small value so less structure in relative weights - par['calibrations']['flatfield']['fit_2d_det_response'] = True # Include the 2D detector response in the pixelflat. + # This is probably not necessary for KCRM + par['calibrations']['flatfield']['fit_2d_det_response'] = False # Include the 2D detector response in the pixelflat. # Sky subtraction parameters par['reduce']['skysub']['bspline_spacing'] = 0.4 diff --git a/pypeit/utils.py b/pypeit/utils.py index fbc6550c15..0a922ce18c 100644 --- a/pypeit/utils.py +++ b/pypeit/utils.py @@ -1380,6 +1380,38 @@ def find_nearest(array, values): return idxs +def linear_interpolate(x1, y1, x2, y2, x): + r""" + Interplate or extrapolate between two points. + + Given a line defined two points, :math:`(x_1,y_1)` and + :math:`(x_2,y_2)`, return the :math:`y` value of a new point on + the line at coordinate :math:`x`. + + This function is meant for speed. No type checking is performed and + the only check is that the two provided ordinate coordinates are not + numerically identical. By definition, the function will extrapolate + without any warning. + + Args: + x1 (:obj:`float`): + First abscissa position + y1 (:obj:`float`): + First ordinate position + x2 (:obj:`float`): + Second abscissa position + y3 (:obj:`float`): + Second ordinate position + x (:obj:`float`): + Abcissa for new value + + Returns: + :obj:`float`: Interpolated/extrapolated value of ordinate at + :math:`x`. + """ + return y1 if np.isclose(x1,x2) else y1 + (x-x1)*(y2-y1)/(x2-x1) + + def replace_bad(frame, bpm): """ Find all bad pixels, and replace the bad pixels with the nearest good pixel