Wavelength Calibration¶

In spectroscopy, the raw data from a detector typically records the intensity of light at different pixel positions. For scientific analysis, we need to know the physical wavelength (or frequency, or energy) corresponding to each pixel. Wavelength calibration is the process of determining this mapping, creating a model that converts pixel coordinates to wavelength values.

This is often achieved by observing a calibration source with well-known emission lines (e.g., an arc lamp). By identifying the pixel positions of these lines, we can fit a mathematical model describing the spectrograph’s dispersion.

The specreduce library provides two complementary classes for 1D calibration:

WavelengthCalibration1D: a high-level calibration workflow to detect and match arc lines and fit a dispersion model.
WavelengthSolution1D: a lightweight wavelength-solution container that holds the pixel→wavelength transform with its inverse, a WCS solution, and a method for flux-conserving resampling.

Tools for calibrating 2D spectra will be added in later versions.

1D Wavelength Calibration¶

The WavelengthCalibration1D class encapsulates the data and methods needed to compute a 1D wavelength solution, and the WavelengthSolution1D class provides the tools to use the computed solution. A typical workflow involves:

Initialization: Create an instance of the class with either arc spectra or pre-identified observed line positions, and provide one or more line lists of known wavelengths. Multiple arcs and multiple catalogs are supported.
Line Detection (if using arc spectra): Identify line positions in the arc spectra.
Fitting the Wavelength Solution: Use either direct fitting with known pixel– wavelength pairs (fit_lines()) or automated global optimization (fit_dispersion()).
Refinement and Inspection: Optionally refine the solution and check fit quality using RMS statistics and plots.
Applying the Solution: Access the fitted solution as a WavelengthSolution1D to resample spectra, or attach the solution to a specutils.Spectrum as a WCS object.

Quickstart¶

1. Initialization¶

You instantiate the WavelengthCalibration1D by providing basic information about your setup and data. You must provide either a list of arc spectra (or a single arc spectrum) or a list of known observed line positions:

Using an Arc Spectrum: Provide the arc spectrum as a specutils.Spectrum object via the arc_spectra argument. You also need to provide a line_lists argument, which can be a list of known catalog wavelengths or the name(s) of standard line lists recognized by specreduce (e.g., "HeI"). You can query the available line lists by using the get_available_line_catalogs function.
```
wc = WavelengthCalibration1D(arc_spectra=arc_he, line_lists=["HeI"])
```

Using multiple Arc Spectra:

wc = WavelengthCalibration1D(arc_spectra=[arc_he, arc_hg_ar],
                             line_lists=[["HeI"], ["HgI", "ArI"]])

Using Observed Line Positions: If you have already identified the pixel centroids of lines in your calibration spectrum, you can provide them directly via the obs_lines argument (as a list of NumPy arrays). In this case, you must also provide the detector’s pixel boundaries using pix_bounds (a tuple like (min_pixel, max_pixel)) and a reference pixel, ref_pixel, which serves as the anchor point for the polynomial fit. You also need to provide the line_lists containing the potential matching catalog wavelengths.
```
obs_he = np.array([105.3, 210.8, 455.1, 512.5, 680.2])
pixel_bounds = (0, 1024)

wc = WavelengthCalibration1D(ref_pixel=512,
                             obs_lines=obs_he,
                             line_lists=["HeI"],
                             pix_bounds=pixel_bounds)
```

Using Observed Line Positions From Multiple Arcs:

obs_he = np.array([105.3, 210.8, 455.1, 512.5, 680.2])
obs_hg_ar = np.array([234.2, 534.1, 768.2, 879.6])
pixel_bounds = (0, 1024)

wc = WavelengthCalibration1D(ref_pixel=512,
                             obs_lines=[obs_he, obs_hg_ar],
                             line_lists=[["HeI"], ["HgI", "ArI"]],
                             pix_bounds=pixel_bounds)

2. Finding Observed Lines¶

If you initialized the class with arc_spectra, you need to detect the lines in it. Use the find_lines() method:

wc.find_lines(fwhm=3.5, noise_factor=5)

This populates the observed_line_locations attribute.

print(wc.observed_line_locations)

3. Matching and Fitting the Solution¶

The core of the calibration process is fitting a model that maps pixels to wavelengths.

