In [1]:
%pylab inline
In [2]:
from mpl_toolkits.basemap import Basemap
from matplotlib.collections import PolyCollection
In [3]:
import astropy.units as u
from astropy.coordinates import SkyCoord
In [4]:
import healpy as hp
print(hp.version.__version__)
In [5]:
import bossdata.meta
print(bossdata.__version__)
Ignore expected (harmless) warnings.
In [6]:
import warnings, matplotlib.cbook
warnings.filterwarnings('ignore', category=matplotlib.cbook.mplDeprecation)
Get list of ra, dec, z, and warning flags for quasar observations.
In [7]:
quasar_catalog = bossdata.meta.Database(quasar_catalog=True)
In [8]:
quasar_table = quasar_catalog.select_all(what='RA,DEC,Z_VI,BAL_FLAG_VI,ZWARNING', max_rows=0)
print('Found {0} total quasars'.format(len(quasar_table)))
In [9]:
bal_flagged = quasar_table[quasar_table['BAL_FLAG_VI'] != 0]
print('Found {0} with BAL identified from visual inspection'.format(len(bal_flagged)))
In [10]:
zwarning_flagged = quasar_table[quasar_table['ZWARNING'] != 0]
print('Found {0} with ZWARNING from pipepline'.format(len(zwarning_flagged)))
In [11]:
hiz_quasars = quasar_table[(quasar_table['Z_VI'] > 2.1) & (quasar_table['Z_VI'] < 3.5)]
print('Found {0} in high redshift sample (2.1 < z < 3.5)'.format(len(hiz_quasars)))
In [12]:
np.min(quasar_table['Z_VI']), np.max(quasar_table['Z_VI'])
Out[12]:
Plot redshift distribution of quasar observations
In [13]:
plt.figure(figsize=(8,6))
dr12_survey_area = 10400.0 # square degrees
wgt = 1.0/dr12_survey_area
zbins = np.linspace(0,6.5,66)
plt.hist(quasar_table['Z_VI'], weights=wgt*np.ones(len(quasar_table)), label='DR12Q', bins=zbins, histtype='step')
plt.hist(bal_flagged['Z_VI'], weights=wgt*np.ones(len(bal_flagged)), label='BAL_FLAG_VI != 0', bins=zbins, histtype='step')
plt.hist(zwarning_flagged['Z_VI'], weights=wgt*np.ones(len(zwarning_flagged)), label='ZWARNING != 0', bins=zbins, histtype='step')
plt.xlabel(r'Redshift $z$')
plt.ylabel(r'Number of quasars per $\Delta z = %.1f$ per sq degree' % (zbins[1]-zbins[0]))
plt.axvline(2.1, c='k', ls='--', lw=2)
plt.axvline(3.5, c='k', ls='--', lw=2)
plt.xlim(0,6.5)
plt.grid()
plt.legend()
plt.show()
Plot an "all sky map" using healpix binning. We use the Eckert IV projection as it is the best area preserving projection (see this link for discussion of "best" map projections).
In [14]:
def plot_sky(ra, dec, data=None, nside=16, label='', projection='eck4', cmap=plt.get_cmap('jet'), norm=None,
hide_galactic_plane=False):
# get pixel area in degrees
pixel_area = hp.pixelfunc.nside2pixarea(nside, degrees=True)
# find healpixels associated with input vectors
pixels = hp.ang2pix(nside, 0.5*np.pi-np.radians(dec), np.radians(ra))
# find unique pixels
unique_pixels = np.unique(pixels)
# count number of points in each pixel
bincounts = np.bincount(pixels)
# if no data provided, show counts per sq degree
# otherwise, show mean per pixel
if data is None:
values = bincounts[unique_pixels]/pixel_area
else:
weighted_counts = np.bincount(pixels, weights=data)
values = weighted_counts[unique_pixels]/bincounts[unique_pixels]
# find pixel boundaries
corners = hp.boundaries(nside, unique_pixels, step=1)
corner_theta, corner_phi = hp.vec2ang(corners.transpose(0,2,1))
corner_ra, corner_dec = np.degrees(corner_phi), np.degrees(np.pi/2-corner_theta)
# set up basemap
m = Basemap(projection=projection, lon_0=90, resolution='l', celestial=True)
m.drawmeridians(np.arange(0, 360, 30), labels=[0,0,1,0], labelstyle='+/-')
m.drawparallels(np.arange(-90, 90, 15), labels=[1,0,0,0], labelstyle='+/-')
m.drawmapboundary()
# convert sky coords to map coords
x,y = m(corner_ra, corner_dec)
# regroup into pixel corners
verts = np.array([x.reshape(-1,4), y.reshape(-1,4)]).transpose(1,2,0)
# Make the collection and add it to the plot.
coll = PolyCollection(verts, array=values, cmap=cmap, norm=norm, edgecolors='none')
plt.gca().add_collection(coll)
plt.gca().autoscale_view()
if not hide_galactic_plane:
# generate vector in galactic coordinates and convert to equatorial coordinates
galactic_l = np.linspace(0, 2*np.pi, 1000)
galactic_plane = SkyCoord(l=galactic_l*u.radian, b=np.zeros_like(galactic_l)*u.radian, frame='galactic').fk5
# project to map coordinates
galactic_x, galactic_y = m(galactic_plane.ra.degree, galactic_plane.dec.degree)
m.scatter(galactic_x, galactic_y, marker='.', s=2, c='k')
# Add a colorbar for the PolyCollection
plt.colorbar(coll, orientation='horizontal', pad=0.01, aspect=40, label=label)
return m
Plot the number of quasars per square degree:
In [15]:
plt.figure(figsize=(12,9))
plot_sky(quasar_table['RA'].data, quasar_table['DEC'].data, label='Number of quasars per square degree')
plt.show()
Plot the number of high redshift quasars per square degree:
In [16]:
plt.figure(figsize=(12,9))
plot_sky(hiz_quasars['RA'].data, hiz_quasars['DEC'].data, label='Number of quasars per square degree')
plt.show()
We can show the mean value of a quantity on the sky as well:
In [17]:
plt.figure(figsize=(12,9))
plot_sky(quasar_table['RA'].data, quasar_table['DEC'].data, data=quasar_table['Z_VI'].data,
label='Mean redshift', nside=16,
norm=mpl.colors.Normalize(vmin=1, vmax=3))
plt.show()
Using smaller bins (adjust using the nside parameter), we can see finer details:
In [18]:
plt.figure(figsize=(12,9))
plot_sky(quasar_table['RA'].data, quasar_table['DEC'].data, data=quasar_table['Z_VI'].data,
label='Mean redshift', nside=32,
norm=mpl.colors.Normalize(vmin=1, vmax=3))
plt.show()
In [ ]: