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 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 = #opens the FITS, assigns it to 'hdu'
print( #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.

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')
ax1.tick_params(direction = 'in')

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 ='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)

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.