Fitting Models to Spectra

Whew! This was a little tougher than expected, and I struggled with it quite a bit. The next phase of my research will be focusing on determining the relative strengths of emission and absorption lines in order to characterize the classification parameters of observed galaxies. For this, I will eventually be using a “Baldwin-Phillips-Terlevich” Diagram (BPT).  By measuring the ratios of the observed selected emission/absorption lines, we can empirically calculate indicators of the primary excitation mechanism in the observed galaxies. A great, and easy-to-read rundown of the application of BPT diagrams can be found here.  But, first, you have to determine the location of the targeted emission/absorption lines and measure their intensity.  This is where fitting models to spectra becomes handy.

Fig 1: Line Fitting For Hydrogen Alpha Line. The observed flux, on the right, is nearly covered by the model generated, visually indicating an accurate fit. The continuum spectrum uses a median-smoothing algorithm to produce a generic continuum. An unfortunate drawback to this method results in a less accurate fit when approaching the bounds of the spectrum. By confining the measured Hydrogen Alpha region to within 2-sigma, a more accurate calculation of the equivalent width can be performed.
Fig 2: Line Fitting for [OIII] Line. Again, this turned out to be a really good fit.
Fig 3: Line Fitting for Hydrogen Beta. Here, you can see that there is some deviation in the fit of the model, but is pretty consistent within 2-Sigma.

I thought about writing a tutorial, but a great one already exists. I found this “Quick Fit” modeling tutorial to be extremely useful when working with model parameters, especially when performing error analysis. Once you feel comfortable with the basics, I recommend moving on to the “User-defined Model” tutorial. In this tutorial, you’ll learn the tools necessary to fit Gaussian lines to emission/absorption lines using more complex models. Using this method, you can fit a Gaussian line to the emission/absorption lines, AND fit a one-degree polynomial to estimate the local continuum. Super useful! One thing to mention that wasn’t included in the tutorial: the model fitters give you the ability to have certain information returned to you, including a covariance matrix for the error in your fit parameters. Let’s look at In[22] from the User-defined Model tutorial.

Screen Grab of In[22] from the tutorial

Using this as a starting point, we can add a few lines of code to extract more information. Add the following lines to your code.

print(fitter.fit_info)

The output from ‘fitter.fit_info’.

There are an awful lot of numbers there. Let’s just focus on the parameter covariance matrix. We can call the keyword ‘param_cov’ from the fitter to retrieve the matrix. Next, we can “diagonalize” the matrix, and take the square-root of the values to determine the uncertainty in the associated parameters and then print the parameters associated with those errors..

print(np.sqrt(np.diag(fitter.fit_info['param_cov']))
print(compound_fit_bounded.param_names)
Errors in the fit are show in Out[9] and their associated parameters are shown in Out[10].

Armed with this information, you can ensure that your model is statistically sound.

Posted by Will

Painting A Picture With Plots

These figures represent a selection of the host galaxies within the SAGA Survey. Plotted are the positions of each host galaxy and their respective satellite galaxies. The common name of the galaxy, the virial radius in arcminutes, and the number of satellite galaxies are included in the title of the plot. The virial radius is defined by the red circle, while the smaller circles represent the location of the satellite galaxies and their recessional velocities in comparison to their hosts. Of note, there are several “clumps” of satellite galaxies that appear to have similar recessional velocities, and this is worth further investigating.

But, let’s get down to the basics. How do we go about generating these plots? Luckily for us, we can do this using completely open-source software and publicly available data. To begin, we’ll need a few things:

That’s it. We these five things, you can generate a usable, simple figure. Let’s begin by retrieving a FITS file. Head over to the DSS website, where you’ll be met with a homepage that is built for those astronomers who wax nostalgic for the 1990s.

Image above: Landing page from the archaeological record.

For the purpose of this example, we’ll be using galaxy ‘NSA 32’. As shown above, enter the galaxy designation in the ‘Object name’ field and click ‘GET COORDINATES’. At this point, the galaxy will have its designator cross-identified by either the SIMBAD or NED astronomical database, and then return the selected galaxy’s name as it appears in the database along with the Right Ascension and Declination automatically filled out in the form. Next, change the ‘Height’ and ‘Width’ fields to 60 arcminutes. Ensure the ‘File Format’ is set to ‘FITS’. At this point, your form should be the same as the above image. Now, click ‘RETRIEVE IMAGE’ and your download should begin. I HIGHLY RECOMMEND saving it with its galaxy designation as the file name. Make note of where you save it, because we’re going to need it.

Now, we can begin with the coding part. First we’ll import the necessary libraries into our Python Environment. Explanations for each library are included in the comments. I find that one can never have enough comments, and it has served me well so far. Especially for code that I haven’t looked at in ages.

import os #allows for the code to use system-dependent functionalities
import matplotlib.pyplot as plt #used for generating graphs/plots/figures
from astropy.io import fits #used to read in/out FITS formatted data
%autosave 45 #when used in a Jupyter Notebook, this will autosave your work every 45 seconds.

Next, we’ll open our downloaded FITS file using the imported FITS module. I usually use the following method.

file = ("E:\Use_Absolute_Path") #location of FITS file
hdu = fits.open(file) #opens the FITS, assigns it to 'hdu'
print(hdu.info()) #prints the contents of the file
Contents of the file are printed to the output.

It is absolutely necessary to use the absolute path when dealing with a singular FITS file. *HOT TIP* For Windows users, hold “Shift” and “Right-Click” the file, then select “Copy as path”. This will copy the absolute path to your clipboard and you can paste that into the code. Now that we’ve loaded the FITS file, we can take a look at what it contains. Depending on the type of FITS file you encounter, there can be several segments called Header/Data Units (HDUs) that are contained within the structure. For this particular file, we will be dealing with the Primary HDU, and as when dealing with arrays, the Primary HDU is at position [0]. One of the first things I like to do is to view the header of the FITS file’s Primary HDU.

hdu[0].header
Snippet of Header Output

The header contains a lot of useful information. We’ll return to this shortly. But, for now, let’s take a look at the image.

image = hdu[0].data #extracts the image data array from the HDU
fig = plt.figure() #initializes a new figure
plt.imshow(image, cmap = 'Greys') #displays data as an image. cmap is 'colormap'
The image extracted from the FITS file

We can see NGC 2967 near the center of the field. As it stands, this figure doesn’t provide much information. Looking at the X and Y axes, the numbers shown there are actually the pixel dimensions of the image. You can verify this by looking at ‘NAXIS1’ and ‘NAXIS2’ in the header of the file. Let’s spruce this figure up a bit by adding a title, axis labels, inward-facing tick marks, and the standard convention for measuring coordinates: Right Ascension and Declination. To get the RA/DEC from the FITS file, we’ll import an additional module from Astropy: WCS. We’ll modify the previous block of code to achieve this.

from astropy import wcs #to use WCS header information

image = hdu[0].data #extracts the image data array from the HDU
w = wcs.WCS(hdu[0].header) #extracts the WCS keywords from the header
fig = plt.figure() #initializes a new figure

#The following portion initiates a 1 by 1 subplot, with the image in the
#first position (1,1,1) and plotted with the RA/DEC using 'projection' and the WCS keywords from the header.

ax1 = fig.add_subplot(1,1,1, projection = w) #initiates subplot with WCS data
ax1.tick_params(direction = 'in') #changes the direction of the tick marks

plt.imshow(image, cmap = 'Greys') #displays data as an image. cmap is 'colormap'
plt.title('NGC 2967') #title for the plot
plt.ylabel('DEC') #label for the y-axis
plt.xlabel('RA') #label for the x-axis.
 

We’re making progress here, but the DEC and RA are using different units of measure. It’s not necessarily required for them to be uniform, as you’ll usually see these units being used for their respective coordinate. If you would like to have both axes use the same units, just add the following lines to your code.

from astropy import units as u #imports common astronomical units

#the following line forces degrees(u.deg) on the x-axis(.coords[0]) or y-axis (.coords[1])
ax1.coords[0].set_format_unit(u.deg) #forces degrees on x-axis

Lastly, we can add a circle to the plot to readily identify the object of interest. We’ll first import a new module that allows for super-imposing shapes on top of a plot, and then we’ll use information from the file header to correctly place the shape. Using the Circle package from Matplotlib, we have to provide three pieces of information: the x-coordinate of the center, the y-coordinate of the center, and the desired radius of the circle. In order to get the (x,y) position of the center, we’ll use two keywords from our header. There are two methods of doing this, but I’ll only be using the ‘pixel method’ for simplicity. Returning to the header of the file, as you scroll down you will find information designating the reference pixels for the object of interest.

from matplotlib.patches import Circle #allows for use of circle shape

ra,dec = hdu[0].header['CRPIX1'], hdu[0].header['CRPIX2'] #calls RA/DEC using keywords 'CRPIX1' and 'CRPIX2'

#the following line initiates the circle centered on the galaxy, with a edge color of red, and no face color to allow for the galaxy to be seen.
circle =Circle((ra,dec),1000, edgecolor='r',facecolor = 'none') 

fig = plt.figure()
w = wcs.WCS(hdu[0].header)
ax1 = fig.add_subplot(1,1,1,projection = w)
plt.imshow(hdu[0].data, cmap = 'Greys')
plt.title('NGC 2967')
plt.ylabel('DEC')
plt.xlabel('RA')
ax1.tick_params(direction = 'in')
ax1.coords[0].set_format_unit(u.deg)

ax1.add_patch(circle) #adds the circle patch to the image;


The plotting of the satellites proved to be a little tricky when using the colormap and colorbar plotting features. I used the following code to plot the satellites galaxies and their additional information. Unfortunately, I do not believe some of the information is publicly available yet, so I will be using some stand-in variables. In the following example, ‘z’ is the redshift for each satellite galaxy, ‘vHost’ is the recessional velocity of the host galaxy, ‘satRA’ is the satellite galaxy right ascension, ‘satDEC’ is the satellite galaxy declination,

#calculating the recessional velocity using the redshift, z, in km/s
zVel = (z*3e5- Host)

#Sets the color of the colormap
 
#For continuous color
cMap = 'winter'
    
#For discrete colors
#cMapN = 10 #number of discrete colors(bins)
#cMap = plt.cm.get_cmap('winter',cMapN)
    
      
#plots the location of the satellites

#marker = "$\u25EF$" produced the hollow circles with cmap.
#'vmin' and 'vmax' set the range of values for the colorbar

showPlot = ax1.scatter(satRA, satDEC,s = 225,marker = "$\u25EF$",c = zVel, cmap = cMap,transform=ax1.get_transform('fk5'),vmin = -250,vmax=250)

    
#Colorbar
cbar = fig.colorbar(showPlot)
cbar.set_label(r'$\Delta V (\frac{km}{s})$')

That concludes this portion of the plotting. I tried to be as general as possible, but with the code and examples here, you can generate some informative, simple plots.

-Will

Posted by Will

Welcome!

This summer, I have the great fortune of being a Research Experience for Veteran Undergraduates Fellow, where I will be working under the direction of Dr. Marla Geha.

With the advent of multi-fiber spectroscopy, a multitude of data have become available to astronomers. By observing and identifying emission lines of the resulting spectra, one can determine the abundance of elements that exist within the dwarf galaxy populations. For the purpose of this research, particular attention will be paid to the Hydrogen Alpha emission line. By comparing the emission line to the continuum of the spectra, one can calculate the star formation rates (SFRs) that exist within these populations. Furthermore, an aperture correction can be performed in order to reduce the aperture bias that exists due to the geometry of the spectrometer fibers. This results in a more accurate measurement of SFRs. Using these SFRs, I will attempt to determine if these rates are influenced by the morphological features of the dwarf galaxies.

During this period I hope to learn, or improve on, the following:

  • Statistical Analysis for Astronomical Data
  • Modeling and fitting data
  • Aperture Correction Methods for Fiber Spectroscopy
  • Extracting Data from Stellar Spectroscopy
  • Technical Writing and Communication
  • Scientific Presentations

As the summer continues, and as I progress with this research, I intend to use this website to document my progress and pitfalls. My hope is that the information included here will help other aspiring astronomers in their scientific journies.

Posted by Will