# 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 3
• astropy
• numpy
• scipy
• matplotlib
• healpy
• ligo.skymap

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.

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'

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

Out[18]:
Table length=225709
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

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

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

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

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]