Global Fitting for Automated Pipelines: If you have observed_lines (either found automatically or provided initially) and catalog_lines (from line_lists), but don’t know the exact pixel-wavelength pairs, you can use fit_dispersion(). This method applies the differential evolution optimization algorithm to find the polynomial coefficients that minimize the distance between observed line positions (transformed to wavelength space) and the nearest catalog lines.

The fit is anchored around the reference pixel (ref_pixel), which defines the centre of the polynomial model. If not explicitly set at initialization, it defaults to the middle of the detector when arc spectra are supplied. You must provide estimated bounds for both the wavelength and the dispersion (dλ/dx) at ref_pixel.

The polynomial degree is set with the degree argument. A higher degree can capture more complex dispersion behaviour but risks overfitting if too few calibration lines are available.

The higher_order_limits argument optionally constrains the absolute values of the higher-order polynomial coefficients, reducing the risk of unrealistic solutions.

The popsize argument controls the population size in the differential evolution search. Larger values generally improve robustness at the expense of longer computation times.
```
ws = wc.fit_dispersion(wavelength_bounds=(7450, 7550),
                      dispersion_bounds=(1.8, 2.2),
                      higher_order_limits=[0.001, 1e-5],
                      degree=4,
                      popsize=30,
                      refine_fit=True)
```
The fitted WavelengthSolution1D is returned by the method and also stored in WavelengthCalibration1D.solution. Setting refine_fit=True automatically performs a least-squares refinement after the global fit finds an initial solution and matches lines.
Fitting Known Pairs for an Interactive Workflow: If you have already established explicit pairs of observed pixel centers and their corresponding known wavelengths, you can use fit_lines() to perform a direct least-squares fit.
```
ws = wc.fit_lines(pixels=[105.3, 512.5, 780.1],
                 wavelengths=[6965.43, 7503.87, 7723.76],
                 degree=2,
                 refine_fit=True)
```
When refine_fit=True is set, the method automatically identifies matching pairs between observed and catalog lines, then performs a least-squares refinement using all matching lines. This goes beyond the subset of lines provided to fit_lines(), resulting in an improved wavelength solution.

4. Inspecting the Fit¶

Several tools help assess the quality of the wavelength solution:

RMS Error: Calculate the root-mean-square error of the fit in wavelength or pixel units using rms().

rms_wave = wc.rms(space='wavelength')
rms_pix = wc.rms(space='pixel')
print(f"Fit RMS (wavelength): {rms_wave}")
print(f"Fit RMS (pixel): {rms_pix}")

Plotting: Visualize the fit and residuals:
- plot_fit(): Shows the observed line positions mapped to the wavelength axis, overlaid with the catalog lines and the fitted solution. Also shows the fit residuals (observed - fitted wavelength) vs. pixel.
- plot_residuals(): Plots residuals vs. pixel or vs. wavelength.
- plot_observed_lines(): Plots the identified observed line positions (in pixels or mapped to wavelengths). Can optionally overlay the arc spectrum.
- plot_catalog_lines(): Plots the catalog line positions (in wavelengths or mapped to pixels).

5. Using the Wavelength Solution¶

The fitted wavelength solution is stored as a WavelengthSolution1D instance that you can use to transform and resample spectra. The fitting methods return the solution, but it is also stored in WavelengthCalibration1D.solution.

Convert Coordinates: Use pix_to_wav() and wav_to_pix() to convert between pixel and wavelength coordinates.
```
pixels = np.array([100, 500, 900])
wavelengths = ws.pix_to_wav(pixels)
print(wavelengths)
```
Get WCS Object: Access the WCS object representing the solution via the gwcs attribute. This is particularly useful for attaching the calibration to a Spectrum object.
Rebin Spectrum: Resample a spectrum onto a new wavelength grid using resample(). The rebinning is flux-conserving, meaning the integrated flux in the output spectrum matches the integrated flux in the input spectrum.
```
resampled_arc = ws.resample(arc_spectrum, nbins=1000)
print(resampled_arc.spectral_axis)
```
By default, resample() produces a spectrum on a linear wavelength grid. Alternatively, you can pass an arbitrary 1D array of wavelength bin edges via bin_edges to define a custom grid.
```
resampled_arc = ws.resample(arc_spectrum, bin_edges=np.geomspace(4500, 7500, 1000))
print(resampled_arc.spectral_axis)
```

Navigation