Gaussian Fitting

Gaussian Fitting#

This notebook shows how to use astropy and specutils to fit a Gaussian line profile to a Spectrum. The data comes from the calseq tutorial.

Dysh commands#

The following dysh commands are introduced (leaving out all the function arguments):

  filename = dysh_data()
  ta = Spectrum.read()
  p1 = ta.plot()
  p1.oshow(ta1)

Loading Modules#

We start by loading the modules we will use for the data reduction.

# These are the modules we will use for loading the data and fitting.
from dysh.spectra import Spectrum
from specutils.fitting import fit_lines
from astropy.modeling import models
from astropy import units as u
from dysh.log import init_logging

# These modules are used for file I/O
from dysh.util.files import dysh_data
from pathlib import Path

Setup#

We start the dysh logging, so we get more information about what is happening. This is only needed if working on a notebook. If using the CLI through the dysh command, then logging is setup for you.

init_logging(2)

# also create a local "output" directory where temporary notebook files can be stored.
output_dir = Path.cwd() / "output"
output_dir.mkdir(exist_ok=True)

Data Retrieval#

Download the example FITS spectrum, if necessary.

filename = dysh_data(example="nod-W/outputs/M82_ifnum_3_polavg.fits")

23:05:45.124 I url: http://www.gb.nrao.edu/dysh//example_data/nod-W/outputs/M82_ifnum_3_polavg.fits

Downloading M82_ifnum_3_polavg.fits from http://www.gb.nrao.edu/dysh//example_data/nod-W/outputs/M82_ifnum_3_polavg.fits

23:05:45.404 I Starting download...



23:05:45.656 I Saved M82_ifnum_3_polavg.fits to M82_ifnum_3_polavg.fits





Retrieved M82_ifnum_3_polavg.fits













Data Loading#


We use
Spectrum.read
to load the spectrum.






spec = Spectrum.read(filename, format="fits")









The loaded spectrum is now a
Spectrum,
with all its methods.
We can plot it, using
plot.






sp_plt = spec.plot(xaxis_unit="GHz")








../_images/11d64b5c65333afab826015f546361de0b4a7323b194b0bb8044e31311156058.png





Before fitting, lets remove the edges of the spectrum.
We only keep channels between 88.2 and 89.5 GHz.






spec = spec[88.2*u.GHz:89.5*u.GHz]









Plot again.






spec_plt = spec.plot(xaxis_unit="GHz")








../_images/8f94a7bb03281c253e6735a763b162ee44a924ee4273cd8d94c609c5e134fced.png





Now we are ready to start fitting.






Gaussian Fitting#


Here we show how to fit a Gaussian profile using astropy and specutils.
There are other options available, but we won’t cover them all.




Fitting a Single Gaussian#


We start fitting a single Gaussian line profile to the line on the left, at about 88.55 GHz.
We create a model, with starting values, and then use the specutils.fitting.fit_lines function to fit the model to the spectrum.
After the fit is done, we evaluate the best fit model over the spectral axis of the spectrum so we can plot it on top of the data.






g_init = models.Gaussian1D(amplitude=0.5*u.K, mean=88.55*u.GHz, stddev=0.1*u.GHz)
g_fit = fit_lines(spec, g_init)
y_fit = g_fit(spec.spectral_axis)









Now plot it.
We assign a “group id” (gid), so we can remove the line after.






spec_plt.axis.plot(spec.spectral_axis.to("GHz"), y_fit, gid="1gauss")
spec_plt.figure # This will show the figure here.








../_images/cb57a2aacd6e07983803c626000b7ddae67dad52cb2db788fffa4d7103d7f59c.png







Extracting Best Fit Parameters#


The best fit parameters are now properties of the output best fit model.
We can access them through the parameter names, for the Gaussian, amplitude, mean and stddev.
For example, the line amplitude






g_fit.amplitude









Parameter('amplitude', value=0.3670704696308813, unit=K)









and its formal error






g_fit.amplitude.std









np.float64(0.006236258674343534)













Fitting Two Gaussians#


We can also combine models.
In this case, we add a second Gaussian line profile to fit the second line at about 89.15 GHz.






g_init2 = models.Gaussian1D(amplitude=0.6*u.K, mean=89.15*u.GHz, stddev=0.1*u.GHz)
g_fit = fit_lines(spec, g_init + g_init2)
y_fit = g_fit(spec.spectral_axis)









Before plotting the new model with two Gaussians, we remove the previous line.






spec_plt.clear_lines("1gauss")









We repeat the plotting with the new model results.
This time, we add a label to the best fit model line (label="best fit").






spec_plt.axis.plot(spec.spectral_axis.to("GHz"), y_fit, label="best fit")
spec_plt.axis.legend()
spec_plt.figure
spec_plt.show()








../_images/ea741e6f478caa4cfd9ae0b832c67c0b86190e2bbce21d53e28484790b75a582.png





You can find more examples of how to use specutils.fitting.fit_lines in the specutils documentation.


In the case of a compound model, the model parameters get assigned an underscore and a number to differentiate them.
So the amplitude of the first Gaussian is stored in g_fit.amplitude_0 and that of the second Gaussian in g_fit.amplitude_1






g_fit.amplitude_0









Parameter('amplitude', value=0.3670154759025362, unit=K)













g_fit.amplitude_1









Parameter('amplitude', value=0.5278759700867988, unit=K)













Other Profiles and Further Examples#


Additional profiles and models are available through the astropy.modeling submodule.


The astropy documentation also provides examples of how to fit models with constraints, and examples of how to fit models, without using specutils.










Final Stats#


Finally, at the end we compute some statistics over a spectrum, merely as a checksum if the notebook is reproducible.






spec.check_stats(0.11606252 * u.K)









23:05:47.789 I Note: found 27 NaN (masked) values





23:05:47.789 I rms is OK 














                

              

              
              
              
              
                
              
            

            
            
              
                



  

  

     Contents