pypeit.core.skysub module

Module for sky subtraction

pypeit.core.skysub.convolve_skymodel(input_img, fwhm_map, thismask, subpixel=5, nsample=10)[source]

Convolve an input image with a Gaussian function that ensures all pixels have the same FWHM (i.e. the returned image has a spectral resolution corresponding to the maximum FWHM specified in the input fwhm_map). To speed up the computation, the input image is uniformly convolved by a grid of FWHM values, and the final pixel-by-pixel convolved map is an interpolation of the grid point values.

Args:

input_img (numpy.ndarray):
The science frame, shape = nspec, nspat

fwhm_map (numpy.ndarray):
An array (same shape as input_img), that specifies the FWHM at every pixel in the image.

thismask (numpy.ndarray):
A boolean mask (True = good), same shape as input_img, of the detector pixels that fall in a slit and have a measured FWHM.

subpixel (int, optional):
Divide each pixel into this many subpixels to improve the accuracy of the convolution (a higher number gives a more accurate result, at the expense of computational time).

nsample (int, optional):
Number of grid points that will be used to evaluate the convolved image.

Returns:: The convolved input_img, same shape as input_img.
Return type:: numpy.ndarray

pypeit.core.skysub.ech_local_skysub_extract(sciimg, sciivar, fullmask, tilts, waveimg, global_sky, left, right, slitmask, sobjs, spat_pix=None, fit_fwhm=False, min_snr=2.0, bsp=0.6, trim_edg=(3, 3), std=False, prof_nsigma=None, use_2dmodel_mask=True, niter=4, sigrej=3.5, bkpts_optimal=True, force_gauss=False, sn_gauss=4.0, model_full_slit=False, model_noise=True, debug_bkpts=False, show_profile=False, show_resids=False, show_fwhm=False, adderr=0.01, base_var=None, count_scale=None)[source]

Perform local sky subtraction, profile fitting, and optimal extraction slit by slit. Objects are sky/subtracted extracted in order of the highest average (across all orders) S/N ratio object first, and then for a given object the highest S/N ratio orders are extracted first. The profile fitting FWHM are stored and progressively fit as the objects are extracted to properly ensure that low S/N orders use Gaussian extracted FWHMs from higher-S/N orders (i.e. in the regime where the data is too noisy for a non-parametric object profile fit). The FWHM of higher S/N ratio objects are used for lower S/N ratio objects (note this assumes point sources with FWHM set by the seeing).

Note on masking: This routine requires that all masking be performed in the upstream calling routine (Extract) and thus the left and right slit edge arrays must only contain these slits. Similarly, the sobjs object must only include the unmasked (good) objects that are to be extracted. The number of sobjs objects must equal to an integer multiple of the number of good slits/orders. The routine will fault if any of these criteria are not met.

Parameters:

sciimg (numpy.ndarray) – science image, usually with a global sky subtracted. shape = (nspec, nspat)
sciivar (numpy.ndarray) – inverse variance of science image. shape = (nspec, nspat)
fullmask (ImageBitMaskArray) – Image bitmask array.
tilts (numpy.ndarray) – spectral tilts. shape=(nspec, nspat)
waveimg (numpy.ndarray) – 2-d wavelength map
global_sky (numpy.ndarray) – Global sky model produced by global_skysub
left (numpy.ndarray) – Spatial-pixel coordinates for the left edges of each order. Shape = (nspec, norders)
right (numpy.ndarray) – Spatial-pixel coordinates for the right edges of each order. Shape = (nspec, norders)
slitmask (numpy.ndarray) – Image identifying the 0-indexed order associated with each pixel. Pixels with -1 are not associatead with any order.
sobjs (SpecObjs object) – Object containing the information about the objects found on the slit/order from objfind or ech_objfind
spat_pix (numpy.ndarray, optional) – Image containing the spatial location of pixels. If not input, it will be computed from spat_img = np.outer(np.ones(nspec), np.arange(nspat)). This option should generally not be used unless one is extracting 2d coadds for which a rectified image contains sub-pixel spatial information. shape=(nspec, nspat)
fit_fwhm (bool, optional) – if True, perform a fit to the FWHM of the object profiles to use for non-detected sources
min_snr (float, optional) – FILL IN
bsp (float, default = 0.6) – Break point spacing in pixels for the b-spline sky subtraction.
trim_edg (tuple of ints of floats, default = (3,3)) – Number of pixels to be ignored on the (left,right) edges of the slit in object/sky model fits.
std (bool, default = False) – This should be set to True if the object being extracted is a standards star so that the reduction parameters can be adjusted accordingly.
prof_nsigma (int or float, default = None) – Number of sigmas that the object profile will be fit, i.e. the region extending from -prof_nsigma to +prof_nsigma will be fit where sigma = FWHM/2.35. This option should only be used for bright large extended source with tails in their light profile like elliptical galaxies. If prof_nsigma is set then the profiles will no longer be apodized by an exponential at large distances from the trace.
use_2dmodel_mask (bool, optional) – Use the mask made from profile fitting when extracting?
niter (int, optional) – Number of iterations for successive profile fitting and local sky-subtraction
sigrej (float, optional) – Outlier rejection threshold for sky and object fitting Set by par[‘scienceimage’][‘skysub’][‘sky_sigrej’]
bkpts_optimal (bool, optional) –
Parameter governing whether spectral direction breakpoints for b-spline sky/object modeling are determined optimally. If bkpts_optima=True, the optimal break-point spacing will be determined directly using the optimal_bkpts function by measuring how well we are sampling the sky using piximg = (nspec-1)*yilyd. The bsp parameter in this case corresponds to the minimum distance between breakpoints which we allow. If bkpts_optimal = False, the break-points will be chosen to have a uniform spacing in pixel units sets by the bsp parameter, i.e. using the bkspace functionality of the bspline class:
```
bset = bspline.bspline(piximg_values, nord=4, bkspace=bsp)
fullbkpt = bset.breakpoints
```
force_gauss (bool, default = False) – If True, a Gaussian profile will always be assumed for the optimal extraction using the FWHM determined from object finding (or provided by the user) for the spatial profile.
sn_gauss (int or float, default = 4.0) – The signal to noise threshold above which optimal extraction with non-parametric b-spline fits to the objects spatial profile will be performed. For objects with median S/N < sn_gauss, a Gaussian profile will simply be assumed because there is not enough S/N to justify performing a more complicated fit.
model_full_slit (bool, default = False) – Set the maskwidth of the objects to be equal to the slit width/2 such that the entire slit will be modeled by the local skysubtraction. This mode is recommended for echelle spectra with reasonably narrow slits.
model_noise (bool, default = True) – If True, construct and iteratively update a model inverse variance image using variance_model(). Construction of the model variance requires base_var, and will use the provided values or defaults for the remaining variance_model() parameters. If False, a variance model will not be created and instead the input sciivar will always be taken to be the inverse variance. Note that in order for the noise model to make any sense one needs to be subtracting the sky and not the sky residuals. In other words, for near-IR reductions where difference imaging has been performed and this algorithm is used to fit out the sky residuals (but not the sky itself) one should definitely set model_noise=False since otherwise the code will attempt to create a noise model using sky residuals instead of the sky, which is incorrect (does not have the right count levels). In principle this could be improved if the user could pass in a model of what the sky is for near-IR difference imaging + residual subtraction
debug_bkpts (bool, optional) – debug
show_profile (bool, default=False) – Show QA for the object profile fitting to the screen. Note that this will show interactive matplotlib plots which will block the execution of the code until the window is closed.
show_resids (bool, optional) – Show the model fits and residuals.
show_fwhm (bool, optional) – show fwhm
adderr (float, default = 0.01) – Error floor. The quantity adderr**2*sciframe**2 is added to the variance to ensure that the S/N is never > 1/adderr, effectively setting a floor on the noise or a ceiling on the S/N. This is one of the components needed to construct the model variance (this is the same as noise_floor in variance_model()); see model_noise.
base_var (numpy.ndarray, shape is (nspec, nspat), optional) – The “base-level” variance in the data, set by the detector properties and the image processing steps. See base_variance().
count_scale (float, numpy.ndarray, optional) – A scale factor that has already been applied to the provided science image. It accounts for the number of frames contributing to the provided counts, and the relative throughput factors that can be measured from flat-field frames. For example, if the image has been flat-field corrected, this is the inverse of the flat-field counts. If None, set to 1. If a single float, assumed to be constant across the full image. If an array, the shape must match base_var. The variance will be 0 wherever this array is not positive, modulo the provided adderr. This is one of the components needed to construct the model variance; see model_noise.

Returns:

skyimage (numpy.ndarray) – Model sky flux where thismask is true.
objimage (numpy.ndarray) – Model object flux where thismask is true.
ivarmodel (numpy.ndarray) – Model inverse variance where thismask is true.
outmask (numpy.ndarray) – Model mask where thismask is true.
sobjs (SpecObjs object) – Same object as passed in

pypeit.core.skysub.generate_mask(pypeline, skyreg, slits, slits_left, slits_right, spat_flexure=None)[source]

Generate the mask of sky regions

Parameters:

pypeline (str) – Name of the pypeline being used (e.g. MultiSlit, Echelle, IFU, …)
skyreg (list) – A list of size nslits. Each element contains a numpy array (dtype=bool) where a True value indicates a value that is part of the sky region.
slits (SlitTraceSet) – Data container with slit trace information
slits_left (numpy.ndarray) – A 2D array containing the pixel coordinates of the left slit edges
slits_right (numpy.ndarray) – A 2D array containing the pixel coordinates of the right slit edges
resolution (int, optional) – The percentage regions will be scaled to the specified resolution. The resolution should probably correspond to the number of spatial pixels on the slit.

Returns:

mask – Boolean mask containing sky regions

Return type:

numpy.ndarray

pypeit.core.skysub.global_skysub(image, ivar, tilts, thismask, slit_left, slit_righ, inmask=None, bsp=0.6, sigrej=3.0, maxiter=35, trim_edg=(3, 3), pos_mask=True, max_mask_frac=0.8, show_fit=False, no_poly=False, npoly=None)[source]

Perform global sky subtraction on an input slit THIS NEEDS MORE DESCRIPTION

Parameters:

image (numpy.ndarray) – Frame to be sky subtracted. float ndarray, shape (nspec, nspat)
ivar (numpy.ndarray) – Inverse variance image. float ndarray, shape (nspec, nspat)
tilts (numpy.ndarray) – Tilts indicating how wavelengths move across the slit. float ndarray, shape (nspec, nspat)
thismask (numpy.ndarray) – Specifies pixels in the slit in question. boolean array, shape (nspec, nspat)
slit_left (numpy.ndarray) – Left slit boundary in floating point pixels. shape (nspec, 1) or (nspec)
slit_righ (numpy.ndarray) – Right slit boundary in floating point pixels. shape (nspec, 1) or (nspec)
inmask (numpy.ndarray) – boolean ndarray, shape (nspec, nspat) Input mask for pixels not to be included in sky subtraction fits. True = Good (not masked), False = Bad (masked)
bsp (float, optional) – break point spacing in pixel units
sigrej (float, optional) – sigma rejection threshold
trim_edg (tuple, optional) – floats (left_edge, right_edge) that indicate how many pixels to trim from left and right slit edges for creating the edgemask. These pixels are excluded from sky subtraction fits.
pos_mask (bool, optional) – First do a prelimnary fit to the log of the sky (i.e. positive pixels only). Then use this fit to create an input mask from the residuals lmask = (res < 5.0) & (res > -4.0) for the full fit. NOTE: pos_mask should be False for near-IR sky residual subtraction, since fitting the log(sky) requires that the counts are positive which will not be the case for i.e. an A-B image. Thus the routine will fail if pos_mask is not set to False.
max_mask_frac (float, optional) – Maximum fraction of total pixels that can be masked by the input masks. If more than this threshold is masked the code will return zeros and throw a warning.
show_fit (bool, optional) – If true, plot a fit of the sky pixels and model fit to the screen. This feature will block further execution until the screen is closed.
no_poly (bool, optional) – If True, do not incldue polynomial basis
npoly (int, optional) – Order of polynomial to use for the polynomial in the bspline Only used if no_poly=False

Returns:

The model sky background at the pixels where thismask is True:

>>>  skyframe = np.zeros_like(image)
>>>  thismask = slitpix == thisslit
>>>  skyframe[thismask] = global_skysub(image,ivar, tilts, thismask, slit_left, slit_righ)

Return type:

numpy.ndarray

pypeit.core.skysub.local_skysub_extract(sciimg, sciivar, tilts, waveimg, global_sky, thismask, slit_left, slit_righ, sobjs, ingpm=None, fwhmimg=None, spat_pix=None, adderr=0.01, bsp=0.6, trim_edg=(3, 3), std=False, prof_nsigma=None, niter=4, extract_good_frac=0.005, sigrej=3.5, bkpts_optimal=True, debug_bkpts=False, force_gauss=False, sn_gauss=4.0, model_full_slit=False, model_noise=True, show_profile=False, show_resids=False, use_2dmodel_mask=True, no_local_sky=False, base_var=None, count_scale=None)[source]

Perform local sky subtraction and extraction

IMPROVE THIS DOCSTRING

Parameters:

sciimg (numpy.ndarray) – science image, usually with a global sky subtracted. shape = (nspec, nspat)
sciivar (numpy.ndarray) – inverse variance of science image. shape = (nspec, nspat)
tilts (numpy.ndarray) – spectral tilts. shape=(nspec, nspat)
waveimg (numpy.ndarray) – 2-d wavelength map
global_sky (numpy.ndarray) – Global sky model produced by global_skysub
thismask (numpy.ndarray) – Specifies pixels in the slit in question
slit_left (numpy.ndarray) – Left slit boundary in floating point pixels. shape (nspec, 1) or (nspec)
slit_righ (numpy.ndarray) – Right slit boundary in floating point pixels. shape (nspec, 1) or (nspec)
sobjs (SpecObjs object) – Object containing the information about the objects found on the slit/order from objfind or ech_objfind
ingpm (numpy.ndarray, optional) – Input mask with any non-zero item flagged as False using pypeit.images.imagebitmask.ImageBitMask shape=(nspec, nspat)
fwhmimg (numpy.ndarray, None, optional:) – Floating-point image containing the modeled spectral FWHM (in pixels) at every pixel location. Must have the same shape as sciimg, \((N_{\rm spec}, N_{\rm spat})\).
spat_pix (numpy.ndarray, optional) – Image containing the spatial location of pixels. If not input, it will be computed from spat_img = np.outer(np.ones(nspec), np.arange(nspat)). This option should generally not be used unless one is extracting 2d coadds for which a rectified image contains sub-pixel spatial information. shape (nspec, nspat)
adderr (float, default = 0.01) – Error floor. The quantity adderr**2*sciframe**2 is added to the variance to ensure that the S/N is never > 1/adderr, effectively setting a floor on the noise or a ceiling on the S/N. This is one of the components needed to construct the model variance (this is the same as noise_floor in variance_model()); see model_noise.
bsp (float, default = 0.6) – Break point spacing in pixels for the b-spline sky subtraction.
trim_edg (tuple of ints of floats, default = (3,3)) – Number of pixels to be ignored on the (left,right) edges of the slit in object/sky model fits.
std (bool, default = False) – This should be set to True if the object being extracted is a standards star so that the reduction parameters can be adjusted accordingly.
prof_nsigma (int or float, default = None) – Number of sigmas that the object profile will be fit, i.e. the region extending from -prof_nsigma to +prof_nsigma will be fit where sigma = FWHM/2.35. This option should only be used for bright large extended source with tails in their light profile like elliptical galaxies. If prof_nsigma is set then the profiles will no longer be apodized by an exponential at large distances from the trace.
niter (int, default = 4) – Number of iterations for successive profile fitting and local sky-subtraction
extract_good_frac (float, default = 0.005) – Minimum fraction of pixels along the spectral direction with good optimal extraction
sigrej (float, optional) – Outlier rejection threshold for sky and object fitting Set by par[‘scienceimage’][‘skysub’][‘sky_sigrej’]
bkpts_optimal (bool, optional) –
Parameter governing whether spectral direction breakpoints for b-spline sky/object modeling are determined optimally. If bkpts_optimal=True, the optimal break-point spacing will be determined directly using the optimal_bkpts function by measuring how well we are sampling the sky using piximg = (nspec-1)*yilyd. The bsp parameter in this case corresponds to the minimum distance between breakpoints which we allow. If bkpts_optimal = False, the break-points will be chosen to have a uniform spacing in pixel units sets by the bsp parameter, i.e. using the bkspace functionality of the bspline class:
```
bset = bspline.bspline(piximg_values, nord=4, bkspace=bsp)
fullbkpt = bset.breakpoints
```
debug_bkpts (bool, default=False) – Make an interactive plot to the screen to indicate how the breakpoints are being chosen.
force_gauss (bool, default = False) – If True, a Gaussian profile will always be assumed for the optimal extraction using the FWHM determined from object finding (or provided by the user) for the spatial profile.
sn_gauss (int or float, default = 4.0) – The signal to noise threshold above which optimal extraction with non-parametric b-spline fits to the objects spatial profile will be performed. For objects with median S/N < sn_gauss, a Gaussian profile will simply be assumed because there is not enough S/N to justify performing a more complicated fit.
model_full_slit (bool, default = False) – Set the maskwidth of the objects to be equal to the slit width/2 such that the entire slit will be modeled by the local skysubtraction. This mode is recommended for echelle spectra with reasonably narrow slits.
model_noise (bool, default = True) – If True, construct and iteratively update a model inverse variance image using variance_model(). Construction of the model variance requires base_var, and will use the provided values or defaults for the remaining variance_model() parameters. If False, a variance model will not be created and instead the input sciivar will always be taken to be the inverse variance. Note that in order for the noise model to make any sense one needs to be subtracting the sky and not the sky residuals. In other words, for near-IR reductions where difference imaging has been performed and this algorithm is used to fit out the sky residuals (but not the sky itself) one should definitely set model_noise=False since otherwise the code will attempt to create a noise model using sky residuals instead of the sky, which is incorrect (does not have the right count levels). In principle this could be improved if the user could pass in a model of what the sky is for near-IR difference imaging + residual subtraction
show_profile (bool, default=False) – Show QA for the object profile fitting to the screen. Note that this will show interactive matplotlib plots which will block the execution of the code until the window is closed.
show_resids (bool, optional) – Show the model fits and residuals.
use_2dmodel_mask (bool, optional) – Use the mask made from profile fitting when extracting?
no_local_sky (bool, optional) – If True, do not fit local sky model, only object profile and extract optimally The objimage will be all zeros.
base_var (numpy.ndarray, shape is (nspec, nspat), optional) – The “base-level” variance in the data set by the detector properties and the image processing steps. See base_variance().
count_scale (float, numpy.ndarray, optional) – A scale factor, \(s\), that has already been applied to the provided science image. It accounts for the number of frames contributing to the provided counts, and the relative throughput factors that can be measured from flat-field frames. For example, if the image has been flat-field corrected, this is the inverse of the flat-field counts. If None, set to 1. If a single float, assumed to be constant across the full image. If an array, the shape must match base_var. The variance will be 0 wherever \(s \leq 0\), modulo the provided adderr. This is one of the components needed to construct the model variance; see model_noise.

Returns:

skyimage (numpy.ndarray) – Model sky flux where thismask is true.
objimage (numpy.ndarray) – Model object flux where thismask is true.
modelivar (numpy.ndarray) – Model inverse variance where thismask is true.
outmask (ImageBitMaskArray) – Copy of fullmask but with added flags were the image was extracted.

pypeit.core.skysub.optimal_bkpts(bkpts_optimal, bsp_min, piximg, sampmask, samp_frac=0.8, skyimage=None, min_spat=None, max_spat=None, debug=False)[source]

Generate an array of optimally spaced breakpoints for the global sky subtraction algorithm.

Parameters:

bkpts_optimal (bool) – If True, then the breakpoints are optimally spaced. If False, then the breakpoints are spaced uniformly.
bsp_min (float) – Desired B-spline breakpoint spacing in pixels
piximg (numpy.ndarray) – Image containing the pixel sampling, i.e. (nspec-1)*tilts. shape = (nspec, nspat)
sampmask (numpy.ndarray) – Boolean array indicating the pixels for which the B-spline fit will actually be evaluated. True = Good, False=Bad
samp_frac (float) – The fraction of spectral direction pixels required to have a sampling difference < bsp_min in order to instead adopt a uniform break point spacing, rather adopting the optimally spaced breakpoints.
skyimage (numpy.ndarray) – Sky model image used only for QA. shape = (nspec, nspat)
min_spat (float, optional) – Minimum spatial pixel used for local sky subtraction fitting. Only used for title of QA plot.
max_spat (float, optional) – Maximum spatial pixel used for local sky subtraction fitting. Only used for title of QA plot.
debug (bool, optional) – Show QA plot to debug breakpoint placing.

Returns:

fullbkpt – Locations of the optimally sampled breakpoints

Return type:

numpy.ndarray

pypeit.core.skysub.read_userregions(skyreg, nslits, maxslitlength)[source]

Parse the sky regions defined by the user. The text should be a comma separated list of percentages to apply to all slits.

Example

The string ':10,35:65,80:' would select (in all slits):

the leftmost 10% of the slit length,

the inner 30% (from 35-65% of the slit length), and

the final 20% of the slit length (from 80-100% of the slit length)

Parameters:

skyreg (str) – The sky region definition.
nslits (int) – Number of slits on the detector
maxslitlength (float) – The maximum slit length (in pixels).

Returns:

status (int) – Status of the region parsing (0 = Successful, 1,2 = fail)
regions (list) – A list of size nslits. Each element contains a numpy array (dtype=bool) of size resolution. A True value indicates a value that is part of the sky region.

pypeit.core.skysub.skyoptimal(piximg, data, ivar, oprof, sigrej=3.0, npoly=1, spatial_img=None, fullbkpt=None)[source]

Utility routine used by local_skysub_extract that performs the joint b-spline fit for sky-background and object flux.

Parameters:

piximg (numpy.ndarray) – piximg is tilts*(nspec-1) where nspec is the number of pixels in the spectral direction of the raw image. This is a wavelength in image coordinates which acts as the independent variable for sky and object model fits. This is 1d array (flattened in the calling routine) with shape= (nflat,).
data (numpy.ndarray) – science data that is being fit. Same shape as piximg.
ivar (numpy.ndarray) – inverse variance of science data that is being fit. Same shape as piximg.
oprof (numpy.ndarray) – Flattened object profiles for the data that is being fit. Shape = (nflat, nobj) where nobj is the number of objects being simultaneously fit. In other words, there are nobj object profiles.
sigrej (float, optional) – Sigma threshold for outlier rejection.
npoly (int, optional) – Order of polynomaial for the sky-background basis function. If spatial_img is passed in a fit with two independent variables will be performed (spectral described by piximg, and spatial direction described by spatia_img) and a legendre polynomial basis of order npoly will be used for the spatial direction. If npoly=1 or if spatial_img is not passed, a flat spatial profile basis funciton will instead be used.
spatial_img (numpy.ndarray, optional) – Image of the spatial coordinates of each pixel in the image used for 2d fitting. Same shape as piximg.
fullbkpt (numpy.ndarray, optional) – A 1d float array containing the breakpoints to be used for the B-spline fit. The breakpoints are arranged in the spectral direction, i.e. along the directino of the piximg independent variable.

Returns:

sky_bmodel (numpy.ndarray) – Array with same shape as piximg containing the B-spline model of the sky.
obj_bmodel (numpy.ndarray) – Array with same shape as piximg containing the B-spline model of the object flux.
gpm (numpy.ndarray) – Boolean good pixel mask array with the same shape as piximg indicating whether a pixel is good (True) or was masked (False).

pypeit.core.skysub.skysub_npoly(thismask)[source]

Utility routine used by global_skysub and local_skysub_extract. Determine the order for the spatial polynomial for global sky subtraction and local sky subtraction.

Parameters:: thismask (numpy.ndarray) – bool mask of shape (nspec, nspat) which specifies pixels in the slit in question
Returns:: Order of polynomial
Return type:: int