SAFFRON: Spectral Analysis Fitting Framework, Reduction Of Noise.
Project description
SAFFRON
Spectral Analysis, Fitting Framework, Reduction Of Noise
(Unpackaged version)
Data fitting pipeline adapted to SPICE instrument onboard SolaOrbiter. I will apreceate any suggestions to improve the quality of the content If you have any question you can contact the author: Mzerguat Slimane
Requirements
- Python<=3.9.15 (I suggest to use pyenv to change versions easily )
- Chianti database (Download)
- Add the variable
XUVTOP
to the path of the database (After extraction)
Linux :
export XUVTOP=/home/../pathTo/Chianti_Database
Windows :
search Environment Variables in start panel
In the Environment Variables window, under the "System variables" section, scroll down to find the "Path" variable.
Click on "New..." to add a new system variable.
Set Variable Name and Value:
In the "New System Variable" dialog, set the variable name as XUVTOP.
In the "Variable value" field, enter the path or database value you want to assign to XUVTOP.
Save the Variable:
Click "OK" to save the new system variable.
python :
import os
# Set the XUVTOP environment variable, will be gone by the end of the script
os.environ['XUVTOP'] = 'path/to/your/database'
Install
- Install the right version of python
pyenv install 3.9.15
- Create and activate your environment
pyenv virtualenv 3.9.15 SPICE_SlimPy
pyenv activate SPICE_SlimPy
MAC: MAC is upset with me using shared memory. The package is not working on MAC machines :(. Oh,the quirks of using MAC machines for multi-processing tasks within this package! It's like expecting a cat to enjoy water. But seriously, who embarks on scientific adventures with a MAC anyway? 😄
- install requirement
The library hasn't bee packaged yet so it's better to put in the parent folder to be able to use it. This is going to be changed eventually.
Inside SlimPy folder do:
pip install -r requirements.txt
Tutorial
Basic Fitting
from SlimPy.manager.manager import Manager
from SlimPy.utils import get_input_template
get_input_template(where='./input.json')
#Go inside JSON file and put the list of L2 files to use
session = Manager("./input.json")
# print(session) if you want
session.build_files_list()
#you can do session.fitinits = 4 before calling the next command if you want to save the graphs of initial parameters preparation in side ./tmp/ as png
session.build_rasters()
# change sessions.rasters[i].init_params if not good or
# change the interval if there are artifacts in the borders of spectral windows (you can see if there are by using SAFFRAN.utils.quickview(L2_path))
session.run_preparations() #clean,convolve,denoise,despike,estimate error
session.fit_manager() #fit and save
# the fitting will be saved in the directory that is specified in input.json under the keyword "data_save_dir" and the naming structure is under the keyword "data_filename"
Code details
The odds are that you will need no major adjustments for the fitting as the input parameters are most likely tuned out.
JSON Input In the beginning you will need a input json file that contains all the parameters needed for the analysis part.
from SlimPy.utils import get_input_template
get_input_template(where='./input.json')
Inside the json the most important keys are:
- "SELECTION_MODE": list, folder, date intervale
- "files" : list of files in case SELECTION_MODE is a list
- "data_filename" : "use ::SAMENAME to replace the placeholder, ::SAMENAMEL2.5 is to change L2 to L2.5 by the name of the fits input file or ::CONV to replace by the convolution level(s) with - as separator or ::TIME to replace it by fitting start time"
- "data_save_dir" : The folder to save data in ('.'=in the current directory)
- "window_size" : the size of the window to fitted in y,x (not x,y) in pixel coordinates. [0,Null] = all the pixel in a given direction.
- "Jobs" : Number jobs while fitting (It's 30 by default so make sure to change it to the right number). Else your CPU will pop out of the machine and sit next to you.
Initiating and preparing Manager
from SlimPy.manager import Manager
session = Manager("./input.json")
session.build_files_list()
- You can visualize and change this list by accessing the variable session.selected_fits
Building rasters (<RasterFit>)
Building rasters mean you will have a list of rasters that each have their parameters passed from the <Manager>. It's also auto calculating the initial parameters based on a catalogue the library. If you want to see the results of initial parameters set the verbose for geninits_verbose >3 or <-2. the plots are going to be in a ./tmp
folder with dates of creation in their names.
session.build_rasters()
Now you have a new variable in your <Manager> (session.rasters) that have a list of <RasterFit> objects. if you are not satisfied of the initial parameters algorithm you can adjust it manually.
you only have to access these 2 parameters:
(raster.init_params): list of array for each window contains all the parameter for a gaussian [int,wvl,wid,int,wvl,wid....,Bg]
(raster.quantities) the same list of the same size as init_params however, it's a list of characters for each parameter. 'I' for intensity, "x" for wavelength, "s" for width
Inside each raster there is a list of windows as <WindowFit> object, but, it's going to be explained later
data cleaning denoising, error estimation (for fitting weights), despising (remove cosmics), convolving (improving S2N ratio) is not yet applied so we are going to do so by calling.
session.run_preparations()
this is also possible by calling the same method inside one of the rasters or even one of the windows
session.run_preparations() #Loop over all rasters of the session that will loop all over the windows
session.rasters[i].run_preparations() #loop over all the windows of the raster
session.rasters[i].windows[j].run_preparations() #Apply on the windows only
Data preparation is a little time consuming so it will not run if it's called twice in the second time. unless if you set redo=True as an argument for .run_preparations()
Fitting data
session.fit_manager()
you can also run fitting for only one raster or even one window
session.rasters[i].fit_raster() #i index of the raster you want to fit
session.rasters[i].windows[j].fit_window() #j index of the window inside raster of index i you want to fit
I tried a progress bar buuut Hmmm it's not yet well coded -_- . So it looks stupid. But, it runs independently and doesn't affect the fitting processes.
Managing Line Catalog
SAFFRON employs a dynamic approach to identify line parameters essential for initiating the fitting process. This functionality is anchored in an integral component of SAFFRON: the line catalog. Our module not only comes with a pre-configured internal catalog but also offers the flexibility to customize or even replace it with your catalog, catering to all lines you deem necessary, whether they currently exist or not. The existing catalog predominantly features main lines optimized for composition lines ...
Getting Started with the LineCatalog Class
The LineCatalog
class is your gateway to interacting with the line catalog. Here's how you can get started:
from SAFFRON.catalog import LineCatalog
# Initialize the catalog; you can specify a custom file location for your catalog
catalog = LineCatalog(file_location="path/to/your/catalog.json") #if not specified it will upload the default catalog
Loading and Viewing the Catalog
Upon initialization, the catalog is automatically loaded. You can directly interact with it to view the lines and spectral windows it contains.
# Access the catalog data
lines = catalog.get_catalog_lines()
windows = catalog.get_catalog_windows()
Modifying the Catalog
SAFFRON's LineCatalog
facilitates various operations to tailor the catalog to your needs, including adding or removing lines and spectral windows.
You can change the catalogue manually by going to SAFFRON.catalog.SPICE_SpecLines.json
or with code
-
Adding a New Line
catalog.add_line(name="NewLine", wvl=123.45)
-
Removing a Line
catalog.remove_line(name="OldLine")
-
Adding a Spectral Window Important Note: If there is a line that doesn't exist yet you can not use it to create a new window. the line must be first added to the list of lines with its wavelength value
catalog.add_window(lines=["Line1", "Line2"], max_line="Line2")
-
Removing a Spectral Window
catalog.remove_window(lines=["Line1", "Line2"])
Saving Changes
Any modifications made can be persisted back to the JSON file, ensuring your catalog stays up-to-date.
catalog.dump(new_path="path/to/save/your/updated_catalog.json")
Using your catalog
In Manager.build_raster
add the argument catalog_location
that has the location of your catalog. If you modified the default catalog, no need to specify the loaction
Important to know
- Ensure all new lines added have a unique name within the catalog.
- When adding a spectral window, the
max_line
should be the most intense line within that window. - Use descriptive names for lines to maintain clarity and ease of identification.
With SAFFRON's line catalog management, you're equipped to customize your spectral line analysis to fit your unique requirements, enhancing the flexibility and accuracy of your scientific exploration.
Locking Technique of Spectral Lines
Overview
The locking technique can be an essential part of the spectral line fitting in some cases in SAFFRON module, designed to enhance the accuracy of distinguishing between blended lines. This method is particularly useful in complex spectra where lines such as S_iv 750 and mg_ix 749 are blended in each other. Thus, we try to lock the position of mg_ix 706 with mg_ix 749 (and sometimes: S_iv 750 with S_iv 748) to increase the accuracy on the position. This would also cause a better estimation on the other parameters
Application
if the lines to lock are in two different windows we should first fuse those windows Example: (Approximate wavelengths)
-
session.raster[i],window[0] have
0 . O iii ~702.7 Å.
1 . O iii ~703.2 Å.
2 . Mg ix ~706 Å.
-
session.raster[i],window[1] have
0 . S iv ~748 Å.
1 . Mg ix ~749.5 Å.
2 . S iv ~750 Å.
If we want to lock the two sulfur lines or oxygen lines we don't need to fuse the windows but locking mg_ix lines need to fuse window[0] with window[1]
Windows Fusion
session.fuse_windows([0,1]) #0 is windiow 0 and 1 is window 1
if shape of the argument is ($N_{selected _files}$x2) there will be fusions different each file in the session
you can also fuse by selecting a raster directly
session.rasters[i].fuse_windows(0,1) #0 is windiow 0 and 1 is window 1
Now that you have fused the windows you can access this fused window by calling session.rasters[i].fused_windows[0]
Important note: Once you fuse two windows they are automatically excluded from fitting and the resulting fused window is automatically included. However the component windows are not deleted you can specially fit them explicitly by running session.rasters[i].windows[0<or >].fit_window()
Window locking
After Fusion the index of each line in the fused window will be:
-
session.raster[i],window[0] have 0 . O iii ~702.7 Å.
1 . O iii ~703.2 Å.
2 . Mg ix ~706 Å.
3 . S iv ~748 Å.
4 . Mg ix ~749.5 Å.
5 . S iv ~750 Å.
and now we use
Manager.set_lock_protocol(window_type,window_index,lock_line1_index,lock_line2_index,loc_distance)
to specify the locking protocol we want to use.
window_type
: "fuse" if the protocol in a fused window, "solo" if the protocol is in a single window.
window_index
: once selected window type. you select the order of the selected window in their respective list.
lock_line1_index
: the index of the leading line for locking.
lock_line2_index
: the index of the following line for locking into the leading one.
loc_distance
: The distance in between the two lines that it will remain the same during fitting.
session.set_lock_protocol("fuse", 0, 2, 4, (749.54-706.02))
session.set_lock_protocol("fuse", 0, 3, 5, (750.22-748.40))
Adaptive Convolution
TODO
How Denoise works?
TODO Credit: Frédéric Auchère
How Despike works?
TODO Credit: Frédéric Auchère
How Sigma is estimated?
Credit: Eric Buchlin
- Estimate errors on data: Ask Eric Buchlin.
- Estimation after convolution: $\frac{\sqrt{\Sigma_{Conv_pixels}{\Delta I_{val}^2 }}}{N_Conv_pixels}$
- with after denoise: Nothing done. That means we overshoot error values a bit because of denoise
Postprocessing (Composition maps)
This part is adapted from fiplcr module to SAFFRON output data. The FIP maps are generated using Linear Combination Method (LCR)
Using <SPICEL3Raster>
L3data_dir = './<output_data>/'
LFLines = ('s_4',750),('s_5',786), #Low Fip lines (needed to compute the FIP maps using LCR the choice is conditioned read the paper for more infos)
HFLines = ('n_4',765.15),('n_3',991.59) #High Fip lines (needed to compute the FIP maps using LCR the choice is conditioned read the paper for more infos)
from SAFFRON.postprocessing import SPICEL3Raster
L3_raster = SPICEL3Raster(folder_path = _con_L3_data_dir) #Generating L3 raster: loading outputs and computing the radiance maps and the error.
L3_raster.gen_compo_LCR(LFLines =LFLines,HFLines =HFLines)# optimizing the linear combination and computing FIP maps and the error values. Add "suppressOutput=True" if you want no graphs
Finally you have FIP maps
L3_raster.FIP
2D numpy array FIP map.
L3_raster.FIP_err
2D numpy array FIP error map.
L3_raster.show_all_wvls()
return all available lines' wavelengths in order.
L3_raster.find_line(wvl)
return the a line as SAFFRON.postprocessing.SPECLine
object with the closest wavelength to wvl.
L3_raster.lines
list of lines of type SAFFRON.postprocessing.SPECLine
in the raster.
L3_raster.lines[i].wavelength
Return the wavelength of the line i.
L3_raster.lines[i].observatory
Return the observatory of the line i.
L3_raster.lines[i].instrument
Return the instrument of the line i.
L3_raster.lines[i].ion
Return the ion of the line i.
L3_raster.lines[i].line_id
Return the line_id of the line i.
L3_raster.lines[i][par]
2D Return array parameter of line i with par in ['int' and/or 'wav' and/or 'wid' and/or 'rad' and/or 'int_err' and/or 'wav_err' and/or 'wid_err' and/or 'rad_err' ]
.
L3_raster.lines[i].header[par]
Return astropy header object of line i with par in ['int' and/or 'wav' and/or 'wid' and/or 'int_err' and/or 'wav_err' and/or 'wid_err' ]
(no "rad" neither "rad_err").
L3_raster.lines[i].plot(params='all',axes =None,add_keywords = False)
Plot a parameter or a set of parameters ['int' and/or 'wav' and/or 'wid' and/or 'rad' and/or 'int_err' and/or 'wav_err' and/or 'wid_err' and/or 'rad_err' ]
. If axes is not None than it should be a 1D Iterable with size equal to the number of parameters to plot.
Citation
If you use SAFFRON in your research, please cite the following:
Slimane MZERGUAT. SAFFRON: Spectral Analysis Fitting Framework, Reduction Of Noise, Version 1.0.0, 2024. Available at: https://github.com/slimguat/SAFFRON.
Project details
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
Built Distribution
Hashes for saffron_spice-0.1.1-py3-none-any.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | afa1ccd8f88ca3b8e81da8463d86d79aea037fbd41d2b994c5e272aeff872ed7 |
|
MD5 | 3db16021c443fdd3958822b57a4de7b2 |
|
BLAKE2b-256 | 7f2427ce09d8d5ae14beac376d0f7572bb427bb4da4940d2f34a7a173667afb1 |