In [1]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from astropy.io import fits
import astropy.table as table
import matplotlib as mpl
import matplotlib.image as mpimg
import os
sys.path.insert(0, '../../ppxf')
objname = 'NGC2562'
im=fits.open(objname + '.msobj.fits')
#plt.imshow(im[0].data, cmap='bone')
l = table.Table.read('../arc_disp.txt', format='ascii')
fibers = table.Table.read('../../fibers.txt', format='ascii')
#get the target coordinates and print them
f = open(objname + '.list')
targ_raw_im_name = f.readline()
h = fits.open(targ_raw_im_name)[0].header
#print h['TARGRA']
#print h['TARGDEC']
ralist = [float(i) for i in str(h['TARGRA']).split(':')]
declist = [float(i) for i in str(h['TARGDEC']).split(':')]
for i in range(0, 3):
ralist[i] = ralist[i] * [360./24., 360./24./60., 360./24./60./60.][i]
declist[i] = declist[i] * [1., 1./60., 1./60./60.][i]
print 'coords:', sum(ralist), sum(declist)
In [2]:
def circles(x, y, s, c='b', ax=None, vmin=None, vmax=None, **kwargs):
"""
Make a scatter of circles plot of x vs y, where x and y are sequence
like objects of the same lengths. The size of circles are in data scale.
Parameters
----------
x,y : scalar or array_like, shape (n, )
Input data
s : scalar or array_like, shape (n, )
Radius of circle in data scale (ie. in data unit)
c : color or sequence of color, optional, default : 'b'
`c` can be a single color format string, or a sequence of color
specifications of length `N`, or a sequence of `N` numbers to be
mapped to colors using the `cmap` and `norm` specified via kwargs.
Note that `c` should not be a single numeric RGB or
RGBA sequence because that is indistinguishable from an array of
values to be colormapped. `c` can be a 2-D array in which the
rows are RGB or RGBA, however.
ax : Axes object, optional, default: None
Parent axes of the plot. It uses gca() if not specified.
vmin, vmax : scalar, optional, default: None
`vmin` and `vmax` are used in conjunction with `norm` to normalize
luminance data. If either are `None`, the min and max of the
color array is used. (Note if you pass a `norm` instance, your
settings for `vmin` and `vmax` will be ignored.)
Returns
-------
paths : `~matplotlib.collections.PathCollection`
Other parameters
----------------
kwargs : `~matplotlib.collections.Collection` properties
eg. alpha, edgecolors, facecolors, linewidths, linestyles, norm, cmap
Examples
--------
a = np.arange(11)
circles(a, a, a*0.2, c=a, alpha=0.5, edgecolor='none')
"""
from matplotlib.patches import Circle
from matplotlib.collections import PatchCollection
import pylab as plt
#import matplotlib.colors as colors
if ax is None:
ax = plt.gca()
if isinstance(c,basestring):
color = c # ie. use colors.colorConverter.to_rgba_array(c)
else:
color = None # use cmap, norm after collection is created
kwargs.update(color=color)
if isinstance(x, (int, long, float)):
patches = [Circle((x, y), s),]
elif isinstance(s, (int, long, float)):
patches = [Circle((x_,y_), s) for x_,y_ in zip(x,y)]
else:
patches = [Circle((x_,y_), s_) for x_,y_,s_ in zip(x,y,s)]
collection = PatchCollection(patches, **kwargs)
if color is None:
collection.set_array(np.asarray(c))
if vmin is not None or vmax is not None:
collection.set_clim(vmin, vmax)
ax.add_collection(collection)
return collection
In [3]:
#test fiber map
fibersize=4.687
fig = plt.figure(figsize=(6,6))
for row in fibers:
plt.scatter(row['ra'], row['dec'], s=200, c='y')
plt.text(row['ra']-2, row['dec']-1.5, s=row['fiber'], color='k')
plt.xlabel(r'RA (arcsec)')
plt.ylabel(r'dec (arcsec)')
plt.title('SparsePak fiber spectrograph rows')
Out[3]:
In [4]:
fig = plt.figure(figsize=(6,6), dpi=200)
for row in fibers:
if row['sky']!=1:
plt.scatter(row['ra'], row['dec'], s=500, c='y')
plt.text(row['ra']-.85, row['dec']-.85, s=row['row'], color='k')
plt.xlabel(r'RA offset from center (arcsec)')
plt.ylabel(r'Dec offset from center (arcsec)')
plt.title('SparsePak fiber data rows')
Out[4]:
In [5]:
#now put the fibers into radius bins
fibers = table.Table.read('../../fibers.txt', format='ascii')
radius = table.Column(name = 'r', data = np.sqrt( (np.asarray(fibers['ra']))**2. + (np.asarray(fibers['dec']))**2. ) )
#print 'radius:', len(radius)
#print 'fibers:', len(fibers)
if 'radius' not in fibers.colnames:
fibers.add_column(radius)
fig = plt.figure(figsize=(6,6), dpi=200)
for row in fibers:
if row['sky']!=1:
plt.scatter(row['ra'], row['dec'], s=500, c='y')
plt.text(row['ra']-.85, row['dec']-.85, s=row['fiber'], color='k')
#print row['ra'], row['dec'], row['row']
circles(np.zeros(len(np.unique(radius))), np.zeros(len(np.unique(radius))), np.unique(radius)[::-1], c='w', edgecolor = 'k', alpha = 0.9, linestyle = '--', lw = 1, zorder = 0)
plt.xlabel(r'RA offset from center (arcsec)')
plt.ylabel(r'Dec offset from center (arcsec)')
plt.title('SparsePak fiber rows')
fibers.sort(['r'])
#fibers.pprint(max_lines=-1, max_width=-1)
In [6]:
#now figure out radius groups
fig = plt.figure(figsize=(16,16), dpi=200)
from sklearn.cluster import MeanShift, estimate_bandwidth
R = np.array(zip(fibers['r'][fibers['sky']!=1], np.zeros(len(fibers['r'][fibers['sky']!=1]))))
bandwidth = estimate_bandwidth(R, quantile=.08)
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(R)
labels = ms.labels_
cluster_centers = ms.cluster_centers_
labels_unique = np.unique(labels)
n_clusters_ = len(labels_unique)
#now match fibers to a radius group by iterating
group = -9999*np.ones(len(fibers['r'])).astype(int)
r_group_dict = {}
for k in range(n_clusters_):
my_members = labels == k
#print "cluster {0}: {1}".format(k, R[my_members, 0])
plt.axhline(np.mean(R[my_members, 0]))
#print k, my_members
for i, fiber in enumerate(my_members):
if fiber==True:
group[i] = k
r_group_dict[k] = np.mean(R[my_members, 0])
#print group
print r_group_dict
r_group = table.Column(name = 'r group', data = group)
if 'r group' not in fibers.colnames:
fibers.add_column(r_group)
for i, row in enumerate(fibers):
if row['sky']!=1:
plt.scatter(i, row['r'])
plt.annotate(str(row['fiber']) + '\n' + str(row['row']), (i, row['r']), xytext=(i, row['r']+4*(i%2)-2.), arrowprops=dict(arrowstyle='-|>'))
#print row['ra'], row['dec'], row['row']
#plt.yscale('log')
#plt.ylim([-5., 60.])
plt.title('radial bins')
plt.annotate('fiber \n data', (8, 35), xytext = (10, 40), fontsize = 16, arrowprops=dict(arrowstyle='-|>'))
print fibers
#fibers now includes a radius-group designation for each fiber, so we can radially bin data. This method should actually be fairly extensible to other IFUs
#we also have a dictionary r_group_dict, which associates each group number with a mean radius
In [7]:
fig = plt.figure(figsize=(12,6), dpi=200)
im=fits.open(objname + '.msobj.fits')
l = table.Table.read('../arc_disp.txt', format='ascii')
ifu=im[0].data
for i in range(0, len(ifu)):
plt.plot(l['Wavelength'], ifu[i])
plt.xlabel('observed wavelength ($\AA$)')
plt.ylabel('counts')
plt.title('All fibers')
plt.show()
In [20]:
fig = plt.figure(figsize=(12,12), dpi=200)
ax = fig.add_subplot(111)
im=fits.open(objname + '.msobj.fits')
l = table.Table.read('../arc_disp.txt', format='ascii')
z_sdss = .017
l_corr = 1 + z_sdss
fiberlist = [47, 42, 22, 41, 27, 60, 38, 56, 39, 36, 0, 21, 43, 46, 40, 56, 37, 18, 19, 24, 20, 29, 72, 64, 58, 74, 57, 50]
fiberlist = [42, 22, 41, 27, 60, 38, 56, 39, 36, 0, 21, 43, 46, 40, 56, 37, 18, 19, 24, 20, 29, 72, 64, 58, 74, 57, 50]
fiberlist = [42, 22, 41, 27, 60, 38, 56, 39, 36, 0, 21, 43, 46, 40]
labely = np.linspace(.1, .9, len(fiberlist))
#print labely
maxes = [ifu[i].max() for i in fiberlist]
maxarray = np.array(zip(fiberlist, maxes))
maxarray = maxarray[maxarray[:,1].argsort()]
#print maxarray
ifu=im[0].data
ax.plot(l['Wavelength'] / l_corr, ifu[i])
#.8*maxes[fiberlist.index(i)]/np.max(maxes)
for i in range(0, len(ifu)):
if i in fiberlist:
spectrum, = ax.plot(l['Wavelength'] / l_corr, ifu[i])
ax.annotate(str(i), (l['Wavelength'][np.argmax(ifu[i])] / l_corr, ifu[i][np.argmax(ifu[i])]), xytext=(.92, labely[np.where(maxarray[:,0]==i)] ), textcoords = 'figure fraction', fontsize = 12, arrowprops=dict(facecolor = 'black', arrowstyle='-|>'), color = spectrum.get_color())
ax.set_xlim(ax.get_xlim()[0], ax.get_xlim()[1] + 100)
ax.grid(linestyle = '--')
ax.set_xlabel('emitted wavelength ($\AA$)')
ax.set_ylabel('counts')
#ax.set_title('Fiber spectra (labeled by data row)')
ax2 = ax.twiny()
ax2.set_xlim(ax.get_xlim()[0] * l_corr, (ax.get_xlim()[1] + 100) * l_corr)
ax2.set_xlabel('observed wavelength ($\AA$)')
ax2.grid(color = 'b', linestyle = ':')
plt.subplots_adjust(hspace=0.1)
In [87]:
#try radial binning
#first plot total flux against by radius
def galprof(r, i0, rd):
return i0*np.exp(-r/rd)
total_light = np.sum(ifu, axis=1)
fibers.sort(['row'])
rlist = []
ilist = []
for row in fibers:
if row['sky']!=1:
rlist.append(r_group_dict[row['r group']])
ilist.append(total_light[row['row']])
plt.errorbar(r_group_dict[row['r group']], total_light[row['row']], marker = 'x', c = 'k')
plt.xlabel('r (arcsec)')
plt.ylabel('counts')
plt.title('light gathered by radius')
rlist = np.asarray(rlist)
ilist = np.asarray(ilist)
popt, pcov = sp.optimize.curve_fit(galprof, rlist, ilist, p0 = [3.7e7, 5])
print popt
plt.plot(np.linspace(0., 50., 100), galprof(np.linspace(0., 50., 100), popt[0], popt[1]))
plt.axvline(popt[1])
plt.show()
In [28]:
#try mapping the log of total light onto an IFU fiber diagram
import matplotlib
fig= plt.figure(figsize=(12,12), dpi=400)
#you can get custom-sized images at http://skyserver.sdss3.org/dr10/en/tools/chart/image.aspx
#explore at http://skyserver.sdss3.org/dr10/en/tools/explore/default.aspx
im=mpimg.imread(objname + '.png')
x = y = np.linspace(-40., 40., shape(im)[0])
X, Y = np.meshgrid(x, y)
galim = plt.imshow(im, extent = [-40, 40, -40, 40], origin = 'lower', interpolation = 'nearest', zorder = 0)
total_light = np.sum(ifu, axis=1)
fibers.sort(['row'])
plt.contour(X, Y, im.sum(axis = -1), zorder = 1, cmap = 'Greens')
#print fibers
cir = circles(fibers['ra'][fibers['sky'] != 1], fibers['dec'][fibers['sky'] != 1], s = fibersize/2, c = total_light, edgecolor = 'w', alpha = 0.8, norm=matplotlib.colors.LogNorm(), cmap = 'hot', zorder = 2)
for row in fibers:
if row['sky']!=1:
plt.text(row['ra']-.85, row['dec']-1.25, s=str(row['fiber']) + '\n' + str(row['row']), color = 'g', size = 14, zorder = 3)
fig.colorbar(cir, shrink=0.8)
plt.xlim([-40, 40])
plt.ylim([-40, 40])
plt.xlabel('RA offset (arcsec)', fontsize = 16)
plt.ylabel('Dec offset (arcsec)', fontsize = 16)
plt.title(objname, fontsize=18)
#this is the correct way to do this!
Out[28]:
In [38]:
from ppxf import ppxf
import ppxf_util as util
import glob
from scipy import ndimage
from time import clock
directory = 'spectra/'
from ppxf import ppxf
import ppxf_util as util
In [42]:
#first log-rebin the spectrum
lamRange1 = np.array([l['Wavelength'][0] / l_corr, l['Wavelength'][-1] / l_corr]) #lower and upper limits of the wavelength
FWHM_gal = 1.4 * 4
galaxy, logLam1, velscale = util.log_rebin(lamRange1, ifu[0])
galaxy = galaxy/np.median(galaxy)
noise = galaxy*0 + 0.0049
# Read the list of filenames from the Single Stellar Population library
# by Vazdekis (1999, ApJ, 513, 224). A subset of the library is included
# for this example with permission. See http://purl.org/cappellari/idl
# for suggestions of more up-to-date stellar libraries.
vazdekis = glob.glob(directory + 'Rbi1.30z*.fits')
FWHM_tem = 1.8 # Vazdekis spectra have a resolution FWHM of 1.8A.
# Extract the wavelength range and logarithmically rebin one spectrum
# to the same velocity scale of the SAURON galaxy spectrum, to determine
# the size needed for the array which will contain the template spectra.
hdu = fits.open(vazdekis[0])
ssp = hdu[0].data
h2 = hdu[0].header
lamRange2 = h2['CRVAL1'] + np.array([0.,h2['CDELT1']*(h2['NAXIS1']-1)])
sspNew, logLam2, velscale = util.log_rebin(lamRange2, ssp, velscale=velscale)
templates = np.empty((sspNew.size,len(vazdekis)))
# Convolve the whole Vazdekis library of spectral templates
# with the quadratic difference between the SAURON and the
# Vazdekis instrumental resolution. Logarithmically rebin
# and store each template as a column in the array TEMPLATES.
# Quadratic sigma difference in pixels Vazdekis --> SAURON
# The formula below is rigorously valid if the shapes of the
# instrumental spectral profiles are well approximated by Gaussians.
FWHM_dif = np.sqrt(FWHM_gal**2 - FWHM_tem**2)
sigma = FWHM_dif/2.355/h2['CDELT1'] # Sigma difference in pixels
for j in range(len(vazdekis)):
hdu = fits.open(vazdekis[j])
ssp = hdu[0].data
ssp = ndimage.gaussian_filter1d(ssp,sigma)
sspNew, logLam2, velscale = util.log_rebin(lamRange2, ssp, velscale=velscale)
templates[:,j] = sspNew/np.median(sspNew) # Normalizes templates
# The galaxy and the template spectra do not have the same starting wavelength.
# For this reason an extra velocity shift DV has to be applied to the template
# to fit the galaxy spectrum. We remove this artificial shift by using the
# keyword VSYST in the call to PPXF below, so that all velocities are
# measured with respect to DV. This assume the redshift is negligible.
# In the case of a high-redshift galaxy one should de-redshift its
# wavelength to the rest frame before using the line below (see above).
c = 299792.458
dv = (logLam2[0]-logLam1[0])*c # km/s
vel = 450. # Initial estimate of the galaxy velocity in km/s
goodPixels = util.determine_goodpixels(logLam1,lamRange2,vel)
# Here the actual fit starts. The best fit is plotted on the screen.
# Gas emission lines are excluded from the pPXF fit using the GOODPIXELS keyword.
start = [vel, 180.] # (km/s), starting guess for [V,sigma]
t = clock()
pp = ppxf(templates, galaxy, noise, velscale, start, goodpixels=goodPixels, plot=True, moments=4, degree=4, vsyst=dv)
print("Formal errors:")
print(" dV dsigma dh3 dh4")
print("".join("%8.2g" % f for f in pp.error*np.sqrt(pp.chi2)))
print('Elapsed time in PPXF: %.2f s' % (clock() - t))
In [ ]: