SDSS: Data I/O and Binned Plotting¶
In this tutorial we will be loading up some SDSS data (mass, metallicity, and star formation rate) for a bunch of galaxies. We will then plot the Mass metallicity relation, and explore it in a bit more depth.
You can work on this tutorial in anyway you like- You can use this html to copy and paste into your own ipython notebook (or just a python script in a plaintext editor) and then make changes, you can follow along working from scratch, or you can download the .ipynb version of this document from the Tutorials page and work directly in it.
Let's go ahead and get started with our import statements:
import numpy as np
import matplotlib.pyplot as plt
import pyfits as pf
from matplotlib.colors import LogNorm
%matplotlib inline
Great. Now we will want to load the data from the fits files we downloaded. let's define a function that will do this for any fits data table file, and then feed it the three file-names we have.
def load_fits(fname):
#load the fits file into python
data = hdu.data #accesses the data table
return data
sfr_full = load_fits('gal_totsfr_dr7_v5_2.fits')
mass_full = load_fits('totlgm_dr7_v5_2b.fit')
z_full = load_fits('gal_fiboh_dr7_v5_2.fits')
We now have three sets of data. Each of our datasets above has multiple columns being stored in what is known as a "record array" (or rec.array). This is a form of array that allows indexing by attribue, much like a dictionary. We can print out what the column names are by using dot notation on our data sets.
sfr_full.columns
As we can see, we have columns like "Avg", "Entropy", "Median", "Flag", etc. For the purposes of this tutorial, we are interested in the "Avg column." Like most datasets of real data, not all of the galaxies or entries in these files are usable- generally the pipeline that creates the datasets will "flag" bad data in an easily programmatically removable way. Sometimes it is using a "flag" column, and other times is is by selecting an arbitrary and non-physical number to enter as the value (this is how the SDSS data is handled). The values of the flags can usually be found in a readme file.
So, we need to restrict our data to just those that don't have any warning flags. In the cell below, (or in your code), find the indices (locations) for which the following conditions are satisfied, and save them to a variable called "restrictions":
- In the SFR array, the value is > -99
- In the Mass array, the value is not equal to -1
- In the metallicity array, the value is > -99.9
Hint: Do not use a for loop to iterate over the arrays and check the conditions- there is a much faster and more efficient method.
restrictions =
Great, so now we know which indices correspond to the data that's actually useable. Out of curiosity, calculate the number of elements this process actually removed below:
print
Next, we need to actually grab that data (remember, we currently only have the indices). Fill in the rest below
sfr = np.array()
mass =np.array()
z = np.array()
Note: These arrays (because we are using our full arrays from above), are record arrays containing multiple values for each galaxy. We can print out sfr below to see what that means (you can re-run the cell below to ensure you got the same thing I did):
sfr
Essentially, sfr is a an array of the same type as the original, only with the bad entries removed. Using the same indexing, create variables sfrs, masses, and metallicities by indexing your sfr, mass, and z arrays for the "AVG" column.
sfrs =
masses =
metallicities =
print sfrs
Now we are getting somewhere. We have a single array for all the good star formation rates, an array of masses, and an array of metallicities. What is the Mass-Metallicity relation for galaxies? Let's find out. Define a function to plot mass vs metallicity as a 2D histogram with 300 bins (plt has a command for this). You'll want to utilize the "LogNorm" function I imported at the top to get a good color gradient- google "plt 2d histogram LogNorm" if you don't know how to do this already.
def plot_mass_vs_metal(masses,metallicities):
#Plot mass against metalicity as a 2d histogram
plt.colorbar()
plt.title('Mass/Metallicity relation for SDSS Galaxies')
plt.xlabel(r'log Mass [$M_\odot$]')
plt.ylabel(r'log Gas Phase Metallicities')
plt.show()
plot_mass_vs_metal(masses,metallicities)
When you run your function, you should get an image identical to the one I have above. As we can see, the gas phase metallicity of a galaxy (that is, it's composition in terms of elements) has a positive correlation with the log Mass. However, maybe this is not quite a direct relation as we might think. We also pulled in Star Formation Rates. Perhaps those are also correlated? Fill in the functions below to make the histograms for SFR v metallicity and SFR v mass.
def plot_sfr_metal(sfrs,metallicities):
plt.title('SFR/Metallicity relation for SDSS Galaxies')
plt.colorbar()
plt.xlabel(r'log SFR')
plt.ylabel(r'log Gas Phase Metallicities')
plt.show()
def plot_mass_sfr(masses=masses,sfrs=sfrs):
plt.title('Mass/SFR relation for SDSS Galaxies')
plt.colorbar()
plt.xlabel(r'log Mass')
plt.ylabel(r'log SFR')
plt.show()
plot_sfr_metal(sfrs,metallicities)
plot_mass_sfr()
Oh dear. It seems like SFR correlates positively with mass, and metallicity correlates positively with SFR. We need to tease out which of these things are actually correlated, and which only look correlated because they depend on something else which is correlated.
One way we can do this is to take slices of one variable. For example, if metallicity truly does depend on mass, then for all galaxies of a single mass, the correlation should dissapear. We can do such a check for all three of our variables- taking slices of single SFR, metallicity, and mass, and see for which the positive correlation dissapears.
In the space below, write 3 functions, which will bin your data by mass, metallicity, and sfr slices. In theory we would select multiple single slices (choosing a specific value), but in practice our bins will have to have a certain width. To make things easier, I have looked at the data and created a bins array for each function- see if you can figure out how it works (look at the bounds on the graphs above and the behavior of linspace).
def mass_bins(masses):
bins = np.linspace(7,12,10)
binned_masses = []
binned_sfrs = []
binned_z = []
for i in range(len(bins)-1):
mass_indices = #check "where" masses are > left edge of bin and < right edge of bin
mass_indices = np.array(mass_indices)
masses_needed = #index masses for the indices found above
sfr_needed = #ditto for sfrs
z_needed = #ditto for metallicities
binned_masses.append(masses_needed)
binned_sfrs.append(sfr_needed)
binned_z.append(z_needed)
return binned_masses, binned_sfrs, binned_z
def sfr_bins(sfrs):
bins = np.linspace(-2,2,10)
binned_masses = []
binned_sfrs = []
binned_z = []
for i in range(len(bins)-1):
sfr_indices =
to_choose = np.array(sfr_indices)
masses_needed = #same process as above
sfr_needed =
z_needed =
binned_masses.append(masses_needed)
binned_sfrs.append(sfr_needed)
binned_z.append(z_needed)
return binned_masses, binned_sfrs, binned_z
def z_bins(metallicities):
bins = np.linspace(8,9.5,10)
binned_masses = []
binned_sfrs = []
binned_z = []
for i in range(len(bins)-1):
z_indices =
to_choose = np.array(z_indices)
masses_needed = #and once more (but for z)
sfr_needed =
z_needed =
binned_masses.append(masses_needed)
binned_sfrs.append(sfr_needed)
binned_z.append(z_needed)
return binned_masses, binned_sfrs, binned_z
Lastly, let's write some simple functions to plot up all the sliced data: (This time use 100 bins since there is less data per plot)
def plot_mbins(masses=masses):
m,s,z = mass_bins(masses)
for i in range(len(m)):
#2d histogram, think about what to index here
plt.xlabel('Log SFR')
plt.ylabel('Log Metallicity')
plt.colorbar()
plt.figure()
plt.show()
return
def plot_sfrbins(sfrs=sfrs):
m,s,z = sfr_bins(sfrs)
for i in range(len(m)):
#2d histogram, think about what to index here
plt.xlabel('Log Mass')
plt.ylabel('Log Metallicity')
plt.colorbar()
plt.figure()
plt.show()
return
def plot_zbins(metallicities=metallicities):
m,s,z = z_bins(metallicities)
for i in range(len(m)):
#2d histogram, think about what to index here
plt.xlabel('Log Mass')
plt.ylabel('Log SFR')
plt.colorbar()
plt.figure()
plt.show()
return
As a final step, lets run our functions below, and see what we get:
plot_mbins()
plot_sfrbins()
plot_zbins()
As we can see, the relationship between SFR and metallicity disappears at a single mass slice (the first ten plots). Thus, SFR and metallicity are not truly correlated, but only appear so when all masses are included because each depends on mass (which can be seen in the second and third set of plots, where for single slice in metallicity, the SFR/mass relation still exists, and for single slices in SFR, the mass/metallicity relation still exists.
Congrats! You made it to the end of the tutorial. I hope you enjoyed it, practiced a little python, and learned something about galaxy properties. As always, feel free to contact me (post an issue on the github http://github.com/prappleizer/prappleizer.github.io ) if anything was messed up, confusing, or poorly explained.