Lecturer: Leo Singer
Jupyter Notebook Authors: Leo Singer, Dave Cook, Shreya Anand & Cameron Hummels
This is a Jupyter notebook lesson taken from the GROWTH Summer School 2019. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-school-2019.html
Learn how to use LIGO/Virgo localizations and match with galaxies
See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with pip install <module>
. The external astromatic packages are easiest installed using package managers (e.g., rpm
, apt-get
).
None
First, some imports: Numpy, Matplotlib, Healpy, and parts of Astropy.
import astropy.utils.data
from matplotlib import pyplot as plt
import numpy as np
import healpy as hp
Here are some extra imports for the galaxy cross matching:
from astropy.table import Table, vstack, hstack, Column
import astropy.units as u
from astropy.coordinates import SkyCoord
import ligo.skymap.plot
from scipy.stats import norm
import scipy.stats
And configure Matplotlib to send plot output directly to the notebook:
%matplotlib inline
This section on using HEALPix localization files is adapted from the LIGO/Virgo Public Alerts User Guide.
Let's start by downloading a sample localization file from the User Guide. We could do this on the command line using curl
:
$ curl -O https://emfollow.docs.ligo.org/userguide/_static/bayestar.fits.gz
But after all, this is a Python lesson, so let's download the file using the handy astropy.utils.data.download_file
function from Astropy.
url = 'https://emfollow.docs.ligo.org/userguide/_static/bayestar.fits.gz'
filename = astropy.utils.data.download_file(url)
Next, let's read in the HEALPix data using Healpy. Note that by default, Healpy only reads the first column, which provides the 2D probability distribution on the sky.
prob = hp.read_map(filename)
NSIDE = 2048 ORDERING = NESTED in fits file INDXSCHM = IMPLICIT Ordering converted to RING
To get a quick look at a HEALPix data set, you can use the hp.mollview
function:
hp.mollview(prob)
What actually is stored in prob
?
prob
array([2.70726059e-66, 1.27374324e-66, 2.62611513e-67, ..., 2.04700874e-40, 1.05781210e-35, 4.44174764e-31])
It's a one-dimensional array! Yet it represents in 2D image. How does that work? HEALPix is a way to index equal-area regions on the unit sphere using integers.
To decode HEALPix indices, you need to know the resolution of the map, which is described by a parameter called nside
. nside
is the number of subdivisions of 12 base HEALPix tiles, so the relation between the length of a HEALPix array, npix
, and its resolution, nside
, is
The functions hp.npix2nside
and hp.nside2npix
convert between length and resolution.
npix = len(prob)
npix
50331648
nside = hp.npix2nside(npix)
nside
2048
The function hp.pix2ang
allow us to convert from (ra, dec) and HEALPix pixel index.
Note: by default, these functions return 'physics' spherical coordinates $(\theta, \phi)$ in radians, but you can switch to 'astronomy' spherical coordinates in degrees by passing the keyword argument lonlat=True
.
Let's look up the right asce3nsion and declination of pixel 123.
ipix = 123
ra, dec = hp.pix2ang(nside, ipix, lonlat=True)
ra, dec
(129.375, 89.81725848475484)
The function hp.ang2pix
does the opposite. Let's find the pixel that contains the point RA=194.95, Dec=27.98.
ra = 194.95
dec = 27.98
hp.ang2pix(nside, ra, dec, lonlat=True)
13361492
What is the most probable sky location? Just find the pixel with the maximum value, and then find its right ascension and declination.
ipix_max = np.argmax(prob)
ipix_max
32883013
hp.pix2ang(nside, ipix_max, lonlat=True)
(194.30419921875, -17.856895095545454)
Finding the most probable sky location within a HEALPix map involves knowing which pixels correspond to a certain probability contour (say, 90%). We can gain insight into how these probability contours are calculated using scipy.stats. Scipy provides a "t" distribution class that we can use to get values from the "t" statistic probability density function (PDF). As a start, we plot the PDF for a "t" statistic with 3 degrees of freedom:
t_dist = scipy.stats.t(3)
t_values = np.linspace(-4, 4, 1000)
plt.plot(t_values, t_dist.pdf(t_values))
plt.xlabel('t value')
plt.ylabel('probability for t value')
plt.show()
The t distribution object t_dist can also give us the cumulative distribution function (CDF). The CDF gives the area under the curve of the PDF at and to the left of the given t value:
plt.plot(t_values, t_dist.cdf(t_values))
plt.xlabel('t value')
plt.ylabel('probability for t value <= t')
plt.title('CDF for t distribution')
plt.show()
Say I have a t value x drawn from a t distribution. The PDF gives the probability for given values of x. Because it is a probability density, the sum of the probabilities of all possible values for x: ∞<x<∞ must be 1. Therefore the total area under the PDF curve is 1, and the maximum value of the CDF is 1.
The CDF gives us the area under the PDF curve at and to the left of a given t value x. Therefore it is the probability that we will observe a value x<=t if we sample a value x from a t distribution.
Let's show relationship of PDF and CDF for three example t values.
example_values = (-1.5, 0, 1.5)
pdf_values = t_dist.pdf(t_values)
cdf_values = t_dist.cdf(t_values)
fill_color = (0, 0, 0, 0.1) # Light gray in RGBA format.
line_color = (0, 0, 0, 0.5) # Medium gray in RGBA format.
fig, axes = plt.subplots(2, len(example_values), figsize=(10, 6))
for i, x in enumerate(example_values):
cdf_ax, pdf_ax = axes[:, i]
cdf_ax.plot(t_values, cdf_values)
pdf_ax.plot(t_values, pdf_values)
# Fill area at and to the left of x.
pdf_ax.fill_between(t_values, pdf_values,
where=t_values <= x,
color=fill_color)
pd = t_dist.pdf(x) # Probability density at this value.
# Line showing position of x on x-axis of PDF plot.
pdf_ax.plot([x, x],
[0, pd], color=line_color)
cd = t_dist.cdf(x) # Cumulative distribution value for this x.
# Lines showing x and CDF value on CDF plot.
x_ax_min = cdf_ax.axis()[0] # x position of y axis on plot.
cdf_ax.plot([x, x, x_ax_min],
[0, cd, cd], color=line_color)
cdf_ax.set_title('x = {:.1f}, area = {:.2f}'.format(x, cd))
# Hide top and right axis lines and ticks to reduce clutter.
for ax in (cdf_ax, pdf_ax):
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
For example, say I have drawn a t value x at random from a t distribution. The probability that x<=1.5 is (i.e., >0.9253):
t_dist.cdf(1.5)
0.8847080673775886
The total area under the PDF is 1, and the maximum value for the CDF is 1. Therefore the area of the PDF to the right of 1.5 must be (i.e., >0.0746):
1 - t_dist.cdf(1.5)
0.11529193262241144
This is the probability that our t value x will be >1.5. In general, when we sample a value x at random from a t distribution, the probability that x>q is given by:
ℙ(x>q)=1−CDF(q), where CDF is the cumulative distribution function for a t value. We can apply the same methodology to HEALpix pixel probabilities in LIGO/VIRGO localization maps.
First, let's get our galaxy catalog that we will later match to the 3D localization of GW170817.
For this Section we will use a galaxy catalog from the CLU project (Census of the Local Universe; paper: https://ui.adsabs.harvard.edu/abs/2017arXiv171005016C/abstract). However, we will only use those galaxies that are publically availble and in NED (NASA/IPAC Extragalactic Database: https://ned.ipac.caltech.edu/). This catalog has already been prepared for you.
load CLU catalog into astropy table.
clu=Table.read('data/CLU_NEDonly.fits')
nclu=np.size(clu)
Add probability columns to the galaxy catalog: probability density and p-value per volume and per area.
probdencol=Column(np.zeros(nclu,dtype='f4'),name='dP_dV')
probcol=Column(np.zeros(nclu,dtype='f4'),name='P')
probdenAcol=Column(np.zeros(nclu,dtype='f4'),name='dP_dA')
probAcol=Column(np.zeros(nclu,dtype='f4'),name='P_A')
clu.add_columns([probdencol,probcol,probdenAcol,probAcol])
Familiarize yourself with the catalog
print the columns in the catalog that will be used in the cross-match.
clu['NEDname','RA','DEC','TYPE_NED','DISTMPC','Z','ZERR','MODELMAG_R','K_M_K20FE','K_MSIG_K20FE','W1MPRO','W1SIGMPRO']
NEDname | RA | DEC | TYPE_NED | DISTMPC | Z | ZERR | MODELMAG_R | K_M_K20FE | K_MSIG_K20FE | W1MPRO | W1SIGMPRO |
---|---|---|---|---|---|---|---|---|---|---|---|
bytes30 | float64 | float64 | bytes3 | float32 | float64 | float64 | float32 | float64 | float64 | float64 | float64 |
2dFGRS S805Z417 | 0.00162 | -56.14106 | G | 43.25686 | 0.010381646454334259 | 0.000297 | nan | 1e+20 | 1e+20 | 15.983 | 0.046 |
2MASS J00000158-3930463 | 0.0068 | -39.51309 | G | 46.027866 | 0.011046688072383404 | 0.000119 | nan | 1e+20 | 1e+20 | 1e+20 | 1e+20 |
KUG 2357+156 | 0.00905 | 15.88188 | G | 85.72995 | 0.020575188100337982 | 3.3e-05 | 14.917557 | 12.31 | 0.067 | 12.178 | 0.023 |
2dFGRS S199Z178 | 0.00996 | -27.74194 | G | 114.32659 | 0.02743838168680668 | 0.000297 | nan | 1e+20 | 1e+20 | 15.268 | 0.04 |
GALEXASC J000003.23-333033.4 | 0.01362 | -33.50967 | G | 98.505714 | 0.023641372099518776 | 0.000297 | nan | 1e+20 | 1e+20 | 16.651 | 0.083 |
GALEXASC J000005.52-011255.7 | 0.02304 | -1.2155 | UvS | 130.27467 | 0.031265921890735626 | 1e+20 | 17.157835 | 1e+20 | 1e+20 | 15.292 | 0.04 |
USGC U001 | 0.02333 | 26.19028 | GTr | 106.429 | 0.025542961433529854 | 1e+20 | 22.110594 | 1e+20 | 1e+20 | 17.719 | 0.213 |
DUKST 349-064 | 0.02492 | -35.65822 | G | 121.63314 | 0.02919195219874382 | 0.000213 | nan | 1e+20 | 1e+20 | 14.699 | 0.046 |
UGC 12890 | 0.02931 | 8.27918 | G | 165.74657 | 0.0397791750729084 | -999.0 | nan | 10.779 | 0.059 | 11.521 | 0.023 |
MCG -01-01-016 | 0.036 | -6.374 | G | 93.52904 | 0.022446969524025917 | 2.7e-05 | 14.141832 | 11.24 | 0.059 | 11.946 | 0.023 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
GALEXASC J234305.34-552311.9 | 355.77191 | -55.38668 | Vis | 85.65714 | 0.020557714626193047 | -999.0 | nan | 1e+20 | 1e+20 | 16.163 | 0.048 |
GALEXASC J234344.33-545403.7 | 355.93484 | -54.90119 | Vis | 42.82857 | 0.010278857313096523 | -999.0 | nan | 1e+20 | 1e+20 | 15.047 | 0.032 |
SDSS J234403.01+003122.6 | 356.01255 | 0.52295 | G | 98.111694 | 0.02354680746793747 | 1e+20 | 21.51647 | 1e+20 | 1e+20 | 1e+20 | 1e+20 |
GALEXASC J234414.08-364635.4 | 356.05836 | -36.77655 | G | 26.48519 | 0.006356445141136646 | 1.3e-05 | nan | 1e+20 | 1e+20 | 15.055 | 0.036 |
2MASX J23454008+0144274 | 356.41704 | 1.74092 | G | 119.54739 | 0.028691373765468597 | 1e+20 | 15.659925 | 14.084 | 0.168 | 14.073 | 0.03 |
AGC 331333 | 356.75349 | 27.2631 | G | 128.07455 | 0.030737893655896187 | 0.00032 | nan | 1e+20 | 1e+20 | 1e+20 | 1e+20 |
GALEXASC J235033.20-363303.2 | 357.63789 | -36.551 | G | 30.973623 | 0.00743366964161396 | 7e-06 | nan | 1e+20 | 1e+20 | 15.104 | 0.034 |
2dFGRS S274Z118 | 357.82371 | -29.26367 | Vis | 2.1414285 | 0.0005139427958056331 | 0.000213 | nan | 1e+20 | 1e+20 | 13.502 | 0.025 |
SNF 20080706-004 HOST | 358.78238 | 7.68939 | G | 169.34418 | 0.04064260423183441 | 5e-05 | nan | 1e+20 | 1e+20 | 15.167 | 0.038 |
SDSS-II SN 03674 | 359.79162 | 0.16298 | Vis | 149.9 | 0.035976000130176544 | -999.0 | 22.53733 | 1e+20 | 1e+20 | 1e+20 | 1e+20 |
'RA'=Right Ascension in degrees
'Dec'=Declination in degrees
'MODELMAG_R'=SDSS r-band magnitude
'MODELMAGERR_R'=SDSS r-band magnitude Error
'K_M_K20FE'=2MASS K-band magnitude
'K_MSIG_K20FE'=2MASS K-band magnitude Error
'W1MPRO'=WISE W1 magnitude (3.6 micron)
'W1SIGMPRO'=WISE W1 magnitude Error
Use the astropy.coordinates package and the SkyCoord function to store all of the galaxy catalog's locations.
The astropy coordinates package provides classes for representing a variety of celestial/spatial coordinates and their velocity components, as well as tools for converting between common coordinate systems in a uniform way. In addition, the astropy coordinates package facilitates fast manipulation and cross-matching. See here for examples: https://docs.astropy.org/en/stable/coordinates/
Create a coordinate object for the entire CLU catalog (hint: use SkyCoord).
clucoord=SkyCoord(ra=clu['RA']*u.deg,dec=clu['DEC']*u.deg)
Now let's read in the LIGO/VIRGO HEALpix map for GW170817.
LIGO/Virgo localization files for compact binary mergers include directional estimates of distance. The distance information is stored in three additional columns. To get the distance estimates, we need to ask for all four columns: PROB
, DISTMU
, DISTSIGMA
, and DISTNORM
.
url = 'https://dcc.ligo.org/public/0146/G1701985/001/preliminary-LALInference.fits.gz'
filename = astropy.utils.data.download_file(url)
prob, distmu, distsigma, distnorm = hp.read_map(filename, field=[0, 1, 2, 3])
npix = len(prob)
nside = hp.npix2nside(npix)
pixarea = hp.nside2pixarea(nside)
NSIDE = 1024 ORDERING = NESTED in fits file INDXSCHM = IMPLICIT Ordering converted to RING Ordering converted to RING Ordering converted to RING Ordering converted to RING
Find the coordinates of the highest probability pixel and put the coordinates into an astropy coordinate object called 'center'
ipix_max = np.argmax(prob)
ra_max, dec_max = hp.pix2ang(nside, ipix_max, lonlat=True)
center = SkyCoord(ra=ra_max*u.deg,dec=dec_max*u.deg)
print('Coordinates (RA,Dec) = %s' %(center))
Coordinates (RA,Dec) = <SkyCoord (ICRS): (ra, dec) in deg (197.35839844, -25.61308322)>
There are many visualization packages for plotting HEALPix maps. Luckily, LIGO has taken the time to provide its own user-friendly wrapper for plotting LIGO/VIRGO localizations.
Let's plot the sky localization using an 'astroglobe' projection centered on the highest highest probability pixel and overplot this location using the ligo.skymap package. (see here: https://lscsoft.docs.ligo.org/ligo.skymap/ligo/skymap/plot/allsky.html)
ax = plt.axes(
[0.05, 0.05, 0.9, 0.9],
projection='astro globe',
center=center)
ax.grid()
ax.imshow_hpx(filename, cmap='cylon')
ax.plot(
center.ra.deg, center.dec.deg,
transform=ax.get_transform('world'),
marker=ligo.skymap.plot.reticle(inner=0,outer=1),
markersize=10,
markeredgewidth=2)
[<matplotlib.lines.Line2D at 0x101e91048>]
/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1] /opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1] /opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1] /opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]
theta=0.5 * np.pi - clucoord.dec.to('rad').value
phi=clucoord.ra.to('rad').value
ipix=hp.ang2pix(nside,theta,phi)
#probability density per area on the sky for each galaxy
dp_dA=prob[ipix]/pixarea
clu['dP_dA']=dp_dA
#probability along radial distance
dp_dr=clu['DISTMPC']**2 * distnorm[ipix] * norm(distmu[ipix],distsigma[ipix]).pdf(clu['DISTMPC'])
#probability density per volume
dp_dV=prob[ipix] * distnorm[ipix] * norm(distmu[ipix],distsigma[ipix]).pdf(clu['DISTMPC'])/pixarea
clu['dP_dV']=dp_dV
clu.sort('dP_dA')
clu.reverse()
cum_sort=np.cumsum(clu['dP_dA'])
cumnorm_sort=cum_sort/np.max(cum_sort)
clu['P_A']=cumnorm_sort
#ID galaxies inside the 90% prob by volume
icutarea90=np.where(clu['P_A'] <= 0.9)
clucutarea90=clu[icutarea90]
#generate astropy coordinate object for this sample
clucutarea90coord=SkyCoord(ra=clucutarea90['RA']*u.deg,dec=clucutarea90['DEC']*u.deg)
print('# of galaxies in 90%% Area = %i' %(np.size(icutarea90)))
#sort the galaxies by P-value and print out top 20
clucutarea90['NEDname','dP_dA','P_A'][0:20].pprint(max_lines=22)
# of galaxies in 90% Area = 850 NEDname dP_dA P_A ---------------------------- ------------------ -------------------- 2MASX J13085466-2526013 328.4708057481784 0.007695115776384712 2MASX J13105449-2605323 313.15080651412734 0.015031328521702097 2MASX J13113020-2619386 300.98559528753657 0.022082545744755145 GALEXASC J131215.12-264630.4 282.127096965202 0.02869196319315675 GALEXASC J131108.75-255359.1 276.5933656229848 0.0351717414084558 ESO 508- G 003 260.04855554117916 0.04126392283630887 2MASX J13123929-2642121 257.9949138461368 0.04730799341129461 2MASX J13121969-2633515 257.30556126262 0.053335914459905154 2MASX J13124243-2721202 249.64494387777017 0.059184369520478404 2MASS J13122568-2628206 247.63513389002884 0.06498574057762256 2MASX J13132823-2713540 245.78170307672855 0.07074369114063976 2MASX J13111030-2655596 239.4362886006441 0.07635298709512517 2MASX J13060489-2351582 237.55153732271935 0.08191812880807826 UGCA 325 233.4818867376859 0.0873879304424818 ABELL 1664_06:[PSE2006] 0908 232.85121160651062 0.09284295719258284 ABELL 1664_12:[PSE2006] 2233 228.22178389917568 0.09818952991400635 2MASX J13130235-2636422 223.73585028023865 0.10343101025595866 2MASX J13130430-2742012 222.10294178363256 0.10863423630037909 AM 1302-232 NED01 221.5527942299063 0.1138245739874602 2MASX J13054598-2412025 218.38838045315967 0.11894077866195595
ax = plt.axes(
[0.05, 0.05, 0.9, 0.9],
projection='astro globe',
center=center)
#Zoomed-in inset window to better view the locations of the galaxies.
ax_inset = plt.axes(
[0.59, 0.3, 0.4, 0.4],
projection='astro zoom',
center=center,
radius=15*u.deg)
for key in ['ra', 'dec']:
ax_inset.coords[key].set_ticklabel_visible(False)
ax_inset.coords[key].set_ticks_visible(False)
ax.grid()
ax.mark_inset_axes(ax_inset)
ax.connect_inset_axes(ax_inset, 'upper left')
ax.connect_inset_axes(ax_inset, 'lower left')
ax_inset.scalebar((0.1, 0.1), 5 * u.deg).label()
ax_inset.compass(0.9, 0.1, 0.2)
ax.imshow_hpx('data/GW170817_prelim.fits.gz', cmap='cylon')
ax_inset.imshow_hpx('data/GW170817_prelim.fits.gz', cmap='cylon')
for coord in clucutarea90coord:
ax_inset.plot(
coord.ra.deg, coord.dec.deg,
transform=ax_inset.get_transform('world'),
marker=ligo.skymap.plot.reticle(inner=0,outer=1),
markersize=10,
markeredgewidth=1)
plt.show()
/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1] /opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1] /opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1] /opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]
Following the examples above, find galaxies in 90% VOLUME probability contour for GW170817, sort by Wise W1 luminosity, and overplot the top 20 sorted galaxies.
Information on WISE zeropoints and flux transformations http://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec4_4h.html
#load in CLU catalog
clu=Table.read('data/CLU_NEDonly.fits')
clucoord=SkyCoord(ra=clu['RA']*u.deg,dec=clu['DEC']*u.deg)
nclu=np.size(clu)
#make astropy coordinate object of CLU galaxies
clucoord=SkyCoord(ra=clu['RA']*u.deg,dec=clu['DEC']*u.deg)
#sky localization colmns to the galaxy catalog: probability density and p-value per volume and per area.
probdencol=Column(np.zeros(nclu,dtype='f4'),name='dP_dV')
probcol=Column(np.zeros(nclu,dtype='f4'),name='P')
probdenAcol=Column(np.zeros(nclu,dtype='f4'),name='dP_dA')
probAcol=Column(np.zeros(nclu,dtype='f4'),name='P_A')
clu.add_columns([probdencol,probcol,probdenAcol,probAcol])
#load in healpix map
prob,distmu,distsigma,distnorm=hp.read_map('data/GW170817_prelim.fits.gz',field=[0,1,2,3],dtype=('f8','f8','f8','f8'))
npix=len(prob)
nside=hp.npix2nside(npix)
pixarea=hp.nside2pixarea(nside)
#get coord of max prob density for plotting purposes
ipix_max = np.argmax(prob)
theta_max, phi_max = hp.pix2ang(nside, ipix_max)
ra_max = np.rad2deg(phi_max)
dec_max = np.rad2deg(0.5 * np.pi - theta_max)
center = SkyCoord(ra=ra_max*u.deg,dec=dec_max*u.deg)
print(center)
#calc hp index for each galaxy and populate CLU Table with the values
theta=0.5 * np.pi - clucoord.dec.to('rad').value
phi=clucoord.ra.to('rad').value
ipix=hp.ang2pix(nside,theta,phi)
#calc probability density per volume for each galaxy
dp_dV=prob[ipix] * distnorm[ipix] * norm(distmu[ipix],distsigma[ipix]).pdf(clu['DISTMPC'])/pixarea
clu['dP_dV']=dp_dV
#use normalized cumulative dist function to calculate Volume P-value for each galaxy
clu.sort('dP_dV')
clu.reverse()
cum_sort=np.cumsum(clu['dP_dV'])
cumnorm_sort=cum_sort/np.max(cum_sort)
clu['P']=cumnorm_sort
#ID galaxies inside the 90% prob by volume
icut90=np.where(clu['P'] <= 0.9)
clucut90=clu[icut90]
#generate an astropy coordinate object for this subset
clucut90coord=SkyCoord(ra=clucut90['RA']*u.deg,dec=clucut90['DEC']*u.deg)
print('# of galaxies in 90%% volume = %i' %(np.size(clucut90)))
NSIDE = 1024 ORDERING = NESTED in fits file INDXSCHM = IMPLICIT Ordering converted to RING Ordering converted to RING Ordering converted to RING Ordering converted to RING <SkyCoord (ICRS): (ra, dec) in deg (197.35839844, -25.61308322)> # of galaxies in 90% volume = 42
Q: Why are there so fewer galaxies in the volumne probability? A: Taking into account distance puts tighter constraints on the sample.
Part II - Sort by WISE W1 Luminosity
#Add a WISE W1 luminosity column to the CLU table
W1lumcol=Column(np.zeros(nclu,dtype='f8'),name='LumW1')
clu.add_column(W1lumcol)
#constansts needed
F0=309.540
clight=2.99792458e18 #Angstoms/sec
lamW1=33526. #Angtroms
fluxJyW1=F0*10**(-0.4*clucut90['W1MPRO']) #in Jy
fluxdenW1=fluxJyW1*1e-23 #erg/s/cm^2/Hz
freqW1=clight/lamW1
#distcm=4*np.pi*(np.float128(clucut90['DISTMPC'])*1.086e24)**2)
clucut90['LumW1']=fluxdenW1 * freqW1 * 4*np.pi*(np.float128(clucut90['DISTMPC'])*1.086e24)**2
#sort by WISE 1 Luminosity (proportional to galaxy stellar mass)
clucut90.sort('LumW1')
clucut90.reverse()
#then print list of prioritized galaxies
clucut90['NEDname','LumW1','dP_dV','P'][0:20].pprint(max_lines=22)
#Q: Is NGC4993 in your list?
#A: It should be #2
NEDname LumW1 ... P ----------------------- ------------------------- ... ------------------- NGC 4970 1.0884088954320000597e+42 ... 0.6995353723555456 NGC 4993 7.5032646027206775e+41 ... 0.7965022626584416 NGC 4968 7.367336809737711801e+41 ... 0.3905416510700604 IC 4197 7.159408041749106282e+41 ... 0.636392013632011 IC 4180 6.191420321536100089e+41 ... 0.20082875617890292 ESO 508- G 033 3.728938562939855379e+41 ... 0.7710150928931331 IC 0874 2.957639646390979745e+41 ... 0.8792942624546624 ESO 508- G 010 1.3404951718770816046e+41 ... 0.6805763434781986 2MASX J13161781-2926415 1.0298766866197851405e+41 ... 0.8733888921183204 ESO 575- G 053 8.7703046690865904305e+40 ... 0.45294243100812964 IC 0879 7.7850962931557685694e+40 ... 0.890283449648996 UGCA 331 6.182574842412703268e+40 ... 0.808623779107441 ESO 575- G 044 NED01 5.705652214096152376e+40 ... 0.895230316115733 ESO 575- G 055 4.918953488412471152e+40 ... 0.8672853532554123 ESO 508- G 003 4.8158691124006157543e+40 ... 0.0654119434763518 ESO 508- G 019 3.995895285635972178e+40 ... 0.5098985022830953 2MASX J13073768-2356181 3.4884431466584360636e+40 ... 0.8532464568147246 2MFGC 10461 3.3076261973941396845e+40 ... 0.5885153130290641 2MASX J13061939-2258491 2.7938869474082211623e+40 ... 0.6592093486395489 UGCA 327 2.653563605837436326e+40 ... 0.7574665230057527
#plot up the sky localization and overplot the galaxies
ax = plt.axes(
[0.05, 0.05, 0.9, 0.9],
projection='astro globe',
center=center)
ax_inset = plt.axes(
[0.59, 0.3, 0.4, 0.4],
projection='astro zoom',
center=center,
radius=10*u.deg)
for key in ['ra', 'dec']:
ax_inset.coords[key].set_ticklabel_visible(False)
ax_inset.coords[key].set_ticks_visible(False)
ax.grid()
ax.mark_inset_axes(ax_inset)
ax.connect_inset_axes(ax_inset, 'upper left')
ax.connect_inset_axes(ax_inset, 'lower left')
ax_inset.scalebar((0.1, 0.1), 5 * u.deg).label()
ax_inset.compass(0.9, 0.1, 0.2)
ax.imshow_hpx('data/GW170817_prelim.fits.gz', cmap='cylon')
ax_inset.imshow_hpx('data/GW170817_prelim.fits.gz', cmap='cylon')
for coord in clucut90coord:
ax_inset.plot(
coord.ra.deg, coord.dec.deg,
transform=ax_inset.get_transform('world'),
marker=ligo.skymap.plot.reticle(inner=0,outer=1),
markersize=10,
markeredgewidth=1)
#where is NGC4993? hint: use ax_inset.text()
c4993=SkyCoord.from_name('NGC 4993')
ax_inset.text(c4993.ra.deg+10.5, c4993.dec.deg,'NGC 4993',transform=ax_inset.get_transform('world'),fontdict={'size':10,'color':'black','weight':'normal'})
plt.show()
/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1] /opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1] /opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1] /opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]