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
sys.path.insert(0, '../../ppxf')
objname = 'NGC4061'
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
Out[6]:
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 [8]:
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 = .024
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 [11]:
#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[11]:
In [11]: