Sky Subtraction

Overview

This document describes how PypeIt performs sky subtraction.

See SkySubPar Keywords for the complete list of options related to sky subtraction.

Global

Phase I of sky subtraction is to perform a fit to the sky across the entire slit. By default, this is done twice: once without any knowledge of objects in the slit and then again after object detection has taken place (these are masked).

Default masking of objects is relatively benign. The FWHM of each object is estimated and then those pixels above a set threshold in the profile are masked.

One can enforce more aggressive masking by setting mask_by_boxcar which will mask each object by the boxcar_radius set in ExtractionPar Keywords:

[reduce]
   [[extraction]]
      boxcar_radius = 2.5  # arcsec
   [[skysub]]
      mask_by_boxcar = True

Local Sky Subtraction

Local sky subtraction refines the sky model in the immediate vicinity of each detected object while simultaneously fitting the object’s spatial profile. This coupled approach provides superior sky subtraction near sources and enables optimal spectral extraction.

The local sky subtraction is performed by local_skysub_extract (for multi-slit spectrographs) and ech_local_skysub_extract (for echelle spectrographs).

Core Algorithm

The local sky subtraction operates through an iterative procedure that alternates between profile fitting and sky modeling:

1. Object Grouping

Objects are first organized into extraction groups based on spatial proximity. Two objects are grouped together if their extraction regions (defined by maskwidth) overlap at any spectral position:

groups = sobjs.get_extraction_groups(model_full_slit=model_full_slit)

For each group, a local mask defines the region to be modeled:

localmask = (spat_img > min_spat) & (spat_img < max_spat) & thismask

where min_spat and max_spat encompass all objects in the group plus their mask widths.

2. Iterative Profile Fitting and Sky Modeling

The algorithm performs niter iterations (default 4) of the following steps:

a. Profile Fitting

For each object in the group, fit the spatial profile using extract.fit_profile. On the first iteration, this uses a boxcar extraction to initialize; subsequent iterations use optimal extraction from the previous iteration.

The profile fitting:

Determines the object’s FWHM as a function of wavelength (FWHMFIT)
Updates the object trace position (TRACE_SPAT)
Constructs a 2D profile model (profile_model)

For objects with median S/N < sn_gauss (default 4.0), a Gaussian profile is assumed. Higher S/N objects receive a non-parametric B-spline profile fit.

b. Breakpoint Determination

Generate optimal breakpoints for the B-spline sky model using optimal_bkpts. When bkpts_optimal=True, breakpoints are placed to ensure adequate sampling of the sky signal based on the pixel coordinate distribution within the local region.

c. Joint Sky and Object Fitting

The skyoptimal function performs a simultaneous B-spline fit to both the sky background and object flux contributions:

sky_bmodel, obj_bmodel, outmask_opt = skyoptimal(
    piximg, sciimg, modelivar * skymask, obj_profiles,
    spatial_img=spatial_img, fullbkpt=fullbkpt, sigrej=sigrej_eff)

This fit uses a basis consisting of:

Object spatial profiles for each source
Polynomial terms (Legendre) for spatial sky variations
B-splines in the spectral direction

The model simultaneously solves for the sky and all object contributions within the local region, accounting for potential overlap between sources.

d. Variance Update

When model_noise=True, the inverse variance image is updated using the current sky model to properly account for Poisson noise from the sky:

modelivar = procimg.variance_model(base_var, sky=skyimage, ...)

3. Final Extraction

After the iterative fitting converges, final optimal and boxcar extractions are performed for each object using the refined sky model and object profiles:

Optimal extraction: Uses the fitted spatial profile as weights
Boxcar extraction: Simple aperture sum within boxcar_radius

The `skyoptimal` Function

The core fitting engine skyoptimal performs the coupled sky-object fit:

Basis Construction: Create a combined basis matrix with columns for:
- Each object’s normalized spatial profile (nobj columns)
- Spatial polynomial terms (npoly columns, default 1)

First B-spline Fit: Initial fit with relatively loose rejection:

sset, gpm, yfit1, _, _ = fitting.bspline_profile(
    piximg, data, ivar, profile_basis,
    ingpm=mask, fullbkpt=fullbkpt, upper=sigrej, lower=sigrej)

Chi-squared Rejection: Compute chi-squared for each pixel and apply an additional rejection threshold based on Gaussian statistics:
```
chi2_sigrej = chi2_srt[sigind]  # Threshold from sorted chi2
mask1 = (chi2 < chi2_sigrej)
```

Second B-spline Fit: Refit with tightened mask:

sset, gpm_good, yfit, _, _ = fitting.bspline_profile(
    piximg, data, ivar, profile_basis,
    ingpm=mask1, fullbkpt=fullbkpt, upper=sigrej, lower=sigrej,
    kwargs_reject={'groupbadpix': True, 'maxrej': 1})

Model Separation: Extract the sky and object models from the fit coefficients:

skyset.coeff = sset.coeff[nobj:, :]  # Sky coefficients
for i in range(nobj):
    objset.coeff = sset.coeff[i, :]  # Object i coefficients

[reduce]
   [[skysub]]
      no_local_sky = True

Echelle-Specific Processing

For echelle spectrographs, ech_local_skysub_extract adds additional handling:

Order-by-Order Processing: Each order is processed independently in order of decreasing S/N
FWHM Propagation: The FWHM determined from high-S/N orders is used as prior information for low-S/N orders on the same object
Cross-Object FWHM: Fainter objects use FWHM information from brighter objects (assuming point sources with seeing-limited profiles)

Key Parameters

bspfloat: B-spline breakpoint spacing in pixels. Default is 0.6.
niterint: Number of profile fitting and sky subtraction iterations. Default is 4.
sigrejfloat: Sigma rejection threshold. Default is 3.5.
sn_gaussfloat: S/N threshold below which Gaussian profiles are assumed. Default is 4.0.
force_gaussbool: If True, always use Gaussian profiles regardless of S/N.
model_full_slitbool: If True, model the entire slit width rather than just regions near objects. Recommended for echelle spectra with narrow slits.
bkpts_optimalbool: If True, use optimal breakpoint spacing. If False, use uniform spacing.
model_noisebool: If True, iteratively update the variance model. Should be False for A-B difference imaging where sky residuals are being fit.
no_local_skybool: If True, skip local sky fitting but still perform profile fitting and optimal extraction. Useful for extended emission lines.

Interactively defining the sky regions

PypeIt has an automatic algorithm (described above) to define the sky regions, but this may not work in your specific science case. There are several ways to define the sky regions. The first option is to define the locations on the slits where there is sky in your PypeIt Reduction File. The command is a comma separated list of regions that represent the locations on the slit (0 is the left edge, 100 is the right edge):

[reduce]
  [[skysub]]
       user_regions = :20,65:

where in the example above, the sky regions are defined as all pixels in all slices that are in the leftmost 20 percent of the slit (i.e. :20), and the rightmost 35 percent of the slit (65:). You can specify as many regions as you like. For example, 45:55 would indicate that the innermost 10 percent of pixels contains sky.

An alternative approach is to set the sky regions interactively. This is the preferred approach if you want to set different sky regions for every slit. Remember, you really should assign some sky regions in every slit, otherwise the relative spectral sensitivity correction will not work. To interactively define the sky regions, you must first run through the reduction once, and then use the following command:

pypeit_skysub_regions spec2d_file.fits

You will need to manually define the sky regions for each spec2d file. You will see a GUI where you can click and drag regions on each slit to define the sky regions. Hover the mouse over the window and press the ? key. This will print a list of options in the terminal window, so that you know how to operate the GUI. A left (right) mouse button click and drag will add (remove) pixels to (from) the sky regions mask. Once you have defined some regions, the red shaded regions represent the sky pixels. If you want to set the sky regions for multiple slits, use the “Assign sky regions to all slits” bar on the right hand side of the GUI. The gray region represents the slit, and the black regions represent outside the slit. You need to click and drag only on the gray regions, or you can click and drag from the gray to the black regions (i.e. you must click and drag within this small window for it to work).

Alternatively, you can click the “Enter regions” button, which will request input from the command line. You should now enter the regions in the same format as above for the user_regions.

If you’re happy with the sky regions, press the “Continue (and save changes)” button. If you do not wish to save the sky regions, press the “Continue (and don’t save changes)” button. The menu bar at the top of the screen will prompt you if you wish to save these sky regions (click on either YES or NO). If you chose to save the regions file, the regions will be saved in your Calibrations/ folder, with a prefix SkyRegions. A given SkyRegions file is linked to a science frame based on the name of the SkyRegions file.

Once you have defined all of the sky regions manually, you will need to explicitly tell PypeIt to use the manually defined sky regions file by adding the following lines to your PypeIt Reduction File:

[reduce]
   [[skysub]]
      user_regions = user

and then re-run the reduction.