Gravitational Wave Localizations and Galaxy Crossmatch Module

Lecturer: Leo Singer
Jupyter Notebook Author: Leo Singer, Dave Cook, & Cameron Hummels

This is a Jupyter notebook lesson taken from the GROWTH Winter School 2018. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.html

Objective

Learn how to use LIGO/Virgo localizations and match with galaxies

Key steps

  • Manipulate HEALPix localization files
  • Cross-match a LIGO localization with a galaxy catalog

Required dependencies

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).

Python modules

  • python 3
  • astropy
  • numpy
  • scipy
  • matplotlib
  • healpy
  • ligo.skymap

External packages

None

Imports

First, some imports: Numpy, Matplotlib, Healpy, and parts of Astropy.

In [1]:
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:

In [2]:
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:

In [3]:
%matplotlib inline

HEALPix Basics

This section on using HEALPix localization files is adapted from the LIGO/Virgo Public Alerts User Guide.

Download and read localization

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.

In [4]:
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.

In [5]:
prob = hp.read_map(filename)
NSIDE = 2048
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING

Manipulating HEALPix Coordinates

To get a quick look at a HEALPix data set, you can use the hp.mollview function:

In [6]:
hp.mollview(prob)

What actually is stored in prob?

In [7]:
prob
Out[7]:
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

$$ \mathsf{npix} = 12 \cdot \mathsf{nside}^2. $$

The functions hp.npix2nside and hp.nside2npix convert between length and resolution.

In [8]:
npix = len(prob)
npix
Out[8]:
50331648
In [9]:
nside = hp.npix2nside(npix)
nside
Out[9]:
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.

In [10]:
ipix = 123
ra, dec = hp.pix2ang(nside, ipix, lonlat=True)
ra, dec
Out[10]:
(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.

In [11]:
ra = 194.95
dec = 27.98
hp.ang2pix(nside, ra, dec, lonlat=True)
Out[11]:
13361492

What is the most probable sky location? Just find the pixel with the maximum valuem, and then find its right ascension and declination.

In [12]:
ipix_max = np.argmax(prob)
ipix_max
Out[12]:
32883013
In [13]:
hp.pix2ang(nside, ipix_max, lonlat=True)
Out[13]:
(194.30419921875, -17.856895095545468)

Probability distributions with scipy.stats

In [14]:
#Imagine I have a t statistic with 3 degrees of freedom.
#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:

# Make a t distribution object for t with 3 degrees of freedom
t_dist = scipy.stats.t(3)
# Plot the PDF
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 with 20 degrees of freedom. 
#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 of (here).
In [15]:
# Show relationship of PDF and CDF for three example t values.
#give the students most of the code for these plots
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 
#with 20 degrees of freedom. The probability that x<=1.5 is:

# Area of PDF at and to the left of 1.5
t_dist.cdf(1.5)
#>0.9253...

#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:

# Area of PDF to the right of 1.5
1 - t_dist.cdf(1.5)
#>0.0746...

#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.


#have the students fill out the example for 90% or something
Out[15]:
0.11529193262241144

Working with LIGO/Virgo 3D localizations

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.

In [16]:
url = 'https://dcc.ligo.org/public/0146/G1701985/001/bayestar.fits.gz'
filename = astropy.utils.data.download_file(url)
In [17]:
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 = 2048
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING
Ordering converted to RING
Ordering converted to RING
Ordering converted to RING
In [18]:
#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)

#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])
clu
Out[18]:
Table length=225709
CLUIDID_OTHERNAMERADECDMDM_KINBTCARATIO_B2APATYPEM21RA_NEDDEC_NEDZZERRA_NEDB2A_NEDPA_NEDTYPE_NEDSOURCEB_TB_TERRB_R25B_R25ERRFUV_1FUVERR_1NUV_1NUVERR_1RA_SDSSDEC_SDSSPETROMAG_UPETROMAG_GPETROMAG_RPETROMAG_IPETROMAG_ZPETROMAGERR_UPETROMAGERR_GPETROMAGERR_RPETROMAGERR_IPETROMAGERR_ZMODELMAG_UMODELMAG_GMODELMAG_RMODELMAG_IMODELMAG_ZMODELMAGERR_UMODELMAGERR_GMODELMAGERR_RMODELMAGERR_IMODELMAGERR_ZDESIGNATION_WISERA_WISE_1DEC_WISE_1W1MPROW1SIGMPROW1SNRW2MPROW2SIGMPROW2SNRW3MPROW3SIGMPROW3SNRW4MPROW4SIGMPROW4SNRDESIGNATION_2MASSRA_2MASSDEC_2MASSR_K20FEJ_M_K20FEJ_MSIG_K20FEJ_FLG_K20FEH_M_K20FEH_MSIG_K20FEH_FLG_K20FEK_M_K20FEK_MSIG_K20FEK_FLG_K20FECOORD_SOURCEBTC_SOURCESIZE_SOURCEDM_SOURCEDM_KIN_OTHERDM_RANGEZ_SOURCEDM_FLAGRA_GALEXDEC_GALEXNUV_2NUVERR_2FUV_2FUVERR_2FLAGS_FUVRA_WISE_2DEC_WISE_2WISE1WISE1ERRWISE4WISE4ERRFUVNUVFUVERRNUVERRMAGBLUM_BSFR_FUVMSTARDISTMPCHAGALHACANDdP_dVPdP_dAP_A
int64bytes40bytes30float64float64int64float64float64float64float64float64bytes40float64float64float64float64float64float64float64float64bytes3bytes12float64float64float64float64float64float64float64float64float64float64float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32bytes19float64float64float64float64float64float64float64float64float64float64float64float64float64float64bytes16float64float64float64float64float64int64float64float64int64float64float64int64bytes3bytes17bytes17bytes17bytes47float64bytes3int64float64float64float32float32float32float32int16float64float64float32float32float32float32float32float32float32float32float32float64float64float64float32int16int16float32float32float32float32
231192N/A2dFGRS S805Z4170.00162-56.1410699999933.18027481140323nan1e+201e+201e+20N/A1e+200.00162-56.141060.0103816464543342590.0002971e+201e+201e+20GNED_20161120nan1e+20nan1e+2020.34190.15838620.03650.0858811e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J000000.26-560827.50.0011127-56.140987215.9830.04623.715.9180.1248.712.5581e+200.39.0021e+20-0.4N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz33.09060.08967481140322775NED00.0026155762934943604-56.1407478334413920.03650.0858809820.341930.1583855500.0011127-56.140987215.9830.0469.0020.520.341920.03650.158386nannannan0.0164270836265184479255697.8624861443.25686000.00.00.00.0
231193472509PGC5973600.0068-39.5130999999933.3151041868750117.841e+200.81283175.0N/A1e+200.0068-39.513090.0110466880723834040.0001191e+201e+201e+20GNED_20161120nan1e+20nan1e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20N/A1e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedclu.v2clu.v2nedz33.2254/33.18000.13510418687501158NED00.007313433062165586-39.5122984629911520.249180.1194300720.617890.2075267300.00.0nan0.0nan0.020.6178920.249180.20752673nan17.845.307357561598291e+41nannan46.027866000.00.00.00.0
231195497360PGC60.0090515.8818899999934.6656628883578415.081e+200.70794687.6-0.300000009.290.0090515.881880.0205751881003379823.3e-0519.00.7105.0GNED_20161120nan1e+20nan1e+2017.89480.052306517.50150.02428580.008745004302170415.88171271506376816.57777815.44914814.9218714.65134314.4028010.022480890.0061467030.0045977460.00564722620.0108778516.57643315.46315714.91755714.64105614.3857080.0130515050.00294021730.00270334750.00270710440.0052141347J000002.18+155253.80.009100415.881626312.1780.02346.612.0920.02348.18.110.02152.56.0190.05121.500000214+15525390.00895315.881669.513.3270.039012.5840.056012.310.0670nedclu.v2+a.40clu.v2+a.40nedz34.6900/34.57600.11399999999999721NED00.00910112526844493615.8815058617398317.5015050.02428581617.8948290.05230654400.009100415.881626312.1780.0236.0190.05117.894817.50150.0523065nan15.082.339386031800773e+430.879944101456220610355846084.14297185.72995000.00.00.00.0
231196474856PGC7505060.00996-27.7419499999935.29073622314846418.311e+200.47863164.0N/A1e+200.00996-27.741940.027438381686806680.0002971e+201e+201e+20GNED_20161120nan1e+20nan1e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J000002.36-274430.90.0098518-27.741934315.2680.0426.815.1560.09112.012.1731e+201.38.9371e+20-0.4N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedclu.v2clu.v2nedz35.2500/35.20110.08963622314846731NED00.009707923128452632-27.7420348416773721.1672170.278063821.4610790.4363927500.0098518-27.741934315.2680.048.9370.521.46107921.1672170.43639275nan18.312.1238871411927357e+420.095674911511479081069583801.8125993114.32659000.00.00.00.0
231200528573PGC6753910.01362-33.5096799999934.96730712257798518.571e+200.691831144.0N/A1e+200.01362-33.509670.0236413720995187760.0002971e+201e+201e+20GNED_20161120nan1e+20nan1e+2020.66040.15174320.09160.07667591e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J000003.20-333033.70.0133729-33.509369716.6510.08313.117.3030.5152.112.2851e+200.58.3871e+201.5N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedclu.v2clu.v2nedz34.9100/34.87760.08970712257798397NED00.013499105611231244-33.5092903415969320.09160.07667585520.660370.1517431300.0133729-33.509369716.6510.0838.3870.520.660420.09160.151743nan18.571.2409655175467712e+420.12211077330362756222148532.9807124498.505714000.00.00.00.0
231206534792PGC11240390.02304-1.215599999935.574317.591e+200.7413186.9N/A1e+200.02304-1.21550.0312659218907356261e+201e+201e+201e+20UvSNED_20161120nan1e+20nan1e+201e+201e+201e+201e+200.023088624282308956-1.216374038152013218.62214317.7072617.17064717.02219216.6295490.095904490.0325833450.0160633230.0288439540.0696493518.6433717.64032217.15783516.92105516.7080970.043929440.009315410.009245760.0110898110.032711327J000005.57-011257.60.0232444-1.216023615.2920.0427.115.1080.08912.212.3980.4672.38.5541e+200.9N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedclu.v2clu.v2ned_2012112035.66000.08569999999999567N/A00.023040469206642218-1.21549750910660318.9216480.03607401619.0577370.0593869800.0232444-1.216023615.2920.048.5540.519.05773718.9216480.05938698nan17.595.3524990050153495e+420.3095167168516951358440700.9837348130.27467000.00.00.00.0
231207N/AUSGC U0010.0233326.1902899999935.135299907836774nan1e+201e+201e+20N/A1e+200.0233326.190280.0255429614335298541e+201e+201e+201e+20GTrNED_20161120nan1e+20nan1e+201e+201e+201e+201e+200.0169259660992793226.19676295402018724.76850923.12087822.39664821.95252421.497282.75815250.5922760.475191150.69665131.086400424.73939324.27878622.11059421.2519920.5618971.23949440.711362960.170685470.113038240.24993625J000005.21+261122.00.021716326.189456517.7190.2135.116.6621e+201.112.0880.4892.28.4571e+20-1.0N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz35.04560.08969990783677417NED00.00.00.00.00.00.000.021716326.189456517.7190.2138.4570.5nannan0.0nan25.256023.0658317970055924e+39nan96970633.79079935106.429000.00.00.00.0
231208231209DUKST 349-0640.02492-35.6582299999935.42525964272521nan1e+201e+201e+20N/A1e+200.02492-35.658220.029191952198743820.0002131e+201e+201e+20GNED_20161120nan1e+20nan1e+2018.7590.05915218.44430.0156661e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J000006.03-353930.00.0251585-35.658344414.6990.04623.614.5740.0618.111.170.1338.28.7111e+200.6N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz35.3356/35.29590.1293596427252055NED00.024983397999790213-35.6581785808850318.4211650.02931013718.7589950.05915200700.0251585-35.658344414.6990.0468.7110.518.75918.44430.059152nannannan0.295339127255143832044675021.7072265121.63314000.00.00.00.0
231210N/AUGC 128900.029318.2791899999936.097222767584576nan1e+201e+201e+20E1e+200.029318.279180.0397791750729084-999.060.00.91e+20GNED_20161120nan1e+20nan1e+201e+201e+2020.78610.265441e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J000007.03+081645.10.029338.2792211.5210.02348.111.6030.02248.610.4840.1268.68.6381e+200.300000701+08164480.0292258.27911823.611.8620.032611.1010.043610.7790.0596nedN/AN/Anedz36.00750.08972276758457554NED00.029219782815484018.27935205893599920.7861120.265440321.591870.5833424320.029338.2792211.5210.0238.6380.521.5918720.78610.58334243nannannan0.252870931433639370894043111.6906165.74657000.00.00.00.0
231213497361PGC120.036-6.37499999934.8547322507632314.030.05044440.3170.01.20000001e+200.036-6.3740.0224469695240259172.7e-0552.60.3170.0GNED_20161120nan1e+20nan1e+2017.71130.0678117.29680.03684060.03612959593030496-6.374013573347971516.05697414.78959714.1465413.78324113.47853850.037782530.013460680.0072698240.0080855710.01948488716.27774614.896173514.14183213.74109313.3972720.0175067860.00271096080.00248015160.0024356640.004228183J000008.65-062226.40.036046-6.374018411.9460.02348.211.9010.02642.58.5290.02740.56.660.09611.300000865-06222630.036053-6.37399926.312.2740.027011.5240.036011.240.0590nedclu.v2ned_20121120nedz34.8400/34.76510.08963225076323056NED00.03616868643439375-6.374222114506590516.8094670.03452976417.8117030.09128930400.036046-6.374018411.9460.0236.660.09617.711317.29680.06781nan14.037.323686388808671e+430.7293206317194915262072350.89321393.52904000.00.00.00.0
............................................................................................................................................................................................................................................................................................................................................................
654837N/AGALEXASC J234305.34-552311.9355.77191-55.3866899999934.663817920809926nan1e+201e+201e+20N/A1e+20355.77191-55.386680.020557714626193047-999.01e+201e+201e+20VisNED_20161120nan1e+20nan1e+2021.57270.17599220.91720.04285971e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J234305.38-552312.4355.7724516-55.386794216.1630.04822.415.1030.06217.612.3330.4192.68.2661e+201.0N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz34.63400.029817920809925624NED0355.7721764234055-55.38679042625647520.8370130.1489137121.9286040.413019270355.7724516-55.386794216.1630.0488.2660.521.572720.91720.175992nannannan0.09256044687004351263298097.9673837785.65714000.00.00.00.0
654838N/AGALEXASC J234344.33-545403.7355.93484-54.9011999999933.15866794249002nan1e+201e+201e+20N/A1e+20355.93484-54.901190.010278857313096523-999.01e+201e+201e+20VisNED_20161120nan1e+20nan1e+2020.21370.092334119.21480.01651941e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J234344.41-545403.5355.9350627-54.900992115.0470.03233.413.8560.03432.210.9560.1278.58.5580.3812.9N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz33.12900.029667942490021915NED0355.9347096300996-54.9010409786521119.2580.0348957520.2136960.092334130355.9350627-54.900992115.0470.0328.5580.38120.213719.21480.0923341nannannan0.022261448517084066183987110.1889945342.82857000.00.00.00.0
654842N/ASDSS J234403.01+003122.6356.012550.5229599999934.95860381469648nan1e+201e+201e+20N/A1e+20356.012550.522950.023546807467937471e+203.350.97122.0GNED_20161120nan1e+20nan1e+201e+201e+201e+201e+20356.01395395080610.523619933437692522.00231621.42951621.68138322.85100622.8269710.312078770.23191170.357829871.20274581.956809921.8476421.22722821.5164722.61219422.8269080.15591190.0512051250.085009250.358398380.7732769N/A1e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz34.92900.02960381469647899NED0356.01120094140980.52134503880309517.6400990.02826610818.1254670.04627340330.00.0nan0.0nan0.018.12546717.6400990.046273403nan21.0132871.2970715057836036e+41nannan98.111694000.00.00.00.0
654845N/AGALEXASC J234414.08-364635.4356.05836-36.7765599999932.115015347041364nan1e+201e+201e+20N/A1e+20356.05836-36.776550.0063564451411366461.3e-051e+201e+201e+20GNED_20161120nan1e+20nan1e+2018.26180.054834118.18980.0321971e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J234414.03-364634.7356.0584841-36.776308715.0550.03630.514.7790.06117.712.2211e+201.28.2671e+201.6N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz32.08600.029015347041365658NED0356.0586949332161-36.7765091352403818.1897640.0321970418.2617550.0548340680356.0584841-36.776308715.0550.0368.2670.518.261818.18980.0548341nannannan0.0217243493058917969843546.5578939626.48519000.00.00.00.0
654864N/A2MASX J23454008+0144274356.417041.7409299999935.3877005203879116.371e+201e+201e+20N/A1e+20356.417041.740920.0286913737654685971e+2015.40.7165.0GNED_20161120nan1e+20nan1e+2018.6960.080230718.41480.0428454356.417048928266241.740976337558005817.28971716.20464315.69939915.39794915.1932320.0405788420.00642992140.00591019630.0070233010.02013248217.35351416.16076515.65992515.34259915.1289310.0227814730.00413906430.0038023340.00413068660.010281904J234540.10+014427.3356.41711631.740928814.0730.0336.213.890.04623.410.5860.10310.58.0191e+201.923454008+0144274356.4170231.7409597.714.8360.086014.1980.114014.0840.1680nedNED_20150113N/Anedz35.35800.029700520387912377NED0356.417518137019331.741270849657479618.414790.0428454318.6960030.080230680356.41711631.740928814.0730.038.0190.518.69618.41480.0802307nan16.371.3864898197522259e+430.39466617044460213515610064.078043119.54739000.00.00.00.0
654880N/AAGC 331333356.7534927.263199999935.53731436266381nan1e+201e+201e+20N/A1e+20356.7534927.26310.0307378936558961870.000321e+201e+201e+20GNED_20161120nan1e+20nan1e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20N/A1e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz35.50800.02931436266380416NED00.00.00.00.00.00.00356.7533827.260978116.6770.0858.8020.5nannan0.0nannannannan366645225.4569402128.07455000.00.00.00.0
654908N/AGALEXASC J235033.20-363303.2357.63789-36.55199999932.45496002982655nan1e+201e+201e+20N/A1e+20357.63789-36.5510.007433669641613967e-061e+201e+201e+20GNED_20161120nan1e+20nan1e+2018.3810.065753918.15170.03664781e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J235033.13-363301.7357.6380696-36.550484815.1040.03432.115.1380.07514.412.2311e+200.68.3451e+201.8N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz32.42600.02896002982654977NED0357.6383747584509-36.5508962292272818.1516880.03664782318.3810060.0657538850357.6380696-36.550484815.1040.0348.3450.518.38118.15170.0657539nannannan0.02701093033610680591307044.6905354730.973623000.00.00.00.0
654916N/A2dFGRS S274Z118357.82371-29.2636799999926.65351796417011319.131e+201e+201e+20N/A1e+20357.82371-29.263670.00051394279580563310.0002131e+201e+201e+20VisNED_20161120nan1e+20nan1e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J235117.70-291550.4357.8237579-29.264017913.5020.02542.813.3810.0335.712.0551e+201.98.9031e+200.0N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedNED_20150113N/Anedz26.62400.029517964170114652NED00.00.00.00.00.00.00357.8237579-29.264017913.5020.0258.9030.5nannan0.0nan19.133.501418754522525e+38nan1908655.0664323762.1414285000.00.00.00.0
654960N/ASNF 20080706-004 HOST358.782387.6893999999936.14385126737829nan1e+201e+201e+20N/A1e+20358.782387.689390.040642604231834415e-051e+201e+201e+20GNED_20161120nan1e+20nan1e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20nannannan1e+201e+201e+201e+201e+201e+20J235507.69+074122.7358.78208317.689643115.1670.03828.615.0690.0912.012.2520.3992.78.9051e+20-1.3N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz36.11400.029851267378290913NED00.00.00.00.00.00.00358.78208317.689643115.1670.0388.9050.5nannan0.0nannannannan2575495771.518145169.34418000.00.00.00.0
654995N/ASDSS-II SN 03674359.791620.1629899999935.8790081642414nan1e+201e+201e+20N/A1e+20359.791620.162980.035976000130176544-999.01e+201e+201e+20VisNED_20161120nan1e+20nan1e+2023.08130.20435623.70410.331849359.79156614372150.1630237553981678222.2474422.66706522.49953322.4936420.7426760.4827710.32397590.520717740.643079760.539378723.47423722.78463622.5373323.14647722.01590.523076240.142466860.162672880.377172080.53411627N/A1e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+201e+20N/A1e+201e+201e+201e+201e+2011e+201e+2011e+201e+201nedN/AN/Anedz35.85000.029008164241396628NED00.00.00.00.00.00.00359.78896890.163195316.0930.0628.7350.523.081323.7041nannan22.846945.593250283121446e+400.17582424577163125860049955.2722148149.9000.00.00.00.0
In [19]:
#get coord of max prob density for plotting purposes
#ipix_max=np.where(prob == np.max(prob))
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)
center
Out[19]:
<SkyCoord (ICRS): (ra, dec) in deg
    (194.30419922, -17.8568951)>
In [20]:
#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.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)
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy_healpix/core.py:519: RuntimeWarning: invalid value encountered in bilinear_interpolation_weights
  result = _core.bilinear_interpolation_weights(lon, lat, nside)
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy_healpix/core.py:484: RuntimeWarning: invalid value encountered in ring_to_nested
  return _core.ring_to_nested(ring_index, nside)
Out[20]:
[<matplotlib.lines.Line2D at 0x7ffb973f0c18>]
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater
  discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater
  discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater
  discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater
  discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]
In [21]:
#calc 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'])

#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
In [22]:
#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]
clucut90coord=SkyCoord(ra=clucut90['RA']*u.deg,dec=clucut90['DEC']*u.deg)

print('# of galaxies in 90%% volume = %i' %(np.size(icut90)))

#sort by WISE 1 Luminosity (proportional to galaxy stellar mass)
clucut90.sort('P')
clucut90.reverse()

#then print list of prioritized galaxies
clucut90['NAME','dP_dV','P'].pprint(max_lines=20)
# of galaxies in 90% volume = 3047
          NAME                   dP_dV                    P           
------------------------ ---------------------- ----------------------
                PGC47394  6.317992998346537e-64      0.899883254360382
                PGC82565  6.317992998346537e-64     0.8997641503042281
                PGC49468  6.317992998346537e-64     0.8996450462480743
                PGC37565  6.317992998346537e-64     0.8995259421919205
                PGC70582  6.317992998346537e-64     0.8994068381357666
SDSS_J124305.87+242703.8  6.334407010119384e-64     0.8992877340796128
          MCG -02-05-043  6.334407010119384e-64      0.899168320593648
                PGC40483   6.36402347040082e-64     0.8990489071076831
                     ...                    ...                    ...
                UGCA 339 2.2638103395560794e-63  0.0034142113450347576
         HIPASS J1008-33 2.2638525648239792e-63   0.002987447665691437
          [TSK2008] 0687  2.263866444423913e-63   0.002560676026224261
SDSS J105200.35+274532.7 2.2638875317927865e-63  0.0021339017702352706
               AGC748687 2.2639016072753102e-63  0.0017071235389472146
                PGC26218 2.2639016072753102e-63   0.001280342654210411
                PGC43746 2.2639016072753102e-63  0.0008535617694736073
                PGC32259 2.2639016072753102e-63 0.00042678088473680365
Length = 3047 rows

Exercise for students

Following the examples above, find galaxies in 90% prob contour for GW170817, then sort by Wise1 luminosity.

Information on WISE zeropoints and flux transformations http://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec4_4h.html

In [23]:
#load in CLU catalog
clu=Table.read('data/CLU_NEDonly.fits')
#clu=Table.read('data/CLU_20170106_galexwise_DaveUpdate.fits')
clucoord=SkyCoord(ra=clu['RA']*u.deg,dec=clu['DEC']*u.deg)
nclu=np.size(clu)

#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.where(prob == np.max(prob))
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
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]
clucut90coord=SkyCoord(ra=clucut90['RA']*u.deg,dec=clucut90['DEC']*u.deg)

print('# of galaxies in 90%% volume = %i' %(np.size(icut90)))
#sort the galaxies by WISE band 1 luminosity,then print list of prioritized galaxies
W1lumcol=Column(np.zeros(nclu,dtype='f8'),name='LumW1')
clu.add_column(W1lumcol)

#constants 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=n4*np.pi*(cp.float128(lucut90['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['NAME','LumW1','dP_dV','P'].pprint(max_lines=20)
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
            NAME                       LumW1           ...          P         
---------------------------- ------------------------- ... -------------------
                    PGC45466 1.0884088954320000597e+42 ...  0.6995353723555456
                    PGC45657    7.5032646027206775e+41 ...  0.7965022626584416
                    PGC45426  7.367336809737711801e+41 ...  0.3905416510700604
                    PGC45514  7.159408041749106282e+41 ...   0.636392013632011
                    PGC45408  6.191420321536100089e+41 ... 0.20082875617890292
                    PGC46199  3.728938562939855379e+41 ...  0.7710150928931331
                    PGC46410  2.957639646390979745e+41 ...  0.8792942624546624
                    PGC45475 1.3404951718770816046e+41 ...  0.6805763434781986
                         ...                       ... ...                 ...
ABELL 1664_13:[PSE2006] 1313 2.1890173247059815035e+39 ...  0.8851443911817604
                    PGC45679 1.9114414792285373049e+39 ... 0.42337005469909356
ABELL 1664_03:[PSE2006] 2951  6.403713666703690913e+38 ...  0.8999323517709407
                  PGC4075686                       0.0 ...  0.5624329696124779
            HIPASS J1321-31b                       0.0 ...   0.838117879353231
                    PGC45404                       0.0 ...  0.7146971370109703
                    UGCA 325                       0.0 ... 0.15857534419648253
              [TSK2008] 0052                       0.0 ...  0.7293957356412174
Length = 42 rows

Exercise for students

Now, plot up the sky localization and overplot the galaxies on it.

In [24]:
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)
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy_healpix/core.py:519: RuntimeWarning: invalid value encountered in bilinear_interpolation_weights
  result = _core.bilinear_interpolation_weights(lon, lat, nside)
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy_healpix/core.py:484: RuntimeWarning: invalid value encountered in ring_to_nested
  return _core.ring_to_nested(ring_index, nside)
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater
  discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater
  discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater
  discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]
/home/chummels/miniconda3/lib/python3.7/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73: RuntimeWarning: invalid value encountered in greater
  discontinuous = step[1:] > DISCONT_FACTOR * step[:-1]