Core Functionality
This section contains the wrapper for reading datacubes and fitting them using the functions in LuciFit.py.
Core Functions
- class LuciBase.Luci(Luci_path, cube_path, output_dir, object_name, redshift, resolution, ML_bool=True, mdn=False)
This is the primary class for the general purpose line fitting code LUCI. This contains all io/administrative functionality. The fitting functionality can be found in the Fit class (Lucifit.py).
- close()
Functionality to delete Luci object (and thus the cube) from memory
- create_background_subspace(x_min=100, x_max=1900, y_min=100, y_max=1900, bkg_image='deep', n_components=50, n_components_keep=None, sigma_threshold=0.1, npixels=10, bkg_algo='detect_source', interpolation='nn')
This function will create a subspace of principal components describing the background emission. It will then interpolate the background eigenvectors over the entire field. Please see our paper () describing this methodology in detail. After running this code, LUCI will create a fits file containing the PCA coefficients for each pixel in the FOV which can then be used in other fitting functions by setting the optional argument bkg to PCA. LUCI will also save a numpy file with the principal components. We also create 3 plots that are saved in self.output_dir:
Map of background pixels in FOV
Visualization of first 10 principal components and mean
Scree plot
- Parameters:
x_min – Minimum x value (image coordinates) for segmentation region (Default 100)
x_max – Maximum x value (image coordinates) for segmentation region (Default 1900)
y_min – Minimum y value (image coordinates) for segmentation region (Default 100)
y_max – Maximum y value (image coordinates) for segmentation region (Default 1900)
bkg_image – 2D image used for background thresholding (Default deep image). Pass fits file.
n_components – Number of principal components to calculate (Default 50)
n_components_keep – Number of principal components to keep (Default n_components)
sigma_threshold – Threshold parameter for determining the background (default 0.1)
npixels – Minimum number of connected pixels in a detected group (default 10)
bkg_algo – Background algorithm to use (default ‘sourece_detect’; options: ‘source_detect’, ‘threshold’)
interpolation – Scheme for interpolation (default ‘nn’; options: ‘nn’, ‘linear’, ‘nearest’)
- Returns:
Fits file containing the PCA coefficients over the FOV PCA_eigenspectra: Numpy file containing the PCA eigenspectra
- Return type:
PCA_coeffs
- create_deep_image(output_name=None, binning=None)
Create deep image fits file of the cube. This takes the cube and sums the spectral axis. Then the deep image is saved as a fits file with the following naming convention: output_dir+’/’+object_name+’_deep.fits’. We also allow for the binning of the deep image – this is used primarily for astrometry purposes.
- Parameters:
output_name – Full path to output (optional)
binning – Binning number (optional integer; default=None)
- create_snr_map(x_min=0, x_max=2048, y_min=0, y_max=2064, method=1, n_threads=2, lines=[None], binning=1, bkgType=None, pca_coefficient_array=None, pca_vectors=None, pca_mean=None)
Create signal-to-noise ratio (SNR) map of a given region. If no bounds are given, a map of the entire cube is calculated.
- Parameters:
x_min – Minimal X value (default 0)
x_max – Maximal X value (default 2048)
y_min – Minimal Y value (default 0)
y_max – Maximal Y value (default 2064)
method – Method used to calculate SNR (default 1; options 1 or 2)
n_threads – Number of threads to use
lines – Lines to focus on (default None: For SN2 you can choose OIII)
binning – Bin to apply (default 1)
- Returns:
Signal-to-Noise ratio map
- Return type:
snr_map
- create_wvt(x_min_init, x_max_init, y_min_init, y_max_init, pixel_size, stn_target, roundness_crit, ToL, n_threads)
Written by Benjamin Vigneron.
Functionality to create a weighted Voronoi tesselation map from a region and according to arguments passed by the user. It creates a folder containing all the Voronoi bins that can then be used for the fitting procedure.
- Parameters:
x_min_init – Minimal X value
x_max_init – Maximal X value
y_min_init – Minimal Y value
y_max_init – Maximal Y value
pixel_size – Pixel size of the image. For SITELLE use pixel_size = 0.0000436.
stn_target – Signal-to-Noise target value for the Voronoi bins.
roundness_crit – Roundness criteria for the pixel accretion into bins
ToL – Convergence tolerance parameter for the SNR of the bins
n_threads – Number of threads to use
- detection_map(x_min=None, x_max=None, y_min=None, y_max=None, n_threads=1)
Method to call the detection algorithm. The detection algorithm works as follows: For each pixel,
Calculate the median spectrum for a 3x3 pixel region centered on the current pixel
Calculate the median spectrum for a 9x9 pixel region centered on the current pixel
Subtract the 9x9 spectrum from the 3x3 spectrum
Take the maximum value of this subtracted spectrum as the detection map value
In the end, we have a detection map of the maximum values.
If no bounds are added, we calculate over the entire map
Louis-Simon Guité
- Parameters:
x_min – Lower bound in x
x_max – Upper bound in x
y_min – Lower bound in y
y_max – Upper bound in y
n_threads – Number of threads (default 1)
- export_fits()
Export HDF file as fits file. The data will be saved in the same location as the cube (cube_dir) with the same name (object_name).
- extract_spectrum(x_min, x_max, y_min, y_max, bkg=None, binning=None, mean=False)
Extract spectrum in region. This is primarily used to extract background regions. The spectra in the region are summed and then averaged (if mean is selected). Using the ‘mean’ argument, we can either calculate the total summed spectrum (False) or the averaged spectrum for background spectra (True).
- Parameters:
x_min – Lower bound in x
x_max – Upper bound in x
y_min – Lower bound in y
y_max – Upper bound in y
bkg – Background Spectrum (1D numpy array; default None)
binning – Value by which to bin (default None)
mean – Boolean to determine whether or not the mean spectrum is taken. This is used for calculating background spectra.
- Returns:
X-axis (redshifted) and spectral axis of region.
- extract_spectrum_region(region, mean=False)
Extract spectrum in region. This is primarily used to extract background regions. The spectra in the region are summed and then averaged (if mean is selected). Using the ‘mean’ argument, we can either calculate the total summed spectrum (False) or the averaged spectrum for background spectra (True).
- Parameters:
region – Name of ds9 region file (e.x. ‘region.reg’). You can also pass a boolean mask array.
mean – Boolean to determine whether or not the mean spectrum is taken. This is used for calculating background spectra.
- Returns:
X-axis and spectral axis of region.
- static fit_calc(i, x_min, x_max, y_min, fit_function, lines, vel_rel, sigma_rel, cube_slice, spectrum_axis, wavenumbers_syn, transmission_interpolated, interferometer_theta, hdr_dict, step_nb, zpd_index, mdn, mask=None, ML_bool=True, bayes_bool=False, bayes_method='emcee', uncertainty_bool=False, nii_cons=False, bkg=None, bkgType=None, binning=None, spec_min=None, spec_max=None, initial_values=[False], obj_redshift=0.0, n_stoch=1, resolution=1000, Luci_path=None, pca_coefficient_array=None, pca_vectors=None, pca_mean=None)
Function for calling fit for a given y coordinate.
- Parameters:
i – Y coordinate step
lines – Lines to fit (e.x. [‘Halpha’, ‘NII6583’])
fit_function – Fitting function to use (e.x. ‘gaussian’)
vel_rel – Constraints on Velocity/Position (must be list; e.x. [1, 2, 1])
sigma_rel – Constraints on sigma (must be list; e.x. [1, 2, 1])
x_min – Lower bound in x
x_max – Upper bound in x
y_min – Lower bound in y
bkg – Background Spectrum (1D numpy array; default None)
bkgType – default None
binning – Value by which to bin (default None)
ML_bool – Boolean to determione whether or not we use ML priors
bayes_bool – Boolean to determine whether or not to run Bayesian analysis (default False)
bayes_method – Bayesian Inference method. Options are ‘[emcee’, ‘dynesty’] (default ‘emcee’)
uncertainty_bool – Boolean to determine whether or not to run the uncertainty analysis (default False)
nii_cons – Boolean to turn on or off NII doublet ratio constraint (default True)
initial_values – List of files containing initial conditions (default False)
spec_min – Minimum value of the spectrum to be considered in the fit (we find the closest value)
spec_max – Maximum value of the spectrum to be considered in the fit
obj_redshift – Redshift of object to fit relative to cube’s redshift. This is useful for fitting high redshift objects
n_stoch – The number of stochastic runs – set to 50 for fitting double components (default 1)
pca_coefficient_array – Array of PCA Coefficients (default None)
pca_vectors – Vectors corresponding to principal components (default None)
pca_mean – Mean vector from PCA analysis (default None)
- Returns:
all fit parameters for y-slice
- fit_cube(lines, fit_function, vel_rel, sigma_rel, x_min, x_max, y_min, y_max, bkg=None, bkgType=None, binning=None, bayes_bool=False, bayes_method='emcee', uncertainty_bool=False, n_threads=2, nii_cons=True, initial_values=[False], spec_min=None, spec_max=None, obj_redshift=0.0, n_stoch=1, pca_coefficient_array=None, pca_vectors=None, pca_mean=None)
Primary fit call to fit rectangular regions in the data cube. This wraps the LuciFits.FIT().fit() call which applies all the fitting steps. This also saves the velocity and broadening fits files. All the files will be saved in the folder Luci. The files are the fluxes, velocities, broadening, amplitudes, and continuum (and their associated errors) for each linespectrum_axis.
- Parameters:
lines – Lines to fit (e.x. [‘Halpha’, ‘NII6583’])
fit_function – Fitting function to use (e.x. ‘gaussian’)
vel_rel – Constraints on Velocity/Position (must be list; e.x. [1, 2, 1])
sigma_rel – Constraints on sigma (must be list; e.x. [1, 2, 1])
x_min – Lower bound in x
x_max – Upper bound in x
y_min – Lower bound in y
y_max – Upper bound in y
bkg – Background Spectrum (1D numpy array; default None)
bkgType – Type of background (default ‘standard’; options [‘standard’, ‘pca’])
binning – Value by which to bin (default None)
bayes_bool – Boolean to determine whether or not to run Bayesian analysis (default False)
'[emcee' (bayes_method = Bayesian Inference method. Options are)
'dynesty'] (default 'emcee')
uncertainty_bool – Boolean to determine whether or not to run the uncertainty analysis (default False)
n_threads – Number of threads to be passed to joblib for parallelization (default = 1)
nii_cons – Boolean to turn on or off NII doublet ratio constraint (default True)
initial_values – List of files containing initial conditions (default [False])
spec_min – Minimum value of the spectrum to be considered in the fit (we find the closest value)
spec_max – Maximum value of the spectrum to be considered in the fit
obj_redshift – Redshift of object to fit relative to cube’s redshift. This is useful for fitting high redshift objects
n_stoch – The number of stochastic runs – set to 50 for fitting double components (default 1)
pca_coefficient_array – Array of PCA Coefficients (default None)
pca_vectors – Vectors corresponding to principal components (default None)
pca_mean – Mean vector from PCA analysis (default None)
- Returns:
Velocity and Broadening arrays (2d). Also return amplitudes array (3D).
Examples
As always, we must first have the cube initialized (see basic example).
If we want to fit all five lines in SN3 with a sincgauss function and binning of 2 over a rectangular region defined in image coordinates as 800<x<1500; 250<y<1250, we would run the following:
>>> vel_map, broad_map, flux_map, chi2_fits = cube.fit_cube(['Halpha', 'NII6548', 'NII6583', 'SII6716', 'SII6731'], 'sincgauss', [1,1,1,1,1], [1,1,1,1,1], 800, 1500, 250, 750, binning=2)
- fit_entire_cube(lines, fit_function, vel_rel, sigma_rel, bkg=None, binning=None, bayes_bool=False, output_name=None, uncertainty_bool=False, n_threads=1)
Fit the entire cube (all spatial dimensions)
- Parameters:
lines – Lines to fit (e.x. [‘Halpha’, ‘NII6583’])
fit_function – Fitting function to use (e.x. ‘gaussian’)
vel_rel – Constraints on Velocity/Position (must be list; e.x. [1, 2, 1])
sigma_rel – Constraints on sigma (must be list; e.x. [1, 2, 1])
bkg – Background Spectrum (1D numpy array; default None)
binning – Value by which to bin (default None)
bayes_bool – Boolean to determine whether or not to run Bayesian analysis (default False)
output_name – User defined output path/name (default None)
uncertainty_bool – Boolean to determine whether or not to run the uncertainty analysis (default False)
n_threads – Number of threads to be passed to joblib for parallelization (default = 1)
- Returns:
Velocity and Broadening arrays (2d). Also return amplitudes array (3D).
- fit_pixel(lines, fit_function, vel_rel, sigma_rel, pixel_x, pixel_y, binning=None, bkg=None, bayes_bool=False, bayes_method='emcee', uncertainty_bool=False, nii_cons=True, spec_min=None, spec_max=None, obj_redshift=0.0, n_stoch=1, bkgType='standard', pca_coefficient_array=None, pca_vectors=None, pca_mean=None)
Primary fit call to fit a single pixel in the data cube. This wraps the LuciFits.FIT().fit() call which applies all the fitting steps.
- Parameters:
lines – Lines to fit (e.x. [‘Halpha’, ‘NII6583’])
fit_function – Fitting function to use (e.x. ‘gaussian’)
vel_rel – Constraints on Velocity/Position (must be list; e.x. [1, 2, 1])
sigma_rel – Constraints on sigma (must be list; e.x. [1, 2, 1])
pixel_x – X coordinate (physical)
pixel_y – Y coordinate (physical)
binning – Number of pixels to take around coordinate (i.e. bin=1 will take all pixels touching the X and Y coordinates.
bkg – Background Spectrum (1D numpy array; default None)
bkgType – Type of background (default ‘standard’; options [‘standard’, ‘pca’])
bayes_bool – Boolean to determine whether or not to run Bayesian analysis (default False)
bayes_method – Bayesian Inference method. Options are ‘[emcee’, ‘dynesty’] (default ‘emcee’)
uncertainty_bool – Boolean to determine whether or not to run the uncertainty analysis (default False)
nii_cons – Boolean to turn on or off NII doublet ratio constraint (default True)
spec_min – Minimum value of the spectrum to be considered in the fit (we find the closest value)
spec_max – Maximum value of the spectrum to be considered in the fit
obj_redshift – Redshift of object to fit relative to cube’s redshift. This is useful for fitting high redshift objects
n_stoch – The number of stochastic runs – set to 50 for fitting double components (default 1)
bkgType – Type of background (default ‘standard’; options [‘standard’, ‘pca’])
pca_coefficient_array – Array of PCA Coefficients (default None)
pca_vectors – Vectors corresponding to principal components (default None)
pca_mean – Mean vector from PCA analysis (default None)
- Returns:
Returns the x-axis (redshifted), sky, and fit dictionary
- fit_region(lines, fit_function, vel_rel, sigma_rel, region, bkg=None, binning=None, bayes_bool=False, bayes_method='emcee', output_name=None, uncertainty_bool=False, n_threads=1, nii_cons=True, spec_min=None, spec_max=None, obj_redshift=0.0, initial_values=[False], n_stoch=1, pixel_list=False)
Fit the spectrum in a region. This is an extremely similar command to fit_cube except it works for ds9 regions. We first create a mask from the ds9 region file. Then we step through the cube and only fit the unmasked pixels. Although this may not be the most efficient method, it does ensure the fidelity of the wcs system. All the files will be saved in the folder Luci. The files are the fluxes, velocities, broadening, amplitudes, and continuum (and their associated errors) for each line.
- Parameters:
lines – Lines to fit (e.x. [‘Halpha’, ‘NII6583’])
fit_function – Fitting function to use (e.x. ‘gaussian’)
vel_rel – Constraints on Velocity/Position (must be list; e.x. [1, 2])
sigma_rel – Constraints on sigma (must be list; e.x. [1, 2])
region – Name of ds9 region file (e.x. ‘region.reg’). You can also pass a boolean mask array.
bkg – Background Spectrum (1D numpy array; default None)
binning – Value by which to bin (default None)
bayes_bool – Boolean to determine whether or not to run Bayesian analysis (default False)
bayes_method – Bayesian Inference method. Options are ‘[emcee’, ‘dynesty’] (default ‘emcee’)
output_name – User defined output path/name
uncertainty_bool – Boolean to determine whether or not to run the uncertainty analysis (default False)
n_threads – Number of threads to be passed to joblib for parallelization (default = 1)
nii_cons – Boolean to turn on or off NII doublet ratio constraint (default True)
spec_min – Minimum value of the spectrum to be considered in the fit (we find the closest value)
spec_max – Maximum value of the spectrum to be considered in the fit
obj_redshift – Redshift of object to fit relative to cube’s redshift. This is useful for fitting high redshift objects
initial_values – List of files containing initial conditions (default [False])
n_stoch – The number of stochastic runs – set to 50 for fitting double components (default 1)
pixel_list – Boolean indicating if the user passes a 2D list containing pixel IDs for the region (default False)
- Returns:
Velocity and Broadening arrays (2d). Also return amplitudes array (3D).
Examples
As always, we must first have the cube initialized (see basic example).
If we want to fit all five lines in SN3 with a gaussian function and no binning over a ds9 region called main.reg, we would run the following:
>>> vel_map, broad_map, flux_map, chi2_fits = cube.fit_region(['Halpha', 'NII6548', 'NII6583', 'SII6716', 'SII6731'], 'gaussian', [1,1,1,1,1], [1,1,1,1,1],region='main.reg')
We could also enable uncertainty calculations and parallel fitting:
>>> vel_map, broad_map, flux_map, chi2_fits = cube.fit_region(['Halpha', 'NII6548', 'NII6583', 'SII6716', 'SII6731'], 'gaussian', [1,1,1,1,1], [1,1,1,1,1], region='main.reg', uncertatinty_bool=True, n_threads=4)
- fit_spectrum_region(lines, fit_function, vel_rel, sigma_rel, region, initial_values=[False], bkg=None, bayes_bool=False, bayes_method='emcee', uncertainty_bool=False, mean=False, nii_cons=True, spec_min=None, spec_max=None, obj_redshift=0.0, n_stoch=1)
Fit spectrum in region. The spectra in the region are summed and then averaged (if mean is selected). Using the ‘mean’ argument, we can either calculate the total summed spectrum (False) or the averaged spectrum for background spectra (True).
- Parameters:
lines – Lines to fit (e.x. [‘Halpha’, ‘NII6583’])
fit_function – Fitting function to use (e.x. ‘gaussian’)
vel_rel – Constraints on Velocity/Position (must be list; e.x. [1, 2, 1])
sigma_rel – Constraints on sigma (must be list; e.x. [1, 2, 1])
region – Name of ds9 region file (e.x. ‘region.reg’). You can also pass a boolean mask array.
initial_values
bkg – Background Spectrum (1D numpy array; default None)
bayes_bool – Boolean to determine whether or not to run Bayesian analysis
bayes_method – Bayesian Inference method. Options are ‘[emcee’, ‘dynesty’] (default ‘emcee’)
uncertainty_bool – Boolean to determine whether or not to run the uncertainty analysis (default False)
mean – Boolean to determine whether or not the mean spectrum is taken. This is used for calculating background spectra.
nii_cons – Boolean to turn on or off NII doublet ratio constraint (default True)
spec_min – Minimum value of the spectrum to be considered in the fit (we find the closest value)
spec_max – Maximum value of the spectrum to be considered in the fit
obj_redshift – Redshift of object to fit relative to cube’s redshift. This is useful for fitting high redshift objects
n_stoch – The number of stochastic runs – set to 50 for fitting double components (default 1)
- Returns:
X-axis and spectral axis of region.
- fit_wvt(lines, fit_function, vel_rel, sigma_rel, bkg=None, bayes_bool=False, uncertainty_bool=False, n_threads=1, initial_values=[False], n_stoch=1, stn_target=10)
Function that takes the wvt mapping created using self.create_wvt() and fits the bins. Written by Benjamin Vigneron
- Parameters:
lines – Lines to fit (e.x. [‘Halpha’, ‘NII6583’])
fit_function – Fitting function to use (e.x. ‘gaussian’)
vel_rel – Constraints on Velocity/Position (must be list; e.x. [1, 2, 1])
sigma_rel – Constraints on sigma (must be list; e.x. [1, 2, 1])
bkg – Background Spectrum (1D numpy array; default None)
bayes_bool – Boolean to determine whether or not to run Bayesian analysis
uncertainty_bool – Boolean to determine whether or not to run the uncertainty analysis (default False)
n_threads – Number of threads to use
initial_values – Initial values of velocity and broadening for fitting specific lines (must be list)
n_stoch – The number of stochastic runs – set to 50 for fitting double components (default 1)
stn_target – Target signal to noise ratio (default 10)
- Returns:
Velocity, Broadening and Flux arrays (2d). Also return amplitudes array (3D) and header for saving figure.
- heliocentric_correction()
Calculate heliocentric correction for observation given the location of SITELLE/CFHT and the time of the observation
- read_in_cube()
Function to read the hdf5 data into a 3d numpy array (data cube). We also translate the header to standard wcs format by calling the update_header function. Note that the data are saved in several quadrants which is why we have to loop through them and place all the spectra in a single cube.
- skyline_calibration(Luci_path, n_grid, bin_size=30)
Compute skyline calibration by fitting the 6498.729 Angstrom line. Flexures of the telescope lead to minor offset that can be measured by high resolution spectra (R~5000). This function divides the FOV into a grid of NxN spaxel regions of 10x10 pixels to increase the signal. The function will output a map of the velocity offset. The initial velocity guess is set to 80 km/s. Additionally, we fit with a simple sinc function.
- Parameters:
n_grid – NxN grid (int)
Luci_path – Full path to LUCI (str)
bin_size – Size of grouping used for each region (optional int; default=30)
- Returns:
Velocity offset map
- slicing(lines, bkg=None)
Slices the cube along the spectral axis around the specified emission line wavelengths. The input velocity dispersion serves as a rough estimate of the width of the spectral lines to make sure the slices cover a wide enough range on each side of the emission lines. The slices are saved in a new folder for each input emission line.
We currently have the following lines available:
‘Halpha’: 656.280, ‘NII6583’: 658.341, ‘NII6548’: 654.803, ‘SII6716’: 671.647, ‘SII6731’: 673.085, ‘OII3726’: 372.603, ‘OII3729’: 372.882, ‘OIII4959’: 495.891, ‘OIII5007’: 500.684, ‘Hbeta’: 486.133, ‘OH’: 649.873, ‘HalphaC4’: 807.88068, ‘NII6583C4’: 810.417771, ‘NII6548C4’: 806.062493, ‘OIII5007C2’: 616.342, ‘OIII4959C2’: 610.441821, ‘HbetaC2’: 598.429723, ‘OII3729C1’: 459.017742, ‘OII3726C1’: 458.674293, ‘FeXIV5303’: 530.286, ‘NI5200’: 520.026, ‘FeVII5158’: 515.89, ‘HeII5411’: 541.152
- Parameters:
lines – Lines to extract (e.x. [‘Halpha’, ‘NII6583’])
bkg – Background Spectrum (default None)
- visualize()
Wrapper function for LUCI.LuciVisualize()
- wvt_fit_region(x_min_init, x_max_init, y_min_init, y_max_init, lines, fit_function, vel_rel, sigma_rel, stn_target, pixel_size=0.436, roundness_crit=0.3, ToL=0.01, bkg=None, bayes_bool=False, uncertainty_bool=False, n_threads=1, n_stoch=1, initial_values=[False])
Functionality to wrap-up the creation and fitting of weighted Voronoi bins.
- Parameters:
x_min_init – Minimal X value
x_max_init – Maximal X value
y_min_init – Minimal Y value
y_max_init – Maximal Y value
lines – Lines to fit (e.x. [‘Halpha’, ‘NII6583’])
fit_function – Fitting function to use (e.x. ‘gaussian’)
vel_rel – Constraints on Velocity/Position (must be list; e.x. [1, 2, 1])
sigma_rel – Constraints on sigma (must be list; e.x. [1, 2, 1])
stn_target – Signal-to-Noise target value for the Voronoi bins.
pixel_size – Pixel size of the image. For SITELLE use pixel_size = 0.0000436.
roundness_crit – Roundness criteria for the pixel accretion into bins
ToL – Convergence tolerance parameter for the SNR of the bins
bkg – Background Spectrum (1D numpy array; default None)
bayes_bool – Boolean to determine whether or not to run Bayesian analysis
uncertainty_bool – Boolean to determine whether or not to run the uncertainty analysis (default False)
n_threads – Number of threads to use
initial_values – Initial values of velocity and broadening for fitting specific lines (must be list;
[velocity (e.x.)
[False]) (broadening]; default)
n_stoch – The number of stochastic runs – set to 50 for fitting double components (default 1)
- Returns:
Velocity, Broadening and Flux arrays (2d). Also return amplitudes array (3D).