In [15]:
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 = 'UGC04329'
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 [16]:
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 [17]:
#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[17]:
In [18]:
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[18]:
In [19]:
#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 [20]:
#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 [21]:
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 [26]:
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')
dist = 60.
z_sdss = .014
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 [34]:
#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])
scale_angle = popt[1]
scale_length = scale_angle * 4.848e-6 * dist
scale_length_unc = pcov[1,1] * 4.848e-6 * dist
print 'scale length is', str(np.round(scale_length * 1000, 2)), 'kpc +/-', str(np.round(scale_length_unc * 1000, 2)), 'kpc'
plt.plot(np.linspace(0., 50., 100), galprof(np.linspace(0., 50., 100), popt[0], popt[1]))
plt.axvline(scale_angle, c = 'b')
plt.axvline(scale_angle - pcov[1,1], linestyle = '--', c = 'b')
plt.axvline(scale_angle + pcov[1,1], linestyle = '--', c = 'b')
plt.show()
In [25]:
#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[25]: