Light 1 Numerical and Data Analysis Answers


In [1]:
%matplotlib notebook
import matplotlib
from astropy.table import Table
import requests

1. Find and plot optical spectra from SDSS of quasars, galaxies, and stars

There are a number of ways to gather spectra from SDSS. Here, we will first query CasJobs to find good spectra classified in each category. Then we will download their spectra using the SDSS SAS API.


In [2]:
from SciServer import CasJobs
from SciServer import Authentication
from io import StringIO

We will define a convenience function to find different types of objects.


In [3]:
def retrieve_objects(sclass='GALAXY'):
    # We define the columns we want and their 
    columns = ('ra', 'dec', 'z', 'class', 'plate', 'fiberid', 'mjd', 'run2d', 'run1d')
    dtypes= ('f8', 'f8', 'f4', 'S', 'i4', 'i4', 'i4', 'S', 'S')

    # Now define the query
    query = """
SELECT TOP 10 {columns}
FROM specObj 
WHERE class = '{sclass}' and zwarning = 0 and snmedian_r > 8
"""
    query = query.format(columns=', '.join(list(columns)), sclass=sclass)

    # Execute the query (requires internet access)
    responseStream = CasJobs.executeQuery(query, "DR12", format="dict")

    # convert result into astropy table
    result = responseStream['Result'][0]
    data = list(map(list, zip(*result['Data'])))
    objects = Table(data, names=columns, dtype=dtypes)
    return(objects)

Let's find some galaxies and show the information.


In [4]:
galaxies = retrieve_objects(sclass='GALAXY')
galaxies.show_in_notebook()


/Users/blanton/anaconda3/lib/python3.6/site-packages/SciServer-1.10.2-py3.6.egg/SciServer/CasJobs.py:116: Warning: In Authentication.getToken: Authentication token is not defined: the user did not log in with the Authentication.login function, or the token has not been stored in the command line argument --ident.
Out[4]:
<Table length=10>
idxradeczclassplatefiberidmjdrun2drun1d
0215.0474225.3454260.0788727GALAXY21285755380026
1214.6031125.1141030.01650546GALAXY21285775380026
2214.7968525.1733840.07925931GALAXY21285785380026
3215.0730125.0259210.06886992GALAXY21285795380026
4214.8701924.6785060.05472199GALAXY21285815380026
5214.7865324.8529440.1652594GALAXY21285825380026
6214.7277924.9748050.1762732GALAXY21285845380026
7214.6950224.9370020.01809064GALAXY21285855380026
8214.8944224.8122440.3367326GALAXY21285865380026
9214.7557824.9437170.01787376GALAXY21285875380026

Now we will download and plot the spectrum of one of the galaxies. This cell will write a file into your local directory called galaxy.fits.


In [5]:
request_template = 'https://dr13.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid={plate}&mjd={mjd}&fiberid={fiberid}'
igalaxy = 5
request = request_template.format(plate=galaxies['plate'][igalaxy], fiberid=galaxies['fiberid'][igalaxy], mjd=galaxies['mjd'][igalaxy])
r = requests.get(request)
fp = open('galaxy.fits', 'wb')
fp.write(r.content)
fp.close()

Now we plot it. Note that "loglam" is the base-10 logarithm of the wavelength in Angstroms, and flux is in $10^{-17}$ erg/cm$^2$/s/A, for this particular spectrum. The broad underlying continuum light is due to stars. The sharp spikes upward are called emission lines and are emitted by ionized gas in between the stars. This particular spectrum has a fair amount of noise.


In [23]:
data = fitsio.read('galaxy.fits')
plt.plot(10.**data['loglam'], data['flux'])
plt.xlabel('$\lambda$ (Angstrom)')
plt.ylabel('$f_\lambda$ ($10^{-17}$ erg/cm$^2$/s/Angstrom)')


Out[23]:
Text(23.0694,0.5,'$f_\\lambda$ ($10^{-17}$ erg/cm$^2$/s/Angstrom)')

We can smooth this spectrum. Here we choose to smooth the spectrum with a Gaussian with $R\sim 1000$ in terms of full-width half maximum. In $\log_{10}\lambda$ coordinates (i.e. what loglam is in) the FWHM is $(1/R/\ln 10)$. Therefore we should convolve with the Gaussian:

$$\frac{1}{\sqrt{2\pi} \sigma} \exp\left[\frac{1}{2} \frac{\Delta{\rm loglam}^2}{\sigma^2}\right]$$

where:

$$\sigma = \frac{1}{2\sqrt{2\ln 2}} \frac{1}{R\ln 10}$$

The code below creates a kernel with the same log wavelength scale as the data, corresponding to $R=1000$.


In [31]:
R = 1000
sigma = (1. / (2. * np.sqrt(2. * np.log(2)))) * (1./(R * np.log(10)))
dpix = data['loglam'][1] - data['loglam'][0]
kernel_npix = (np.int32(np.floor(10. * sigma / dpix)) // 2) * 2 + 1
kernel_xcen = np.float32(kernel_npix // 2)
kernel_x = (np.arange(kernel_npix) - kernel_xcen) * dpix
kernel = np.exp(- 0.5 * (kernel_x / sigma)**2) / (np.sqrt(2 * np.pi) * sigma)

In [32]:
plt.plot(kernel_x, kernel)
plt.plot(kernel_x, kernel, 'o')


Out[32]:
[<matplotlib.lines.Line2D at 0x127234208>]

We can then convolve the flux with this kernel. We end up seeing the features of the galaxy spectrum at higher signal-to-noise, though also at lower resolution. (Note that in a statistical sense the pixels of the smoother spectrum now have strong off-diagonal covariances).


In [34]:
sflux = np.convolve(data['flux'], kernel, mode='same')
plt.plot(10.**data['loglam'], sflux)
plt.xlabel('$\lambda$ (Angstrom)')
plt.ylabel('$f_\lambda$ ($10^{-17}$ erg/cm$^2$/s/Angstrom)')


Out[34]:
Text(0,0.5,'$f_\\lambda$ ($10^{-17}$ erg/cm$^2$/s/Angstrom)')

Now a quasar. In this case the very blue continuum is due to emission from the accretion disk surrounding the massive black hole. There are emission lines in this quasar but they are broad, because of the high random velocities of the emitting gas, which is within a few parsecs of the black hole.


In [7]:
qsos = retrieve_objects(sclass='QSO')
iqso = 5
request = request_template.format(plate=qsos['plate'][iqso], fiberid=qsos['fiberid'][iqso], mjd=qsos['mjd'][iqso])
r = requests.get(request)
fp = open('qso.fits', 'wb')
fp.write(r.content)
fp.close()
data = fitsio.read('qso.fits')
plt.plot(10.**data['loglam'], data['flux'])
plt.xlabel('$\lambda$ (Angstrom)')
plt.ylabel('$f_\lambda$ ($10^{-17}$ erg/cm$^2$/s/Angstrom)')


/Users/blanton/anaconda3/lib/python3.6/site-packages/SciServer-1.10.2-py3.6.egg/SciServer/CasJobs.py:116: Warning: In Authentication.getToken: Authentication token is not defined: the user did not log in with the Authentication.login function, or the token has not been stored in the command line argument --ident.
Out[7]:
Text(0,0.5,'$f_\\lambda$ ($10^{-17}$ erg/cm$^2$/s/Angstrom)')

Finally a star. This star is an A star. The Balmer line absorption sequence can clearly be seen.


In [8]:
stars = retrieve_objects(sclass='STAR')
istar = 5
request = request_template.format(plate=stars['plate'][istar], fiberid=stars['fiberid'][istar], mjd=stars['mjd'][istar])
r = requests.get(request)
fp = open('star.fits', 'wb')
fp.write(r.content)
fp.close()
data = fitsio.read('star.fits')
plt.plot(10.**data['loglam'], data['flux'])
plt.xlabel('$\lambda$ (Angstrom)')
plt.ylabel('$f_\lambda$ ($10^{-17}$ erg/cm$^2$/s/Angstrom)')


/Users/blanton/anaconda3/lib/python3.6/site-packages/SciServer-1.10.2-py3.6.egg/SciServer/CasJobs.py:116: Warning: In Authentication.getToken: Authentication token is not defined: the user did not log in with the Authentication.login function, or the token has not been stored in the command line argument --ident.
Out[8]:
Text(0,0.5,'$f_\\lambda$ ($10^{-17}$ erg/cm$^2$/s/Angstrom)')

Now convert to $f_\nu$.


In [9]:
cspeed = 2.99792e+18 # Angstrom / s
nu = cspeed / 10.**data['loglam'] / 1.e+12 # into terahertz
fnu = data['flux'] * 10.**(2.*data['loglam']) / cspeed * (1.e-17 / 1.e-23 * 1.e+3) # into mJy
plt.plot(nu, fnu)
plt.xlabel(r'$\nu$ (THz)')
plt.ylabel(r'$f_\nu$ (mJy)')


Out[9]:
Text(0,0.5,'$f_\\nu$ (mJy)')

In [ ]: