import logging
import numpy as np
import pandas as pd
from time import time
from pathlib import Path
from astropy.io import fits
from lmfit.models import PolynomialModel
from inspect import signature
from matplotlib import pyplot as plt
import lime
from lime.tools import ProgressBar, join_fits_files, extract_wcs_header, pd_get, au
from lime.rsrc_manager import lineDB
from lime.fitting.lines import LineFitting, signal_to_noise_rola, sigma_corrections, k_gFWHM, velocity_to_wavelength_band, profiles_computation, linear_continuum_computation
from lime.transitions import Line, lines_frame, _REDSHIFT_DICT, multi_origin_lines_frame
from lime.retrieve.line_bands import determine_line_groups, groupify_lines_df
from lime.io import check_file_dataframe, check_file_array_mask, log_to_HDU, results_to_log, load_frame, LiMe_Error, check_fit_conf, lime_cfg
from lime.fitting.redshift import RedshiftFitting
from lime.plotting.plots import spec_continuum_calculation
from scipy import stats
try:
import aspect
aspect_check = True
except ImportError:
aspect_check = False
_logger = logging.getLogger('LiMe')
def review_bands(spec, line, min_line_pixels=3, min_cont_pixels=2, user_cont_source='central', user_err_from_bands=False):
# Check if the line bands are provided
if line.mask is None:
_logger.warning(f"Line {line} was not found on the input bands database. It won't be measured")
return None
# Check if the line is within the w3, w4 limits
limit_blue, limit_red = spec.wave.compressed()[0], spec.wave.compressed()[-1]
if ((line.mask[2] * (1 + spec.redshift)) < limit_blue) or ((line.mask[3] * (1 + spec.redshift)) > limit_red):
_logger.warning(f"Line {line} bands are outside spectrum wavelengh range: w3 < w_min_rest ({line.mask[2]} < {limit_blue}) or"
f" w4 > w_max_rest ({line.mask[3]} > {limit_red})"
f" it won't be measured")
return None
# Check the units are the same
if line.units_wave != spec.units_wave:
_logger.warning(f"The input {line} line units ({line.units_wave}) are different from the input spectrum units ({spec.units_wave})")
# Check if the spectrum does not have the error arr but the user has requested it
if spec.err_flux is None and user_err_from_bands is False:
_logger.warning(f'The observation does not have an error spectrum but the fit command has requested not to use '
f'the adjacent bands to compute the uncertainty. Please set the "user_err_from_bands=True" to perform'
f' a measurement.')
return None
# Compute the line and adjacent continua indeces:
idcsEmis, idcsCont = line.index_bands(spec.wave, spec.redshift)
if user_cont_source == 'fit':
if spec.cont is None:
raise LiMe_Error(f"The continuum has not been fit. Please run 'Spectrum.fit.continuum' before measuring lines"
f"or change the 'user_cont_source'.")
idcsCont = idcsEmis
# Logic for the bands source
cont_from_bands = True if user_cont_source == 'adjacent' else False
# Check if all the flux entries are masked
emis_flux, cont_flux = spec.flux[idcsEmis], spec.flux[idcsCont]
if np.all(emis_flux.mask):
_logger.warning(f"Line {line} flux is fully masked. It won't be measured")
return None
if np.all(cont_flux.mask) and cont_from_bands:
_logger.warning(f"Line {line} adjacent continua flux is fully masked. It won't be measured")
return None
# Check if all the flux entries are zero
if not np.any(emis_flux):
_logger.warning(f"Line {line} flux entries are all 0. It won't be measured")
return None
if not np.any(cont_flux) and cont_from_bands:
_logger.warning(f"Line {line} continuum flux entries are all 0. It won't be measured")
return None
# Check if the line selection is too narrow
if np.sum(~emis_flux.mask) < min_line_pixels:
_logger.warning(f"Line {line} has only {np.sum(~emis_flux.mask)} pixels. It won't be measured")
return None
# Check if the continua selection is too narrow
if (np.sum(~cont_flux.mask) < min_cont_pixels) and cont_from_bands:
_logger.warning(f"Line {line} continuum bands have only {np.sum(~cont_flux.mask)} pixels. It won't be measured")
return None
return idcsEmis, idcsCont
def import_line_kinematics(line, z_cor, log, fit_conf):
# Check if imported kinematics come from blended component
for idx_child, child_label in enumerate(line.list_comps):
# Check for kinem order
parent_label = fit_conf.get(f'{child_label}_kinem')
if (parent_label is not None) and line.group == 'b':
# Tied kinematics in blended profile
if parent_label in line.list_comps:
idx_parent = line.list_comps.index(parent_label)
factor = f'{line.list_comps[idx_child].wavelength / line.list_comps[idx_parent].wavelength:0.8f}'
fit_conf[f'{child_label}_center'] = {'expr': f'{factor}*{parent_label}_center'}
fit_conf[f'{child_label}_sigma'] = {'expr': f'{factor}*{parent_label}_sigma'}
# Import kinematics from previously measured
elif parent_label in log.index:
mu_parent = log.loc[parent_label, ['center', 'center_err']].to_numpy()
sigma_parent = log.loc[parent_label, ['sigma', 'sigma_err']].to_numpy()
wave_ratio = line.list_comps[idx_child].wavelength/log.loc[parent_label, 'wavelength']
center_child_arr = wave_ratio * (mu_parent / z_cor)
sigma_child_arr = wave_ratio * sigma_parent
# Store the value on the dictionary
fit_conf[f'{child_label}_center'] = {'value': center_child_arr[0], 'vary': False}
fit_conf[f'{child_label}_sigma'] = {'value': sigma_child_arr[0], 'vary': False}
# Error for the propagation
fit_conf[f'{child_label}_center_err'] = center_child_arr[1]
fit_conf[f'{child_label}_sigma_err'] = sigma_child_arr[1]
# Line has not been measured before found
else:
_logger.info(f'\n{parent_label} has not been found on the input lines frame for {child_label} kinematics export.'
f'\n - Please check you are using the right line label and that the line has been measured prior '
f'prior to the current fitting.')
return
def check_cube_bands(input_bands, mask_list, fit_cfg):
if input_bands is None:
# Recover the mask_configuration as a list
for mask_name in mask_list:
mask_fit_cfg = fit_cfg.get(f'{mask_name}_line_fitting')
missing_mask = False
if mask_fit_cfg is not None:
if mask_fit_cfg.get('bands') is None:
missing_mask = True
else:
missing_mask = True
if missing_mask:
error_message = 'No input "bands" provided. In this case you need to include the \n' \
f'you need to specify an "bands=log_file_address" entry the ' \
f'"[{mask_name}_file]" of your fitting configuration file'
raise LiMe_Error(error_message)
return
def recover_level_conf(fit_cfg, mask_key, default_key):
default_cfg = fit_cfg.get(f'{default_key}_line_fitting') if default_key is not None else None
mask_cfg = fit_cfg.get(f'{mask_key}_line_fitting') if mask_key is not None else None
# Case there are not leveled entries
if (default_cfg is None) and (mask_cfg is None):
output_conf = fit_cfg
# Proceed to update the levels
else:
# Default configuration
default_conf = {} if default_cfg is None else default_cfg
default_detect = default_conf.get('line_detection')
# Mask conf
mask_conf = {} if mask_cfg is None else mask_cfg
mask_detect = mask_conf.get('line_detection')
# Update the levels
output_conf = {**default_conf, **mask_conf}
# If no line detection don't add it
if mask_detect is not None:
output_conf['line_detection'] = mask_detect
elif default_detect is not None:
output_conf['line_detection'] = default_detect
else:
pass
return output_conf
def check_compound_line_exclusion(line, lines_df):
# Confirm the dataframe includes the group of lines
group_label = pd_get(lines_df, line, 'group_label', transform='none', nan_to_none=True)
# Confirm if the line is in the group of lines
if group_label is not None:
comp_list = group_label.split('+')
measure_check = False if line in comp_list else True
else:
measure_check = True
return measure_check
def continuum_model_fit(x_array, y_array, idcs, degree):
poly3Mod = PolynomialModel(prefix=f'poly_{degree}', degree=degree)
poly3Params = poly3Mod.guess(y_array[idcs], x=x_array[idcs])
try:
poly3Out = poly3Mod.fit(y_array[idcs], poly3Params, x=x_array[idcs])
cont_fit = poly3Out.eval(x=x_array)
except TypeError:
_logger.warning(f'- The continuum fitting polynomial has more degrees ({degree}) than data points')
cont_fit = np.full(x_array.size, np.nan)
return cont_fit
def res_power_approx(wavelength_arr):
"""
Estimate the spectral resolving power R = λ / Δλ approximation for a wavelength array.
The dispersion per pixel (Δλ/pixel) is computed from the finite differences
of the wavelength array. The resolution element is assumed to be Nyquist-sampled
by 2 pixels, so the FWHM resolution element is 2 * (Δλ/pixel), giving:
R ≈ λ / (2 * Δλ_pixel)
Note: This is an approximation. The true R depends on the slit width,
detector sampling, and optical quality of the spectrograph. For precise
instrumental broadening estimates, an empirical LSF from arc/sky lines
is preferred.
Parameters
----------
wavelength_arr : np.ndarray
1D array of wavelengths, assumed to be in a consistent unit (e.g. Å).
Must be monotonically increasing and uniformly or smoothly sampled.
Returns
-------
res_power : np.ndarray
1D array of resolving power R at each pixel, same shape as wavelength_arr.
Dimensionless.
Notes
-----
- The last pixel is extrapolated by repeating the second-to-last dispersion
value, since np.ediff1d produces N-1 differences for an N-element array.
- If the wavelength array has non-uniform sampling (e.g. from a non-linear
dispersion solution), R will vary across the array accordingly.
- Assumes 2 pixels per resolution element (Nyquist sampling). If your
spectrograph samples the LSF with a different number of pixels, replace
the factor of 2 with the appropriate value.
Examples
--------
>>> wave = np.linspace(4000, 7000, 3000) # 1 Å/pixel
>>> R = res_power_approx(wave)
>>> print(R[0]) # expect ~2000 at 4000 Å with 1 Å/pixel dispersion
2000.0
"""
delta_lambda = np.ediff1d(wavelength_arr, to_end=0)
delta_lambda[-1] = delta_lambda[-2]
return wavelength_arr / (2 * delta_lambda)
def spectrum_resampling(disp_intvl, pixel_width, pixel_number, constant_pixel_width, wave_arr, flux_arr, err_arr, mask_arr):
# Check the input disespersion interval
if disp_intvl is not None:
if disp_intvl[0] < wave_arr[0]:
_logger.warning(f'The input lower dispersion value is below the spectral range: disp_intvl {disp_intvl[0]} < {wave_arr[0]}')
if disp_intvl[-1] > wave_arr[-1]:
_logger.warning(f'The input higher dispersion value is above the spectral range: disp_intvl {disp_intvl[-1]} > {wave_arr[-1]}')
bin_width = np.nanmean(np.diff(disp_intvl))
else:
# Compute the wavelength range based on the pixel width
if pixel_width is not None:
disp_intvl = np.arange(np.round(wave_arr[0]), np.round(wave_arr[-1]), pixel_width)
bin_width = pixel_width
# Compute the wavelength range based on a number of pixels
else:
if pixel_number is not None and pixel_number >= 2:
if not float(pixel_number).is_integer():
_logger.info(f'The input pixel number has been rounded from {pixel_number} to {round(pixel_number)}')
pixel_number = round(pixel_number)
else:
raise ValueError(f'In the number of pixels rebinning the input value must be above 1.')
if constant_pixel_width:
bin_width = np.nanmean(np.diff(wave_arr)) * pixel_number
disp_intvl = np.arange(wave_arr[0], wave_arr[-1], bin_width)
else:
disp_intvl = wave_arr[pixel_number::pixel_number]
bin_width = np.nanmean(np.diff(disp_intvl))
# Compute bin edges from centers — applies to all branches
bin_edges = np.concatenate([[disp_intvl[0] - bin_width / 2], disp_intvl + bin_width / 2])
# Make the binning calculation
flux_binned, edges, binnumber = stats.binned_statistic(wave_arr, flux_arr, statistic='mean', bins=bin_edges)
if err_arr is not None:
err_sum, _, _ = stats.binned_statistic(wave_arr, err_arr ** 2, statistic='sum', bins=bin_edges)
counts, _, _ = stats.binned_statistic(wave_arr, err_arr, statistic='count', bins=bin_edges)
err_binned = np.sqrt(err_sum) / counts
else:
err_binned = None
# nbins = flux_binned.size
# bin_idx = binnumber - 1 # make 0-based
# sum_sq = np.bincount(bin_idx, weights=err_arr ** 2, minlength=nbins)
# counts = np.bincount(bin_idx, minlength=nbins)
# err_binned = np.sqrt(sum_sq) / counts
# # get unique bin numbers #
# uni_bin = np.unique(binnumber)
# err_binned_Svea = []
#
# for binnum in uni_bin[:-1]:
# index_bin = np.where(binnumber == binnum)
# errors_bin = err_arr[index_bin]
# err_bin = np.sqrt(np.sum(errors_bin ** 2)) / (len(errors_bin))
# err_binned_Svea.append(err_bin)
# err_binned_Svea = np.array(err_binned_Svea)
# err_arr = self._spec.err_flux.data
# sum_sq_errors = np.bincount(binnumber, weights=err_arr ** 2)
# bin_counts = np.bincount(binnumber)
# N_bins = flux_binned.size
# sum_sq_errors_filtered = sum_sq_errors[1: N_bins + 1]
# bin_counts_filtered = bin_counts[1: N_bins + 1]
# err_binned = np.sqrt(sum_sq_errors_filtered) / bin_counts_filtered
# # Update the binned wavelength # TODO when do I change this...
# if disp_intvl.size != flux_binned.size:
# disp_intvl = disp_intvl[:-1] + bin_width / 2
return disp_intvl, flux_binned, err_binned
class SpecRetriever:
def __init__(self, spectrum):
self._spec = spectrum
return
[docs]
def lines_frame(self, band_vsigma=70, n_sigma=4, adjust_central_band=True, instrumental_correction=True,
exclude_bands_masked=True, map_band_vsigma=None, grouped_lines=None, automatic_grouping=False,
fit_cfg=None, default_cfg_prefix='default', obj_cfg_prefix=None, update_default=True,
line_list=None, particle_list=None, sig_digits=4, ref_bands=None, vacuum_waves=False,
update_labels=False, update_latex=False, rejected_lines=None, Rayleigh_threshold=2, lines_redshift=None,
map_origin=None, components=None, save_group_label=False):
"""
Return a bands dataframe with the spectral lines within the spectrum wavelength range.
If the user does not provide a ``ref_bands`` This method queries the `LiMe bands database <https://lime-stable.readthedocs.io/en/latest/inputs/n_inputs3_line_bands.html>`_
and returns a :class:`pandas.DataFrame` of transitions visible within the spectrum's observed
wavelength interval, taking into account the spectrum redshift, units, and pixel mask.
The central bands (``w3``–``w4``) are optionally adjusted to match the expected line width,
accounting for the emitting and/or absorbing medium velocity dispersion (``band_vsigma``) and instrumental broadening
(``instrumental_correction``). Lines whose central band falls entirely within a masked
pixel region can be excluded via ``exclude_bands_masked``.
If a fitting configuration (``fit_cfg``) is provided, blended or merged line groups are
resolved and the bands table is updated accordingly.
Parameters
----------
band_vsigma : float, optional
Velocity sigma in km/s used to set the half-width of the central band (``w3``–``w4``).
Default is ``70``.
n_sigma : int, optional
Number of sigma used to compute the band half-width from ``band_vsigma``. Default is ``4``.
adjust_central_band : bool, optional
If ``True`` (default), recompute ``w3`` and ``w4`` from ``band_vsigma``, ``n_sigma``,
and the instrumental broadening.
instrumental_correction : bool, optional
If ``True`` (default), include the instrumental broadening (derived from ``res_power``)
when adjusting the central band width.
exclude_bands_masked : bool, optional
If ``True`` (default), remove lines whose central band pixels are entirely masked.
map_band_vsigma : dict, optional
Per-line overrides for ``band_vsigma``, keyed by line label. Lines not present in the
dict use the global ``band_vsigma`` value.
grouped_lines : dict, optional
Explicit line grouping definitions. If ``None``, grouping is read from ``fit_cfg`` if
available.
automatic_grouping : bool, optional
If ``True``, automatically decide the blended or merged line groups which match the observation. Default is ``False``.
fit_cfg : dict or str or pathlib.Path, optional
Fitting configuration. Can be a dictionary or a path to a configuration file. When
provided, grouped lines and rejected lines are read from this configuration unless
explicitly overridden.
default_cfg_prefix : str, optional
Prefix for default parameter entries in ``fit_cfg``. Default is ``"default"``.
obj_cfg_prefix : str, optional
Prefix for object-specific parameter entries in ``fit_cfg``. Default is ``None``.
update_default : bool, optional
If ``True`` (default), object-specific configuration entries override default entries.
line_list : list or numpy.ndarray, optional
Restrict the output to these line labels. Must follow
`LiMe notation <https://lime-stable.readthedocs.io/en/latest/inputs/n_inputs2_line_labels.html>`_.
particle_list : list or numpy.ndarray, optional
Restrict the output to transitions from these ionic species (e.g. ``["H1", "O3"]``).
sig_digits : int, optional
Number of decimal figures in the line labels. Default is ``4``.
ref_bands : pandas.DataFrame, str, or pathlib.Path, optional
Alternative reference bands database. Defaults to the internal LiMe database.
vacuum_waves : bool, optional
If ``True``, convert wavelengths and band limits to vacuum values. Default is ``False``.
update_labels : bool, optional
If ``True``, recompute line labels from the transition data. Default is ``False``.
update_latex : bool, optional
If ``True``, recompute the ``latex_label`` column. Default is ``False``.
rejected_lines : list, optional
Line labels to exclude from the output. If ``None``, falls back to the value in
``fit_cfg`` if present.
Rayleigh_threshold : float, optional
Minimum wavelength separation (in units of ``band_vsigma``) below which two lines
are considered blended via Rayleigh's criterion. Default is ``2``.
lines_redshift : float, optional
Redshift applied to the transition wavelengths when adjusting the central band, if
no per-line ``z_line`` column is present in the bands table. Falls back to the
spectrum redshift if ``None``.
map_origin : dict, optional
Mapping of origin labels to redshifts for multi-origin line queries.
components : list, optional
List of spectral shape components (e.g. ``["emission", "absorption"]``) used to
filter lines by the predicted profile type from the feature detection algorithm.
Requires ``spec.infer.pred_arr`` to have been computed beforehand.
save_group_label : bool, optional
If ``True``, store the group label in the bands table for grouped lines. Default is
``False``.
Returns
-------
pandas.DataFrame
Bands dataframe with one row per transition, indexed by line label, containing
wavelength, band limits (``w1``–``w6``), and metadata columns.
Notes
-----
- If ``res_power`` is not set on the spectrum, an approximate resolving power is
computed from the wavelength array when ``instrumental_correction=True`` or
``fit_cfg`` is provided.
- ``components`` filtering requires the aspect package and a prior call to the
component detection algorithm; a :exc:`LiMe_Error` is raised otherwise.
- Per-line redshifts in a ``z_line`` column (added by ``map_origin``) take precedence
over ``lines_redshift`` and the spectrum redshift when adjusting central bands.
Examples
--------
Get all lines in the spectrum wavelength range:
>>> bands = spec.retrieve.lines_frame()
Restrict to hydrogen and oxygen transitions with a wider velocity band:
>>> bands = spec.retrieve.lines_frame(particle_list=["H1", "O3"], band_vsigma=120)
Use a fitting configuration to resolve blended lines:
>>> bands = spec.retrieve.lines_frame(fit_cfg="my_cfg.toml")
"""
# Remove the mask from the wavelength array if necessary
wave_intvl = self._spec.wave.compressed()
# Check configuration format
in_cfg = check_fit_conf(fit_cfg, default_cfg_prefix, obj_cfg_prefix, update_default) if fit_cfg else None
# Generate the table of single lines taking into account possible origins
rejected_lines = rejected_lines if rejected_lines is not None else (in_cfg or {}).get('rejected_lines')
bands = multi_origin_lines_frame(map_origin, line_list, self._spec.redshift, wave_intvl=wave_intvl,
particle_list=particle_list, units_wave=self._spec.units_wave,
sig_digits=sig_digits, ref_bands=ref_bands, vacuum_waves=vacuum_waves,
update_labels=update_labels, update_latex=update_latex,
rejected_lines=rejected_lines)
# Compute the resolving power if necessary
if self._spec.res_power is not None:
res_power = self._spec.res_power
else:
res_power = res_power_approx(wave_intvl) if (instrumental_correction or fit_cfg is not None) else None
# Adjust the middle bands to match the line width
if adjust_central_band:
# Input the origin redshift
if 'z_line' in bands.columns:
z_corr = 1 + np.nan_to_num(bands['z_line'].to_numpy(), nan=self._spec.redshift)
# Use the function redshift
elif lines_redshift is not None:
z_corr = 1 + lines_redshift
# Use the object redshift
else:
z_corr = 1 + self._spec.redshift
# Expected transitions in the observed frame # TODO This could be read from the configuration file
lambda_obs = bands.wavelength.to_numpy() * z_corr
# Add correction for the instrumental broadening
if instrumental_correction:
# Indexes for the lines emission peak
idcs = np.searchsorted(wave_intvl, lambda_obs)
# Use the instrumental resolution if available
delta_lambda_inst = lambda_obs / (res_power[idcs] * k_gFWHM)
# Constant velocity width
else:
delta_lambda_inst = 0
# Use unique or specific velocity sigma for the bands
map_band_vsigma = map_band_vsigma if map_band_vsigma else (in_cfg or {}).get('map_band_vsigma')
if map_band_vsigma is not None:
band_vsigma = np.full(lambda_obs.size, band_vsigma)
for idx in bands.index.get_indexer(map_band_vsigma.keys()):
if idx > -1:
band_vsigma[idx] = map_band_vsigma[bands.index[idx]]
# Convert to spectral width
delta_lambda = velocity_to_wavelength_band(n_sigma, band_vsigma, lambda_obs, delta_lambda_inst)
# Add new values to database in the rest frame
bands['w3'] = (lambda_obs - delta_lambda) / z_corr
bands['w4'] = (lambda_obs + delta_lambda) / z_corr
# Remove from the output bands those which have all their pixels masked
if exclude_bands_masked:
idcs_w3_w4 = np.searchsorted(self._spec.wave.data/(1+self._spec.redshift), bands.loc[:, 'w3':'w4'])
idcs_valid = [idx for idx, start, end in zip(bands.index, idcs_w3_w4[:, 0], idcs_w3_w4[:, 1])
if not np.all(self._spec.flux[start:end].mask)]
bands = bands.loc[idcs_valid]
# Combine the blended/merged lines in the bands table
if in_cfg:
grouped_lines = grouped_lines if grouped_lines else (in_cfg or {}).get('grouped_lines')
groups_dict = determine_line_groups(self._spec, bands, in_cfg, grouped_lines, automatic_grouping, n_sigma,
Rayleigh_threshold)
groupify_lines_df(bands, in_cfg, groups_dict, self._spec, save_group_label)
# Filter the table to match the line detections
if components:
if aspect_check:
if self._spec.infer.pred_arr is not None:
# Create masks for all intervals
starts = bands.w3.to_numpy()[:, None] * (1 + self._spec.redshift)
ends = bands.w4.to_numpy()[:, None] * (1 + self._spec.redshift)
# Check if x values fall within each interval
in_intervals = (self._spec.wave.data >= starts) & (self._spec.wave.data < ends)
# Check where y equals the target category
shape_indexes = [aspect.cfg['shape_number'][comp] for comp in components]
is_target_category = np.isin(self._spec.infer.pred_arr, shape_indexes)
# Combine the masks to count target_category occurrences in each interval
counts = np.sum(in_intervals & is_target_category, axis=1)
# Check which intervals satisfy the minimum count condition
idcs = counts >= 4
bands = bands.loc[idcs]
else:
raise LiMe_Error(f'No components data. Please run the components detection algorithm')
else:
raise LiMe_Error(f'Aspect is not installed')
return bands
def spectrum(self, redshift=None, norm_flux=None, crop_waves=None, crop_flux=None, pixel_mask=None, mask_intvls=None,
obj_redshift=False):
# Extract the spectrum data
mask_arr = self._spec.flux.mask
wave_arr = self._spec.wave.data
flux_arr = self._spec.flux.data if self._spec.norm_flux is None else self._spec.flux.data * self._spec.norm_flux
if self._spec.err_flux is not None:
err_arr = self._spec.err_flux.data if self._spec.norm_flux is None else self._spec.err_flux.data * self._spec.norm_flux
else:
err_arr = None
# Use the original pixel mask if none is provided else combine
pixel_mask = mask_arr if pixel_mask is None else mask_arr | pixel_mask
# Add the user masked regions
if mask_intvls is not None:
# LiMe lines frame
if isinstance(mask_intvls, pd.DataFrame):
if {'w3', 'w4'}.issubset(mask_intvls):
mask_intvls = mask_intvls.loc[:, ['w3','w4']].to_numpy()
else:
raise LiMe_Error(f'In order to use a pandas dataframe as "mask_intvls" it must have the a "w3" and'
f' "w4" column to specify the lines location. Alternatively, you may use a matrix '
f'array to define the masked intervals')
# Array of intervals (check it has the rigth dimensions)
else:
# Check the input has the right format
mask_intvls = np.asarray(mask_intvls, dtype=float)
if mask_intvls.ndim != 2 or mask_intvls.shape[1] != 2:
raise ValueError('The argument "exclude_intvls" must be a list of (low, high) pairs')
# Add redshift correction
if obj_redshift:
mask_intvls = mask_intvls * (1 + self._spec.redshift) if redshift is None else mask_intvls * (1 + redshift)
# Loop through the wavelength intervals and add them to the mask
# exclude_intvls = exclude_intvls * z_corr
i0 = np.searchsorted(wave_arr, mask_intvls[:, 0], side="right")
i1 = np.searchsorted(wave_arr, mask_intvls[:, 1], side="left")
for start, stop in zip(i0, i1):
pixel_mask[start:stop] = True
# Recreate the spectrum
out_spec = lime.Spectrum(input_wave=wave_arr,
input_flux=flux_arr,
input_err=err_arr,
redshift=self._spec.redshift if redshift is None else redshift,
res_power=self._spec.res_power,
units_wave=self._spec.units_wave,
units_flux=self._spec.units_flux,
norm_flux=norm_flux,
crop_waves=crop_waves,
crop_flux=crop_flux,
pixel_mask=pixel_mask)
return out_spec
def rebinned(self, disp_intvl=None, pixel_width=None, pixel_number=None, const_pixel_width=True, rest_frame=False,
return_spectrum=False):
"""
Rebin the spectrum onto a new dispersion grid and return binned flux (and errors).
Exactly **one** of ``disp_intvl``, ``pixel_width``, or ``pixel_number`` must be
provided to define the target bins:
- ``disp_intvl``: explicit bin **edges** (1D array, strictly increasing).
- ``pixel_width``: constant bin width (same units as wavelength); bins span the
full observed range ``[wave[0], wave[-1])``.
- ``pixel_number``: group pixels by this integer factor. If ``constant_pixel_width``
is ``True`` (default), uses an average native dispersion to build uniform-width
bins; otherwise, takes every ``pixel_number``-th native wavelength as an edge
(non-uniform bin widths).
The user can choose to perform the rebining in the observation restframe (this will
set redshift=0 on the return spectrum).
The binned flux is the mean flux per bin. If an uncertainty array is available,
per-bin errors are combined as the square-root of the sum of variances divided
by the number of contributing samples.
Parameters
----------
disp_intvl : array-like, optional
Monotonic array of bin **edges** for the target dispersion grid. If provided,
it overrides the other arguments. Values outside the native spectral range
will trigger a warning.
pixel_width : float, optional
Desired constant bin width (in wavelength units). Used to generate uniform
bin edges covering the full native wavelength range.
pixel_number : int, optional
Number of native pixels to aggregate per bin. Must be ``>= 2``. If a
non-integer is given, it is rounded and a message is logged.
constant_pixel_width : bool, optional
When rebinned by ``pixel_number``:
- If ``True`` (default), compute a uniform bin width as the average native
spacing times ``pixel_number``.
- If ``False``, construct edges by taking every ``pixel_number``-th native
wavelength (non-uniform bins).
rest_frame: bool, optional
Perform the rebinning on the observation rest frame.
return_spectrum : bool, optional
If ``False`` (default), return the binned wavelength, flux and flux uncertainty arrays.
If ``True``, return a new :class:`lime.Spectrum` instance containing the binned data.
Returns
-------
disp_centers : numpy.ndarray
Binned wavelength array (same units as disp_intvl if provided).
flux_binned : numpy.ndarray
Binned flux array.
err_binned : numpy.ndarray or None
Binned flux uncertainty array.
spectrum : lime.Spectrum, optional
Returned only if ``return_spectrum=True``. A new spectrum object with the binned data.
Raises
------
ValueError
If none or more than one of ``disp_intvl``, ``pixel_width``, ``pixel_number``
are provided, or if ``pixel_number < 2``.
Notes
-----
- ``disp_intvl`` is interpreted as **edges** (length = ``nbins + 1``); the
returned ``disp_centers`` are computed as edge midpoints.
- Binning uses ``scipy.stats.binned_statistic(..., statistic='mean')`` for flux.
- If the spectrum contains masked or NaN values, ensure they are appropriately
handled upstream; NaNs within a bin will propagate to the mean.
"""
# Check only one approach is necessary
provided = {"dispersion inteval": disp_intvl, "pixel width": pixel_width, "pixel number": pixel_number}
active = [name for name, val in provided.items() if val is not None]
# Enforce only one rebinning array
if len(active) == 0:
raise ValueError("No arguments were provided to compute the new spectral resolution")
if len(active) > 1:
raise ValueError(f"Arguments {active} are mutually exclusive. Please only one.")
# Extract the spectrum data
z_corr = 1 if rest_frame is False else (1 + self._spec.redshift)
mask_arr = self._spec.wave.mask
wave_arr = self._spec.wave.data / z_corr
flux_arr = self._spec.flux.data * z_corr
err_arr = self._spec.err_flux.data * z_corr if self._spec.err_flux is not None else None
# Perform the rebinne
disp_intvl, flux_binned, err_binned = spectrum_resampling(disp_intvl, pixel_width, pixel_number, const_pixel_width,
wave_arr, flux_arr, err_arr, mask_arr)
if np.any(mask_arr):
if mask_arr.size != disp_intvl.size:
_logger.warning('Recalculating input mask')
mask_arr = None
if not return_spectrum:
return disp_intvl, flux_binned, err_binned
else:
return lime.Spectrum(input_wave=disp_intvl,
input_flux=flux_binned * self._spec.norm_flux if self._spec.norm_flux else flux_binned,
input_err=err_binned * self._spec.norm_flux if self._spec.norm_flux else err_binned,
redshift=0 if rest_frame else self._spec.redshift , res_power=self._spec.res_power,
units_wave=self._spec.units_wave, units_flux=self._spec.units_flux,
norm_flux=self._spec.norm_flux)
def normalization(self, return_spectrum=False, crop_waves=None, crop_flux=None, **kwargs):
"""
Normalize the spectrum by its continuum.
If the continuum has not been previously computed (``self._spec.cont is None``),
this method automatically fits it by calling
``self._spec.fit.continuum(**kwargs)`` before performing the normalization.
The normalized flux is defined as ``flux / continuum``. If a flux uncertainty
array is available, the normalized uncertainty is propagated assuming
independent errors in the flux and continuum.
Parameters
----------
return_spectrum : bool, optional
If ``False`` (default), return the normalized flux and uncertainty arrays.
If ``True``, return a new :class:`lime.Spectrum` instance containing the
normalized data.
**kwargs
Additional keyword arguments passed directly to
``lime.Spectrum.fit.continuum`` when the continuum needs to be computed.
Returns
-------
flux_norm : astropy.units.Quantity
Continuum-normalized flux array.
err_norm : astropy.units.Quantity or bool
Uncertainty of the normalized flux. If the original spectrum does not
contain a flux uncertainty array, this returns ``False``.
spectrum : lime.Spectrum, optional
Returned only if ``return_spectrum=True``. A new spectrum object with
dimensionless flux units and ``norm_flux=1``.
Notes
-----
- The uncertainty propagation follows:
``σ(F/C) = |F/C| * sqrt[(σ_F / F)^2 + (σ_C / C)^2]``
where ``F`` is the flux and ``C`` is the continuum.
- The returned spectrum preserves the original wavelength grid, redshift,
spectral resolution, and pixel mask.
- The continuum is computed only once and cached in the parent spectrum.
"""
# Compute the object continuum if not provided
if self._spec.cont is None:
self._spec.fit.continuum(**kwargs)
# Normalize the flux
flux_norm = self._spec.flux / self._spec.cont
# Normalize the flux uncertainty
if self._spec.err_flux is not None:
err_norm = np.abs(flux_norm) * np.sqrt((self._spec.err_flux / self._spec.flux) ** 2 +
(self._spec.cont_std / self._spec.cont) ** 2)
else:
err_norm = False
# Return the normalized flux and uncertainty array
if not return_spectrum:
return flux_norm, err_norm
# Return a LiMe spectrum
else:
return lime.Spectrum(self._spec.wave.data, flux_norm.data, err_norm.data, redshift=self._spec.redshift,
units_wave=self._spec.units_wave, units_flux=au.dimensionless_unscaled, norm_flux=1,
res_power=self._spec.res_power, pixel_mask=flux_norm.mask, crop_waves=crop_waves,
crop_flux=crop_flux)
def upper_line_limit(self, line, bands=None, signal_to_noise=8, err_from_bands=False, continua_sigma=True):
# Use frame if available
bands = self._spec.frame if bands is None else check_file_dataframe(bands)
# Check if line is avaible to calculations
if (bands is None) or (line not in bands.index):
raise LiMe_Error('Upper flux limit requires the line band. Please input the lines frame or include in the '
'spec.lines_frame measurements' if bands is None else f'Line "{line}" not found in the input bands')
line = Line.from_transition(line, data_frame=bands)
idcs_central, idcs_continua = line.index_bands(self._spec.wave, self._spec.redshift)
lambda_central = line.wavelength * (1 + self._spec.redshift)
if err_from_bands:
sigma_cont = np.mean(self._spec.err_flux[idcs_central]) * self._spec.norm_flux
else:
sigma_cont = line.measurements.cont_err
print(line.measurements.intg_flux)
print(line.measurements.profile_flux)
delta_lambda = np.mean(np.diff(self._spec.wave[idcs_central]))
R_line = np.mean(self._spec.res_power[idcs_central])
g_const = np.sqrt(np.pi/(2*np.log(2)))
upper_flux = signal_to_noise * sigma_cont * np.sqrt(g_const * (lambda_central / R_line) * delta_lambda)
print('lambda_central', lambda_central)
print('delta_lambda',delta_lambda)
print('sigma_cont',sigma_cont)
print('R_line',R_line)
print('g_const',g_const)
print('upper_flux',upper_flux)
return upper_flux
class SpecTreatment(LineFitting, RedshiftFitting):
def __init__(self, spectrum):
# Instantiate the dependencies
LineFitting.__init__(self)
# Lime spectrum object with the scientific data
self._spec = spectrum
self.line = None
self._i_line = 0
self._n_lines = 0
self.cov_linear = None
[docs]
def bands(self, label, bands=None, fit_cfg=None, min_method='least_squares', profile=None, shape=None,
cont_source='central', err_from_bands=None, temp=10000.0, default_cfg_prefix='default', obj_cfg_prefix=None,
update_default=True):
"""
Fit a spectral line from defined bands (see :ref:`bands documentation <line-bands-doc>`.).
This method performs a full line measurement and profile fitting. The line is query from the default
line database if no ``bands`` dataframe is provided.
The function will also query the input ``fit_cfg`` dictionary if provided for the configuration settings.
Parameters
----------
label : str or float
Line label in `LiMe notation
<https://lime-stable.readthedocs.io/en/latest/inputs/n_inputs2_line_labels.html>`_,
or a transition wavelength (in the same units as the spectrum).
If a numeric wavelength is given, the corresponding transition is
queried from the ``bands`` table (or falling back to the default database if not provided).
bands : pandas.DataFrame, str, or pathlib.Path, optional
Either:
* A bands DataFrame or file path to one.
* If ``None``, the default LiMe bands database is used.
fit_cfg : dict, str, or pathlib.Path, optional
A dictionary or a path to a .toml file.
min_method : str, optional
Minimization algorithm used by :mod:`lmfit`.
Supported methods are listed in the
`lmfit.minimizer.Minimizer.minimize
<https://lmfit.github.io/lmfit-py/fitting.html#lmfit.minimizer.Minimizer.minimize>`_
documentation. Default is ``"least_squares"``.
profile : str, optional
Profile type for fitting (e.g., ``"g"`` for Gaussian, ``"l"`` Lorentz, ...).
If none is provided, the algorithm will use the default profile from the lines' database.
shape : str, optional
Line shape for the fitted profiles: "emi" for emission or "abs" for absroption.
If none is provided, the algorithm will use the default shape from the lines' database.
cont_source : {'central', 'adjacent', 'fitted'}, optional
Method used to estimate the continuum level for line fitting.
- ``'central'`` — use the edges of the central line band (w₃–w₄).
- ``'adjacent'`` — use the adjacent continuum bands (w₁–w₂ and w₅–w₆).
- ``'fitted'`` — use a previously fitted continuum model from the spectrum.
The default is ``'central'``.
err_from_bands : bool or None, optional
If ``True``, estimate the pixel uncertainty from the continuum bands.
If ``None``, use the spectrum’s ``err_flux`` data or fall back to the
continuum regions if not available (False). The default value is None.
temp : float, optional
Electron temperature in Kelvin used to compute the thermal broadening
correction for the fitted lines. Default is 10,000 K.
default_cfg_prefix : str, optional
Section key prefix for the default configuration in ``fit_cfg``.
Default is ``"default"``.
obj_cfg_prefix : str, optional
Section key prefix for the object-specific configuration in ``fit_cfg``.
update_default : bool, optional
If ``True`` (default), merge parameters from ``obj_cfg_prefix`` into
``default_cfg_prefix``. If ``False``, the object configuration is used falling back to the default if not available.
Returns
-------
None
The results are stored in the spectrum’s internal ``frame`` attribute
and in ``self.line.measurements``.
Notes
-----
- The method performs the following steps:
1. Parse or retrieve the line definition using :meth:`~lime.Line.from_transition`.
2. Select line and continuum regions based on the provided or default ``bands``.
3. Estimate the continuum level and its uncertainty.
4. Compute non-parametric line properties (e.g., integrated flux).
5. Perform optional profile fitting via :mod:`lmfit` using ``min_method``.
6. Apply instrumental and thermal corrections to measured line widths.
7. Recalculate the signal-to-noise ratio and store all results in the spectrum’s log frame.
- Continuum and error estimations are controlled via the
``cont_from_bands`` and ``err_from_bands`` flags.
- Thermal broadening corrections use the provided ``temp`` parameter.
- All results are written to the current spectrum’s measurement log
and accessible through ``Spectrum.frame``.
Examples
--------
Fit a Gaussian emission line using the default configuration:
>>> spec.fit.bands("O3_5007A")
Use a custom bands table and configuration dictionary:
>>> spec.fit.bands("H1_4861A", bands="my_bands.xlsx", fit_cfg=my_fit_cfg)
Change the minimization algorithm and temperature:
>>> spec.fit.bands("O2_3726A", min_method="nelder", temp=12000)
Fit a line providing the central wavelength directly:
>>> spec.fit.bands(5007.0, bands=my_bands_df)
"""
# Reset attributes
self.cov_linear = None
# Make a copy of the fitting configuration
input_conf = check_fit_conf(fit_cfg, default_cfg_prefix, obj_cfg_prefix, update_default)
# User inputs override default behaviour for the pixel error and the continuum calculation
err_from_bands = True if (err_from_bands is None) and (self._spec.err_flux is None) else err_from_bands
# Interpret the input line
if isinstance(label, str):
self.line = Line.from_transition(label, input_conf,
data_frame=lineDB.frame if bands is None else check_file_dataframe(bands, copy_input=False),
def_shape=shape,
def_profile=profile)
else:
self.line = label
# Check the line selection is valid
idcs_selection = review_bands(self._spec, self.line, user_cont_source=cont_source, user_err_from_bands=err_from_bands)
if idcs_selection is not None:
# Unpack the line selections
idcs_line, idcs_continua = idcs_selection
# Compute line continuum and the pixel error
cont_arr, pixel_err_arr = self._cont_level_profile(idcs_line, idcs_continua, cont_source, err_from_bands)
# Non-parametric measurements
self.integrated_properties(self.line, self._spec.wave[idcs_line], self._spec.flux[idcs_line],
pixel_err_arr[idcs_line], cont_arr[idcs_line])
# Import kinematics if requested
import_line_kinematics(self.line, 1 + self._spec.redshift, self._spec.frame, input_conf)
# Profile fitting measurements
idcs_fitting = idcs_selection[0] if cont_source == 'central' else idcs_selection[0] + idcs_selection[1]
self.profile_fitting(self.line,
x_arr=self._spec.wave[idcs_fitting],
y_arr=self._spec.flux[idcs_fitting],
err_arr=pixel_err_arr[idcs_fitting],
cont_arr=cont_arr[idcs_fitting] if cont_source == 'fit' else None,
user_conf=input_conf, fit_method=min_method)
# Instrumental and thermal corrections for the lines
sigma_corrections(self.line, idcs_line, self._spec.wave[idcs_line], self._spec.res_power, temp)
# Recalculate the SNR with the profile parameters
self.line.measurements.snr_line = signal_to_noise_rola(self.line.measurements.amp,
self.line.measurements.cont_err,
self.line.measurements.n_pixels)
# Save the line parameters
results_to_log(self.line, self._spec.frame, self._spec.norm_flux)
return
[docs]
def frame(self, bands, fit_cfg=None, min_method='least_squares', profile=None, shape=None, cont_source='central',
err_from_bands=None, temp=10000.0, line_list=None, default_cfg_prefix='default', obj_cfg_prefix=None,
update_default=True, line_detection=False, plot_fit=False, progress_output='bar'):
"""
Measure multiple spectral lines from a bands dataframe
(see :ref:`bands documentation <line-bands-doc>`).
This method automates the fitting of the lines on the input lines frame.
It iterates through all listed transitions,
performing continuum estimation, line detection (optional), and profile
fitting using the configuration provided in ``fit_cfg``.
Parameters
----------
bands : pandas.DataFrame, str, or pathlib.Path
Bands table defining the line labels and bands limits (w1 to w6) for each line.
fit_cfg : dict, str, or pathlib.Path, optional
Fitting configuration dictionary or a path to a TOML configuration file.
See the `profile fitting documentation
<https://lime-stable.readthedocs.io/en/latest/inputs/n_inputs4_fit_configuration.html>`_.
min_method : str, optional
Minimization algorithm used by :mod:`lmfit`.
See the
`lmfit.minimizer.Minimizer.minimize
<https://lmfit.github.io/lmfit-py/fitting.html#lmfit.minimizer.Minimizer.minimize>`_
documentation for available options. Default is ``"least_squares"``.
profile : str, optional
Profile type for fitting (e.g., ``"g"`` for Gaussian, ``"l"`` for Lorentzian).
Defaults to the line database entry if omitted.
shape : str, optional
Line shape keyword (``"emi"`` for emission or ``"abs"`` for absorption).
Defaults to the line database entry if omitted.
cont_source : {'central', 'adjacent', 'fitted'}, optional
Method used to estimate the continuum level for line fitting.
- ``'central'`` — use the edges of the central line band (w₃–w₄).
- ``'adjacent'`` — use the adjacent continuum bands (w₁–w₂ and w₅–w₆).
- ``'fitted'`` — use a previously fitted continuum model from the spectrum.
The default is ``'central'``.
err_from_bands : bool or None, optional
If ``True``, estimate the pixel uncertainty from the continuum bands.
If ``None``, use the spectrum’s ``err_flux`` data or fall back to the
continuum regions if not available (False). The default value is None.
temp : float, optional
Electron temperature (K) used to compute the thermal broadening correction.
Default is 10,000 K.
line_list : list of str, optional
Subset of line labels from the bands table to process.
If ``None``, all entries in ``bands`` are measured.
default_cfg_prefix : str, optional
Section key prefix for the default configuration in ``fit_cfg``.
Default is ``"default"``.
obj_cfg_prefix : str, optional
Section key prefix for the object-specific configuration in ``fit_cfg``.
update_default : bool, optional
If ``True`` (default), merge parameters from ``obj_cfg_prefix`` into
``default_cfg_prefix``. If ``False``, the object configuration is used falling back to the default if not available.
line_detection : bool, optional
If ``True``, run the continuum fitting and line threshodling to confirme the presence of lines before
measurements. The functions parameters must be specified in the ``fit_cfg`` (e.g., entries under ``"peaks_troughs"``, ``"continuum"``).
Default is ``False``.
plot_fit : bool, optional
If ``True``, display the profile fit after each iteration.
progress_output : {"bar", "counter", None}, optional
Controls progress display in the console.
- ``"bar"`` (default): show a dynamic progress bar.
- ``"counter"``: print current line number and label.
- ``None``: suppress console output.
Returns
-------
None
The resulting measurements are stored in the spectrum’s internal
``frame`` attribute and in ``self.line.measurements``.
Notes
-----
- This method performs the following sequence for each line:
1. Optionally apply continuum and line detection preprocessing steps if enabled via ``line_detection=True`` and the appropicate keys are found at ``fit_cfg``.
2. Retrieve the line list from the ``bands`` table or the default database.
3. Estimate the continuum level and its uncertainty.
4. Perform non-parametric measurements (e.g., flux, EW, FWHM).
5. Run profile fitting using :mod:`lmfit` according to ``min_method``.
6. Apply instrumental and thermal width corrections (via ``Spectrum.res_power`` and ``temp``).
7. Recalculate SNR and store results in the spectrum’s log frame.
- Progress reporting is configurable through ``progress_output``.
- Use ``line_detection=True`` to automatically threshold and select
only detected lines before fitting.
Examples
--------
Measure all lines from a bands file:
>>> spec.fit.frame("my_bands.xlsx", fit_cfg="my_fit_config.toml")
Run the fit with a progress bar:
>>> spec.fit.frame(bands_df, progress_output="bar")
Limit to a subset of lines:
>>> spec.fit.frame(bands_df, line_list=["O3_5007A", "H1_4861A"])
Enable automatic line detection:
>>> spec.fit.frame(bands_df, line_detection=True)
"""
# Check if the lines log is a dataframe or a file address
bands = check_file_dataframe(bands)
if bands is not None:
# Crop the analysis to the target lines
if line_list is not None:
idcs = bands.index.isin(line_list)
bands = bands.loc[idcs]
# Load configuration
input_conf = check_fit_conf(fit_cfg, default_cfg_prefix, obj_cfg_prefix, update_default=update_default,
line_detection=line_detection)
# Line detection if requested
if line_detection:
# Review the configuration entries
cont_fit_conf = input_conf.get('continuum', None)
if cont_fit_conf:
self._spec.fit.continuum(**cont_fit_conf)
else:
_logger.warning(f'No "continuum" entry in input configuration file. No continuum fitting will be applied')
# Perform the line detection
detect_conf = input_conf.get('peaks_troughs', {})
if detect_conf:
bands = self._spec.infer.peaks_troughs(bands, **detect_conf)
else:
_logger.warning(f'No "peaks_troughs" entry in input configuration file. No line thresholding will be applied')
else:
cont_fit_conf, detect_conf = None, None
# Define lines to treat through the lines
label_list = bands.index.to_numpy()
self._n_lines = label_list.size
# Loop through the lines
if self._n_lines > 0:
# On screen progress bar
pbar = ProgressBar(progress_output, f'{self._n_lines} lines')
if progress_output is not None:
print(f'\nLine fitting progress{" (continuum fitting)" if cont_fit_conf is not None else ""}'
f'{" (line detection)" if detect_conf is not None else ""}:')
for self._i_line in np.arange(self._n_lines):
# Ignore line if part of a blended/merge group
line = label_list[self._i_line]
measure_check = check_compound_line_exclusion(line, bands)
if measure_check:
# Progress message
pbar.output_message(self._i_line, self._n_lines, pre_text="", post_text=f'({line})')
# Fit the lines
self.bands(line, bands, input_conf, min_method, profile, shape, cont_source=cont_source,
err_from_bands=err_from_bands, temp=temp, obj_cfg_prefix=None, default_cfg_prefix=None)
if plot_fit:
self._spec.plot.bands()
else:
msg = f'No lines were measured from the input dataframe:\n - line_list: {line_list}\n - line_detection: {line_detection}'
_logger.debug(msg)
else:
_logger.info(f'Not input dataframe. Lines were not measured')
return
[docs]
def continuum(self, degree_list, emis_threshold, abs_threshold=None, smooth_scale=None, exclude_intvls=None,
rest_intvls=False, plot_steps=False, **kwargs):
"""
Fit the spectrum continuum via polynomial clipping.
This routine estimates the continuum by iteratively fitting polynomials and
sigma-clipping outliers above (emission) and below (absorption) a flux threshold.
At each iteration, points outside configurable residual thresholds are excluded
and the polynomial is refitted on the remaining pixels.
The user may optionally provide a list of wavelength intervals to be excluded from
the continuum fitting. By default, these limits are assumed to be defined in the
observed frame.
Parameters
----------
degree_list : list of int
Polynomial degree to use at each iteration.
The number of iterations equals ``len(degree_list)``.
emis_threshold : list of float
Upper (emission-side) clipping factors, in units of number of residual standard
deviation for each iteration. Must have the same length as ``degree_list``.
abs_threshold : list of float, optional
Lower (absorption-side) clipping factors, also in units of number of residual
standard deviation. If ``None``, the values in ``emis_threshold`` are reused
for the lower limit. Must match the length of ``degree_list`` when provided.
smooth_scale : int, optional
Window size (in pixels) for a moving-average smoothing applied to the input
flux before fitting. If ``None``, no smoothing is applied.
exclude_intvls : list of tuple(float, float), optional
List of wavelength intervals (low, high) to exclude from the continuum fitting.
By default, intervals are interpreted in the observed frame.
rest_intvls : bool, optional
If ``True``, the wavelength intervals in ``exclude_intvls`` are assumed to be
defined in the rest frame and are converted to the observed frame using
``λ_obs = λ_rest × (1 + z)`` prior to mask computation.
plot_steps : bool, optional
If ``True``, display a diagnostic plot after each iteration showing the
current fit, clipping limits, and kept/rejected pixels.
**kwargs
Additional keyword arguments forwarded to the plotting helper if ``plot_steps=True``
(e.g., figure size, axis, title customization).
Returns
-------
None
The method updates the spectrum in place, setting:
- ``self._spec.cont`` : masked array of the final continuum model
- ``self._spec.cont_std`` : float, standard deviation of residuals on kept pixels
Notes
-----
- **Initialization:** The first iteration seeds the mask using the 16th–84th
percentile flux range of unmasked pixels, then fits the initial polynomial.
- **Clipping:** After each fit, residuals are computed and the standard
deviation is measured over currently kept pixels. New keep/reject limits are:
``low = model - abs_threshold[i] * std`` and
``high = model + emis_threshold[i] * std``.
- **Masking:** Existing pixel masks are honored; clipping only modifies the
continuum-selection mask on top of the original flux mask.
- **Smoothing:** When ``smooth_scale`` is provided, a boxcar (length
``smooth_scale``) is convolved with the flux prior to fitting; the continuum
itself is always evaluated on the original wavelength grid.
Examples
--------
Fit a three-iteration continuum with increasingly restrictive clipping:
>>> degrees = [1, 2, 2]
>>> thr_hi = [5.0, 3.0, 2.0] # emission-side thresholds (σ)
>>> thr_lo = [5.0, 3.0, 2.0] # absorption-side thresholds (σ)
>>> spec.fit.continuum(degrees, thr_hi, abs_threshold=thr_lo, smooth_scale=11)
Show diagnostic plots at each iteration:
>>> spec.fit.continuum([2, 2], [3.0, 2.0], plot_steps=True, title="Continuum fit")
Exclude known emission-line regions (defined in the rest frame) from the continuum fit:
>>> spec.continuum(degree_list=[3, 2], emis_threshold=[3.0, 2.0], exclude_intvls=[(4861, 4875), (6548, 6584)], rest_intvls=True)
"""
# Create a pre-Mask based on the original mask if available
input_wave = self._spec.wave.data
input_flux = self._spec.flux.data
input_mask = ~self._spec.flux.mask
# Correction if intervals are in the rest frame
z_corr = 1 + self._spec.redshift if rest_intvls else 1
# Adjust the mask to exclude wavelengt intervals
if exclude_intvls is not None:
exclude_intvls = np.asarray(exclude_intvls, dtype=float)
# Check the input has the right format
if exclude_intvls.ndim != 2 or exclude_intvls.shape[1] != 2:
raise ValueError('The argument "exclude_intvls" must be a list of (low, high) pairs')
# Loop through the wavelength intervals and add them to the mask
exclude_intvls = exclude_intvls * z_corr
i0 = np.searchsorted(input_wave, exclude_intvls[:, 0], side="right")
i1 = np.searchsorted(input_wave, exclude_intvls[:, 1], side="left")
for start, stop in zip(i0, i1):
input_mask[start:stop] = False
# Smooth the spectrum
if smooth_scale is not None:
smoothing_window = np.ones(smooth_scale) / smooth_scale
input_flux = np.convolve(input_flux, smoothing_window, mode='same')
# Loop through the fitting degree
abs_threshold = emis_threshold if abs_threshold is None else abs_threshold
for i, degree in enumerate(degree_list):
# First iteration use percentile limits for an initial fit
if i == 0:
low_lim, high_lim = np.nanpercentile(input_flux[input_mask], (16, 84))
mask_cont_0 = input_mask & (input_flux >= low_lim) & (input_flux <= high_lim)
cont_fit = continuum_model_fit(input_wave, input_flux, mask_cont_0, degree)
# Establishing the flux limits
std_flux = np.nanstd((input_flux - cont_fit)[input_mask])
low_lim, high_lim = cont_fit - abs_threshold[i] * std_flux, cont_fit + emis_threshold[i] * std_flux
# Add new entries to the mask
input_mask = input_mask & (input_flux >= low_lim) & (input_flux <= high_lim)
# Fit continuum
cont_fit = continuum_model_fit(input_wave, input_flux, input_mask, degree)
# Compute the continuum and assign replace the value outside the bands the new continuum
if plot_steps is True or (plot_steps == 'last' and (i == len(degree_list) - 1)):
default_ax_cfg = {'title':f'Continuum fitting ({i+1}/{len(degree_list)}), polynomial degree = {degree}'}
input_kwargs = kwargs.get('ax_cfg')
input_kwargs = default_ax_cfg if input_kwargs is None else default_ax_cfg.update(input_kwargs)
spec_continuum_calculation(self._spec, input_wave, input_flux, cont_fit, input_mask, low_lim,
high_lim, smooth_scale, exclude_intvls, ax_cfg=input_kwargs)
# Include the standard deviation of the spectrum for the unmasked pixels
self._spec.cont = np.ma.masked_array(cont_fit, self._spec.flux.mask)
self._spec.cont_std = np.std((input_flux - cont_fit)[input_mask])
return
class CubeTreatment(LineFitting):
def __init__(self, cube):
# Instantiate the dependencies
LineFitting.__init__(self)
# Lime spectrum object with the scientific data
self._cube = cube
self._spec = None
def spatial_mask(self, mask_file, fname, bands=None, fit_cfg=None, mask_list=None, line_list=None,
log_ext_suffix='_LINELOG', min_method='least_squares', profile=None, shape=None, cont_source='central',
err_from_bands=False, temp=10000.0, default_cfg_prefix='default', update_default=True,
line_detection=False, progress_output='bar', plot_fit=False, header=None, join_output_files=True):
"""
Measure lines across an IFS cube using one or more spatial masks.
This routine iterates over spaxels selected by binary masks (from ``mask_file``),
measures the requested lines for each spaxel, and writes the results to one or
more multi-extension FITS logs. Each spaxel’s measurements are saved in a
dedicated extension named ``"{j}-{i}{log_ext_suffix}"`` (e.g., ``"25-30_LINELOG"``).
The bands to use for measurements can be supplied globally via ``bands`` or,
per-mask, via entries in ``fit_cfg``. Line-fitting behavior is configured via
``fit_cfg`` with a multi-level override scheme (default → mask → spaxel).
Parameters
----------
mask_file : str or pathlib.Path or dict or numpy.ndarray
Source of spatial masks. Can be a FITS file produced by
:meth:`~lime.Cube.spatial_masking`, a dictionary mapping mask names to
boolean arrays, or a boolean array.
fname : str or pathlib.Path
Output path for the combined measurements log (or the base name, if generating one .fits file per mask).
bands : pandas.DataFrame, str, or pathlib.Path, optional
Bands table (or path) to use for all masks/spaxels. If ``None``, the
method will look for a per-mask bands path in ``fit_cfg``.
fit_cfg : dict or str or pathlib.Path, optional
Fitting configuration (dict or path to a TOML file). Supports a
three-level override hierarchy:
1) ``default_cfg_prefix`` section (global defaults),
2) mask-level section named after the mask (e.g., ``"MASK_A"``),
3) spaxel-level section named ``"{j}-{i}_line_fitting"``.
mask_list : list of str, optional
Subset of masks in ``mask_file`` to process. If ``None``, all masks are used.
line_list : list of str, optional
Subset of line labels to measure from the bands table. If ``None``,
all lines present in the bands are measured.
log_ext_suffix : str, optional
Suffix appended to each FITS extension name with results.
Default is ``"_LINELOG"``.
min_method : str, optional
Minimization algorithm used by :mod:`lmfit`. See
`lmfit.minimizer.Minimizer.minimize
<https://lmfit.github.io/lmfit-py/fitting.html#lmfit.minimizer.Minimizer.minimize>`_.
Default is ``"least_squares"``.
profile : str, optional
Profile identifier for fitting (e.g., ``"g"`` for Gaussian).
If ``None``, falls back to the line database defaults.
shape : str, optional
Line shape: ``"emi"`` for emission or ``"abs"`` for absorption.
If ``None``, falls back to database defaults.
cont_source : {'central', 'adjacent', 'fitted'}, optional
Method used to estimate the continuum level for line fitting.
- ``'central'`` — use the edges of the central line band (w₃–w₄).
- ``'adjacent'`` — use the adjacent continuum bands (w₁–w₂ and w₅–w₆).
- ``'fitted'`` — use a previously fitted continuum model from the spectrum.
The default is ``'central'``.
err_from_bands : bool or None, optional
If ``True``, estimate the pixel uncertainty from the continuum bands.
If ``None``, use the spectrum’s ``err_flux`` data or fall back to the
continuum regions if not available (False). The default value is None.
temp : float, optional
Electron temperature (K) for thermal broadening corrections.
Default is ``10000.0``.
default_cfg_prefix : str, optional
Section name in ``fit_cfg`` containing global defaults.
Default is ``"default"``.
update_default : bool, optional
If ``True`` (default), apply higher-level overrides by updating lower-level
dictionaries (shared keys only).
line_detection : bool, optional
If ``True``, run the continuum fitting and line threshodling to confirme the presence of lines before
measurements. The functions parameters must be specified in the ``fit_cfg`` (e.g., entries under ``"peaks_troughs"``, ``"continuum"``).
Default is ``False``.
progress_output : {"bar", "counter", None}, optional
Console progress reporting mode. Default is ``"bar"``.
plot_fit : bool, optional
If ``True``, render a plot of each spaxel’s fitted lines during processing.
Default is ``False``.
header : dict, optional
Extra FITS header keywords to add per spaxel page in the output logs.
If a key matching the extension name (e.g., ``"25-30_LINELOG"``) exists,
that dict is used; otherwise ``header`` is treated as a global header.
join_output_files : bool, optional
If multiple masks are processed, merge the per-mask logs into a single
FITS file named after ``fname``. When ``False``, keep one output file
per mask (named ``"{stem}_MASK-{mask}.fits"``). Default is ``True``.
Returns
-------
None
Results are written to disk as one or more FITS files. Each spaxel’s
measurements are stored in a binary table extension.
Notes
-----
- **Configuration hierarchy:** Values are resolved in the order
*default → mask → spaxel*. Higher levels **update** shared keys only if
``update_default=True``. Otherwise, the method applies a fallback protocol
where only the parameters explicitly defined in each section are used.
- **WCS headers:** Spatial WCS keywords are added to each extension header when
available, using the parent cube’s WCS metadata.
Examples
--------
Use a single bands table for all masks and join outputs:
>>> cube.obs.spatial_mask(
... mask_file="O3_masks.fits",
... fname="logs/o3_all_masks.fits",
... bands="my_bands.xlsx",
... fit_cfg="fit_config.toml",
... progress_output="bar",
... )
Provide bands per mask via the configuration and keep files separate:
>>> cube.obs.spatial_mask(
... mask_file="O3_masks.fits",
... fname="logs/o3_base.fits",
... fit_cfg="fit_config.toml",
... join_output_files=False,
... )
Limit measured lines and enable line detection:
>>> cube.obs.spatial_mask(
... mask_file="masks.fits",
... fname="logs/selected_lines.fits",
... line_list=["O3_5007A", "H1_4861A"],
... line_detection=True,
... )
"""
if bands is not None:
bands = check_file_dataframe(bands)
# Check if the mask variable is a file or an array
mask_maps = check_file_array_mask(mask_file, mask_list)
mask_list = np.array(list(mask_maps.keys()))
mask_data_list = list(mask_maps.values())
# Resolve the reference directory with respect to the config file if provided
if fit_cfg is not None and not isinstance(fit_cfg, dict):
cfg_dir = Path(fit_cfg).resolve().parent
else:
cfg_dir = Path.cwd()
# Check the mask configuration is included if there are no masks
input_masks = mask_list if bands is None else None
input_conf = check_fit_conf(fit_cfg, default_key=None, obj_key=None, update_default=update_default,
group_list=input_masks)
# Check if the output log folder exists
fname = Path(fname)
address_dir = fname.parent
if not address_dir.is_dir():
raise LiMe_Error(f'The folder of the output log file does not exist at {fname}')
address_stem = fname.stem
# Determine the spaxels to treat at each mask
spax_counter, total_spaxels, spaxels_dict = 0, 0, {}
for idx_mask, mask_data in enumerate(mask_data_list):
spa_mask, hdr_mask = mask_data
idcs_spaxels = np.argwhere(spa_mask)
total_spaxels += len(idcs_spaxels)
spaxels_dict[idx_mask] = idcs_spaxels
# Header data
hdr_coords = extract_wcs_header(self._cube.wcs, drop_axis='spectral') if self._cube.wcs is not None else None
# Loop through the masks
n_masks = len(mask_list)
mask_log_files_list = [address_dir/f'{address_stem}_MASK-{mask_name}.fits' for mask_name in mask_list]
for i in np.arange(n_masks):
# HDU_container
hdul_log = fits.HDUList([fits.PrimaryHDU()])
# Mask progress indexing
mask_name = mask_list[i]
mask_hdr = mask_data_list[i][1]
idcs_spaxels = spaxels_dict[i]
# Recover the fitting configuration
mask_conf = check_fit_conf(input_conf, default_cfg_prefix, mask_name)
# Load the mask log if provided
if bands is None:
# bands_file = Path(mask_conf['bands']).resolve()
# Get the section bands file when resolving bands or other relative paths from the config:
bands_file = (cfg_dir / mask_conf['bands']).resolve()
if bands_file.exists():
bands_in = load_frame(bands_file)
else:
err_msg = (f'Bands file not found at: {bands_file}.'
f'\n- Resolving from log section - entry: '
f'\n [{mask_name}_line_fitting]'
f'\n bands = {mask_conf["bands"]}')
raise LiMe_Error(err_msg)
else:
bands_in = bands
# Loop through the spaxels
n_spaxels = idcs_spaxels.shape[0]
n_lines, start_time = 0, time()
print(f'\nSpatial mask {i + 1}/{n_masks}) {mask_name} ({n_spaxels} spaxels)')
pbar = ProgressBar(progress_output, f'mask')
for j in np.arange(n_spaxels):
idx_j, idx_i = idcs_spaxels[j]
spaxel_label = f'{idx_j}-{idx_i}'
# Get the spaxel fitting configuration
spaxel_conf = input_conf.get(f'{spaxel_label}_line_fitting')
spaxel_conf = mask_conf if spaxel_conf is None else {**mask_conf, **spaxel_conf}
# Spaxel progress message
pbar.output_message(j, n_spaxels, pre_text="", post_text=f'(spaxel coordinate. {idx_j}-{idx_i})')
# Get spaxel data
spaxel = self._cube.get_spectrum(idx_j, idx_i, spaxel_label)
# Fit the lines
spaxel.fit.frame(bands_in, spaxel_conf, line_list=line_list, min_method=min_method,
line_detection=line_detection, profile=profile, shape=shape,
cont_source=cont_source, err_from_bands=err_from_bands,
temp=temp, progress_output=None, plot_fit=None,
obj_cfg_prefix=None, default_cfg_prefix=None, update_default=update_default)
# Count the number of measurements
n_lines += spaxel.frame.index.size
# Create page header with the default data
hdr_i = fits.Header()
# Add WCS information
if hdr_coords is not None:
hdr_i.update(hdr_coords)
# Add user information
if header is not None:
page_hdr = header.get(f'{spaxel_label}{log_ext_suffix}', None)
page_hdr = header if page_hdr is None else page_hdr
hdr_i.update(page_hdr)
# Save to a fits file
linesHDU = log_to_HDU(spaxel.frame, ext_name=f'{spaxel_label}{log_ext_suffix}', header_dict=hdr_i)
if linesHDU is not None:
hdul_log.append(linesHDU)
# Plot the fittings if requested:
if plot_fit:
spaxel.plot.spectrum(rest_frame=True)
# Save the log at each new mask
hdul_log.writeto(mask_log_files_list[i], overwrite=True, output_verify='ignore')
hdul_log.close()
# Computation time and message
end_time = time()
elapsed_time = end_time - start_time
print(f'\n{n_lines} lines measured in {elapsed_time/60:0.2f} minutes.')
if join_output_files:
output_comb_file = f'{address_dir/address_stem}.fits'
# In case of only one file just rename it
if len(mask_list) == 1:
mask_0_path = Path(mask_log_files_list[0])
mask_0_path.rename(Path(output_comb_file))
else:
print(f'\nJoining spatial log files ({",".join(mask_list)}) -> {output_comb_file}')
join_fits_files(mask_log_files_list, output_comb_file, delete_after_join=join_output_files)
return
