Wavelength Calibration API#

Arc-lamp analysis, line detection via wavelets, cross-correlation with reference catalogues, polynomial fitting, and spectral rebinning (scrunching). See the arc calibration theory for background.

Reduce Arc Module

This module implements the reduction of arc frames to produce wavelength calibrated spectra.

kspecdr.extract.reduce_arc.reduce_arc(args: Dict[str, Any], get_diagnostic: bool | None = False, diagnostic_dir: Path | None = None) → None[source]#: Reduces a raw arc file to produce im(age), ex(tracted) and red(uced) arc files.

kspecdr.extract.reduce_arc.reduce_arcs(args_list: List[Dict[str, Any]], get_diagnostic: bool | None = False, diagnostic_dir: Path | None = None) → None[source]#

Reduces multiple arc frames together by creating a combined landmark register and fitting a global wavelength solution.

Preprocesses all frames (make_im, make_ex).
Identifies a common reference fiber across all frames.
Extracts landmarks from each frame using the common reference.
Merges landmarks and templates to form a master dataset.
Fits a global wavelength solution.
Applies the global solution to all frames using the combined landmarks.

Arc file I/O operations.

kspecdr.wavecal.arc_io.read_arc_file(nx: int, xvec: ndarray, lamp: str, max_entries: int = 20000, arc_dir: str | None = None) → tuple[ndarray, ndarray, list, int][source]#

Reads an ASCII text ‘arc’ file containing wavelengths and intensities.

Parameters:

nx (int) – Size of xvec
xvec (np.ndarray) – Wavelength axis vector (predicted)
lamp (str) – Lamp name
max_entries (int) – Max entries to read
arc_dir (str) – Directory containing arc files

Returns:

wlist (np.ndarray) – Wavelengths
ilist (np.ndarray) – Intensities
labels (list) – Labels
nlist (int) – Number of entries

Main calibration routine.

kspecdr.wavecal.calibrate.apply_calibration_model(coeffs: ndarray, npix: int, nfib: int, goodfib: ndarray, ref_fib: int, lmr: ndarray, nlm: int) → ndarray[source]#: Propagates the master calibration to all fibers using landmark shifts. Returns pixcal_dp (NPIX+1, NFIB).

kspecdr.wavecal.calibrate.calibrate_spectral_axes(npix: int, nfib: int, spectra: ndarray, variance: ndarray, pred_axis: ndarray, goodfib: ndarray, lamb_tab: ndarray, flux_tab: ndarray, size_tab: int, maxshift: int, diagnostic: bool | None = False, diagnostic_dir: Path | None = None, use_blends: bool = False, poly_order: int = 3) → tuple[ndarray, int, dict][source]#

Calibrate the pixels of extracted arclamp spectra.

Parameters:

use_blends (bool, optional) – If True, blended lines (lines closer than 3-sigma to a neighbour) are included in the calibration rather than being excluded. Downstream quality checks (cross-correlation score, Gaussian fit quality) still filter out lines that cannot be reliably centroided, so only lines that can be individually measured will contribute to the solution. This is useful in spectral regions where lines are densely packed. Default is False (blended lines are excluded as before).

Returns:

pixcal_dp (np.ndarray) – Calibrated pixels (NPIX+1, NFIB)
status (int) – Status code (0 = OK)
resolution_info (dict) – Spectral resolution measurements from arc line fits. Empty dict if calibration fails. See compute_resolution_stats.

kspecdr.wavecal.calibrate.compute_resolution_stats(x_pts: ndarray, y_pts: ndarray, sigma_pix: ndarray, coeffs: ndarray, outliers: ndarray, fwhm_poly_order: int = 1) → dict[source]#

Compute spectral resolution statistics from per-line Gaussian fits.

Uses the wavelength solution polynomial to convert per-line sigma (pixels) to FWHM (Angstrom) via the local dispersion at each line position.

Parameters:

x_pts (pixel positions of matched lines)
y_pts (wavelengths of matched lines)
sigma_pix (per-line Gaussian sigma in pixels)
coeffs (wavelength solution polynomial coefficients (np.polyval order))
outliers (boolean mask from the wavelength fit (True = outlier))
fwhm_poly_order (order of the FWHM(lambda) polynomial fit)

Returns:

Dictionary with keys: sigma_pix_median, fwhm_angstrom_median, resolving_power_median, fwhm_poly_coeffs, and per_line.

Return type:

dict

kspecdr.wavecal.calibrate.extract_template_spectrum(spectra: ndarray, nfib: int, npix: int, goodfib: ndarray, ref_fib: int, cen_axis: ndarray, diagnostic: bool | None = False, diagnostic_dir: Path | None = None) → tuple[ndarray, ndarray, ndarray, float, int][source]#

Extracts a high S/N template spectrum by aligning and stacking fibers.

Returns:: 1D array template_mask: 1D boolean array (True = masked/bad) lmr: Landmark register array (shifts) sigma_inpix: Estimated sigma in pixels nlm: Number of landmarks found
Return type:: template_spectra

kspecdr.wavecal.calibrate.find_arc_line_matches(template_spectra: ndarray, template_mask: ndarray, sigma_inpix: float, cen_axis: ndarray, npix: int, muv: ndarray, av: ndarray, mask: ndarray, maxshift: int, diagnostic: bool | None = False, diagnostic_dir: Path | None = None) → tuple[ndarray, ndarray, ndarray, ndarray][source]#

Identifies arc lines in the template spectrum.

Returns:: Measured pixel positions valid_waves: True wavelengths valid_sigmas: Per-line Gaussian sigma in pixels from the fit final_mask: Boolean mask of lamp lines (True = bad/unused)
Return type:: valid_pixels

kspecdr.wavecal.calibrate.find_reference_fiber(nfib: int, goodfib: ndarray) → int[source]#: Finds a suitable reference fiber (middlemost good fiber).

kspecdr.wavecal.calibrate.fit_calibration_model(x_pts: ndarray, y_pts: ndarray, poly_order: int = 3) → tuple[ndarray, ndarray, ndarray][source]#

Fits a robust polynomial to the points.

Returns:: Polynomial coefficients residuals: Residuals of the fit outliers: Boolean mask of outliers
Return type:: coeffs

kspecdr.wavecal.calibrate.refine_peak_gaussian_fast(spectrum: ndarray, idx_guess: int, sigma0: float, hw: int, *, max_iter: int = 6, clip_sigma: Tuple[float, float] = (0.3, 8.0), min_amp_snr: float = 2.0) → Tuple[float, float, bool][source]#

Fast Gaussian(+linear background) peak refinement around idx_guess.

Model:: y = A * exp(-(x-x0)^2/(2*s^2)) + B + C*(x-xmean)

Returns:: x0_refined (float), sigma_pix (float), ok (bool)

Cross-correlation analysis for wavelength calibration.

kspecdr.wavecal.crosscorr.cross_corr_greedy_quad_path_search(crs_cgm: ndarray, nrows: int, npix: int) → ndarray[source]#: Greedy quadratic path search through cross-correlogram.

kspecdr.wavecal.crosscorr.crosscorr_analysis(template_spectra: ndarray, template_mask: ndarray, npix: int, muv: ndarray, av: ndarray, mask: ndarray, m: int, sigma_inpix: float, cen_axis: ndarray, maxshift: int, diagnostic: bool = False) → ndarray[source]#: Main cross-correlation analysis routine. Returns: fshiftv (fractional shifts)

kspecdr.wavecal.crosscorr.determine_saturated_lines(template_mask: ndarray, cen_axis: ndarray, npix: int, sigma_inpix: float, muv: ndarray, av: ndarray, mask: ndarray, m: int, maxshift: int) → ndarray[source]#: Update mask for saturated lines.

kspecdr.wavecal.crosscorr.gen_cross_corr_gram(s1: ndarray, s2: ndarray, n: int, hw: int, max_shift: int) → ndarray[source]#: Generate cross-correlogram. s1: Template s2: Model

kspecdr.wavecal.crosscorr.generate_spectra_model(muv: ndarray, av: ndarray, mask: ndarray, m: int, sig: float, t: ndarray, n: int) → ndarray[source]#: Generate model spectra from line list.

Landmark registration and signal synchronization.

kspecdr.wavecal.landmarks.landmark_register(spectra: ndarray, npix: int, nfib: int, maskv: ndarray, ref_fib: int, scale: float, ztol: float, diagnostic: bool = False) → tuple[ndarray, int][source]#

Parameters:

spectra (np.ndarray) – Input spectra (npix, nfib) - Fortran uses (NPIX, NFIB)
npix (int) – Number of pixels
nfib (int) – Number of fibers
maskv (np.ndarray) – Mask vector (True if masked/bad)
ref_fib (int) – Reference fiber index
scale (float) – Wavelet scale
ztol (float) – Zero tolerance for peak finding
diagnostic (bool) – Whether to print diagnostic info

Returns:

lmr (np.ndarray) – Landmark Register Array (nfib, nlm) Note: Fortran defines LMR(NFIB, NPIX) but fills it sparsely? The python return should be (nfib, nlm) where nlm is the number of found landmarks. Wait, Fortran: REAL, INTENT(OUT) :: LMR(NFIB,NPIX). And: LMR(FIBNO,USEIDX)=TRACKA(I,SEQIDX). So it returns a 2D array where the 2nd dimension is the landmark index (up to NPIX max). We will return (nfib, nlm) sized array.
nlm (int) – Number of landmarks

kspecdr.wavecal.landmarks.robust_polyfit(x, y, order)[source]#: Robust polynomial fitting using RANSAC.

kspecdr.wavecal.landmarks.synchronise_calibration_last(cal_axis: ndarray, npix: int, nfib: int, maskv: ndarray, ref_fib: int, lmr: ndarray, nlm: int) → ndarray[source]#

Synchronise calibration from ref fibre to others. Iterates outwards from ref_fib to propagate calibration on failure.

cal_axis: Calibration of reference fibre (wavelengths).

kspecdr.wavecal.landmarks.synchronise_signals(spectra: ndarray, npix: int, nfib: int, maskv: ndarray, ref_fib: int, lmr: ndarray, nlm: int) → ndarray[source]#: Rebin spectra to align landmarks. Iterates outwards from ref_fib to propagate calibration on failure.

kspecdr.wavecal.landmarks.windsor_istats(ivec: ndarray, n: int, cutoff_percent: float) → tuple[float, float, bool][source]#

Perform Windsor Statistics analysis on an integer vector. Returns (mean, sd, check_flag).

Calculates the mean and standard deviation of the data, but first truncates the data to the range [cutoff_percent, 100-cutoff_percent] (actually just 100-cutoff_percent logic here to match Fortran description approx).

Actually, standard Windsor statistics usually involves replacing tails. The Fortran code likely does a robust mean/sd calculation assuming X% of data is good. Given “assumption that 75% of the data is good” in comments.

Implementation based on typical Windsorized Mean/SD or Truncated Mean/SD logic. Here we implement a simplified version consistent with common robust stats: Sort data, take the middle cutoff_percent (e.g. 75%), calculate mean/std of that.

Returns:

mean (float)
sd (float)
check_flag (bool) – True if quality seems okay (sd < mean + some tolerance?), or just always True unless empty? Fortran code says: “Output warning if sanity check fails.”

kspecdr.wavecal.scrunch.rebin_data(spectra, variance, input_wave, output_wave)[source]#

Rebins spectra and variance from input_wave to output_wave using linear interpolation.

Parameters:

spectra (np.ndarray) – Input spectra (NSPEC, NFIB)
variance (np.ndarray) – Input variance (NSPEC, NFIB)
input_wave (np.ndarray) – Input wavelength axis. Can be 1D (NSPEC,) or 2D (NSPEC, NFIB).
output_wave (np.ndarray) – Output wavelength axis. Can be 1D (NSPEC_OUT,) or 2D (NSPEC_OUT, NFIB).

Returns:

(out_spec, out_var) - Rebinned spectra and variance

Return type:

tuple

kspecdr.wavecal.scrunch.scrunch_from_arc_id(obj_filename, arc_filename, args, reverse=False)[source]#

Scrunches (rebins) the object file using wavelength solution from the arc file.

If reverse=True, performs the reverse operation: Assumes the object file is currently on the linear grid defined by the arc, and rebins it back to the original pixel grid of the arc.

Parameters:

obj_filename (str) – Path to the object FITS file to scrunch (modified in place)
arc_filename (str) – Path to the arc FITS file containing WAVELA extension.
args (dict) – Arguments dictionary (unused in logic but kept for interface compatibility).
reverse (bool) – If True, un-scrunch (Linear -> Pixel). Default False (Pixel -> Linear).

kspecdr.wavecal.scrunch.scrunch_open_arc(args)[source]#

Returns the filename of the arc file from args. Corresponds to SCRUNCH_OPEN_ARC in Fortran.

Parameters:: args (dict) – Arguments dictionary containing ‘WAVEL_FILENAME’.
Returns:: The filename of the arc file, or None if not found/empty.
Return type:: str or None

Wavelet functions for peak detection and signal analysis. Implements Mexican Hat, Haar, and N2BSpline wavelets and convolution routines.

kspecdr.wavecal.wavelets.analyse_arc_signal(arc_sig: ndarray) → tuple[float, float, float, float, float][source]#

Analyse arc signal to estimate noise and PSF parameters.

Return type:: mn_noise, sd_noise, al_sigma, ares, ztol

kspecdr.wavecal.wavelets.calc_medmad(data: ndarray) → tuple[float, float][source]#: Calculate median and MAD of data.

kspecdr.wavecal.wavelets.find_resonant_peaks2(signal: ndarray, t: ndarray, ztol: float) → ndarray[source]#: Find peaks using zero crossings logic (WAVELET_FIND_RES_PEAKS2). Returns peak locations (interpolated float positions).

kspecdr.wavecal.wavelets.mexican_hat_wavelet(t: ndarray, sigma: float = 1.0) → ndarray[source]#: Calculate the Mexican Hat wavelet at points t. Formula: A * (1 - z^2) * exp(-0.5 * z^2) where z = t / sigma, A = 1 / (sqrt(2*pi) * sigma^3)

kspecdr.wavecal.wavelets.wavelet_convolution(signal: ndarray, t: ndarray, scale: float) → ndarray[source]#

Perform continuous wavelet transform convolution using Mexican Hat wavelet.

Parameters:

signal (np.ndarray) – Input signal.
t (np.ndarray) – Time/position axis for the signal. Assumed to be regularly spaced for convolution.
scale (float) – Wavelet scale parameter (a).

Returns:

Convolved signal (wavelet coefficients).

Return type:

np.ndarray

kspecdr.wavecal.wavelets.wavelet_find_res_peaks_ztol(signal: ndarray, t: ndarray, ztol: float) → ndarray[source]#: Find resonant peaks in signal above zero tolerance. Returns indices of peaks.

kspecdr.wavecal.wavelets.wavelet_find_zero_crossings2(signal: ndarray, t: ndarray, peaks: ndarray, ztol: float) → tuple[ndarray, ndarray][source]#: Find LHS and RHS zero crossings for each peak.