GW Example

v2 – Refactor with localization

v3 – Refactor with PATH

[1]:

%matplotlib notebook

[16]:

# imports
from importlib import reload
import os
from pkg_resources import resource_filename
import numpy as np

import healpy as hp
import pandas

from astropy.io import fits
from astropy.table import Table
from astropy.coordinates import SkyCoord

from astropath import bayesian
from astropath import localization
from astropath import path

LIGO example – GW170817

[3]:

lfile = os.path.join(resource_filename('astropath', 'data'), 'gw_examples',
                     'GW170817_skymap.fits.gz')

[4]:

gw170817 = hp.read_map(lfile)

/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/fitsfunc.py:369: UserWarning: If you are not specifying the input dtype and using the default np.float64 dtype of read_map(), please consider that it will change in a future version to None as to keep the same dtype of the input file: please explicitly set the dtype if it is important to you.
  "If you are not specifying the input dtype and using the default "
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/fitsfunc.py:391: UserWarning: NSIDE = 1024
  warnings.warn("NSIDE = {0:d}".format(nside))
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/fitsfunc.py:400: UserWarning: ORDERING = NESTED in fits file
  warnings.warn("ORDERING = {0:s} in fits file".format(ordering))
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/fitsfunc.py:428: UserWarning: INDXSCHM = IMPLICIT
  warnings.warn("INDXSCHM = {0:s}".format(schm))
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/fitsfunc.py:486: UserWarning: Ordering converted to RING
  warnings.warn("Ordering converted to RING")

[5]:

header = fits.open(lfile)[1].header
header

[5]:

XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                   32 / length of dimension 1
NAXIS2  =             12582912 / length of dimension 2
PCOUNT  =                    0 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    4 / number of table fields
TTYPE1  = 'PROB    '
TFORM1  = 'D       '
TUNIT1  = 'pix-1   '
TTYPE2  = 'DISTMU  '
TFORM2  = 'D       '
TUNIT2  = 'Mpc     '
TTYPE3  = 'DISTSIGMA'
TFORM3  = 'D       '
TUNIT3  = 'Mpc     '
TTYPE4  = 'DISTNORM'
TFORM4  = 'D       '
TUNIT4  = 'Mpc-2   '
PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation
ORDERING= 'NESTED  '           / Pixel ordering scheme: RING, NESTED, or NUNIQ
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)
NSIDE   =                 1024 / Resolution parameter of HEALPIX
INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT
OBJECT  = 'GW170817'           / Unique identifier for this event
DATE-OBS= '2017-08-17T12:41:04.429464' / UTC date of the observation
MJD-OBS =    57982.52852348908 / modified Julian date of the observation
DATE    = '2019-05-08T21:51:56.613850' / UTC date of file creation
CREATOR = 'ligo-skymap-from-samples' / Program that created this file
ORIGIN  = 'LIGO/Virgo'         / Organization responsible for this FITS file
DISTMEAN=    38.03408225450563 / Posterior mean distance (Mpc)
DISTSTD =    7.499686641911211 / Posterior standard deviation of distance (Mpc)
VCSVERS = 'ligo.skymap 0.0.15' / Software version
VCSREV  = '65fc6500a1e117fec2e27550a8e9b10c9792ffca' / Software revision (Git)
DATE-BLD= '2018-09-04T14:19:20' / Software build date
HISTORY
HISTORY Generated by running the following script:
HISTORY ligo-skymap-from-samples --seed 170817 --outdir . --objid GW170817 --fit
HISTORY soutname MCMC_TF2_LowSpin_AllSky.fits.gz MCMC_TF2_LowSpin_AllSky_post.da
HISTORY t

[6]:

hp.mollview(
    gw170817,
    coord=["C"],
    title="GW170817",
    #unit="mK",
    #norm="hist",
    #min=-1,
    #max=1,
)
hp.graticule()

/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/projaxes.py:920: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  newcm.set_over(newcm(1.0))
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/projaxes.py:921: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  newcm.set_under(bgcolor)
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/projaxes.py:922: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  newcm.set_bad(badcolor)
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/projaxes.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
  **kwds
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/projaxes.py:543: UserWarning: 0.0 180.0 -180.0 180.0
  pmin / dtor, pmax / dtor, mmin / dtor, mmax / dtor
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/projaxes.py:658: UserWarning: The interval between parallels is 30 deg -0.00'.
  vdeg, varcmin
/home/xavier/anaconda3/lib/python3.7/site-packages/healpy/projaxes.py:666: UserWarning: The interval between meridians is 30 deg -0.00'.
  vdeg, varcmin

Galaxies

https://vizier.u-strasbg.fr/viz-bin/VizieR?-source=VII/275

[7]:

galfile = os.path.join(resource_filename('astropath', 'data'), 'gw_examples',
                     'GW170817_galaxies.csv')

[8]:

cut_galaxies = pandas.read_csv(galfile, index_col=0)

[9]:

cut_galaxies.head()

[9]:
RAJ2000 DEJ2000 Dist Bmag BmagHyp ImagHyp modzHyp mod0Hyp logd25Hyp logr25Hyp ... BmagGWGC maj min PAGWGC DistGWGC Kmag2MPZ Bmag2MPZ zsp2MPZ zph2MPZ Flag
667013 357.169220 -2.362345 92.028 15.33 15.81 14.02 34.92 NaN 0.85 0.19 ... 15.33 0.708 0.457 136.2 92.028 12.761 15.622 0.02219 0.031973 0
667018 150.132019 -31.244976 31.583 12.60 12.95 10.66 32.62 NaN 1.31 0.15 ... 12.60 2.042 1.446 179.9 31.583 8.792 12.490 0.00813 0.008867 0
667030 76.325119 -9.147248 44.347 12.53 13.03 10.86 33.33 NaN 1.43 0.35 ... 12.53 1.816 0.890 112.0 44.347 8.818 12.735 0.01090 0.010972 0
667044 338.249908 -20.094194 80.264 14.13 14.43 13.35 34.57 NaN 1.05 0.03 ... 14.13 0.830 0.772 NaN 80.264 12.792 14.575 0.01905 0.011343 0
667076 138.796539 29.252758 88.167 16.19 17.13 NaN 34.83 NaN 0.72 0.46 ... 16.19 0.525 0.182 2.4 88.167 13.127 16.320 0.02080 0.062994 0

5 rows × 31 columns

Coordinates

[10]:

cut_gal_coord = SkyCoord(ra=cut_galaxies.RAJ2000, dec=cut_galaxies.DEJ2000, unit='deg')

[11]:

ngc_4993 = SkyCoord('13h09m47.706s -23d23m01.79s', frame='icrs')

[12]:

np.min(ngc_4993.separation(cut_gal_coord).to('arcmin'))

[12]:
$6.96812e-05\mathrm{{}^{\prime}}$
[13]:

np.argmin(ngc_4993.separation(cut_gal_coord))

[13]:

11

[14]:

cut_galaxies.iloc[11]

[14]:

RAJ2000      197.448776
DEJ2000      -23.383831
Dist          33.806000
Bmag          12.870000
BmagHyp       13.450000
ImagHyp       11.310000
modzHyp       33.130000
mod0Hyp             NaN
logd25Hyp      1.190000
logr25Hyp      0.060000
logdcHyp       1.260000
PAHyp        173.200000
BmagHypC      12.870000
ImagHypC      11.080000
U-BHypC             NaN
B-VHypC             NaN
Jmag2         10.292000
Hmag2          9.599000
Kmag2          9.330000
a_b            0.900000
PAK          -20.000000
BmagGWGC      12.870000
maj            0.932000
min            0.792000
PAGWGC       173.200000
DistGWGC      33.806000
Kmag2MPZ       9.285000
Bmag2MPZ      13.160000
zsp2MPZ        0.009680
zph2MPZ        0.012307
Flag           0.000000
Name: 667145, dtype: float64

PATH time

Instantiate

[17]:

Path = path.PATH()

Candidates

[19]:

Path.init_candidates(cut_galaxies.RAJ2000.values,
                     cut_galaxies.DEJ2000.values,
                     cut_galaxies.maj.values,
                     mag=cut_galaxies.Bmag.values)

Priors

[20]:

# Candidates
Path.init_cand_prior('inverse', P_U=0.)
# Offsets
Path.init_theta_prior('exp', 6.)

Localization

[21]:

Path.init_localization('healpix',
               healpix_data=gw170817,
               healpix_nside=header['NSIDE'],
               healpix_ordering='NESTED',
               healpix_coord='C')

Calculate Priors

[23]:

P_O = Path.calc_priors()

Calculate Posteriors

[24]:

P_Ox, P_Ux = Path.calc_posteriors('local', box_hwidth=30.)

[25]:

Path.candidates.sort_values('P_Ox', ascending=False)

[25]:
ra dec ang_size mag P_O p_xO P_Ox
11 197.448776 -23.383831 0.932 12.87 0.006478 2.783514e-04 9.999998e-01
144 193.861099 -20.892786 0.776 15.73 0.000128 3.112096e-09 2.201734e-07
64 194.469193 -22.638481 1.047 15.52 0.000167 3.205848e-11 2.972815e-09
347 206.735901 -39.884365 0.757 14.88 0.000388 2.893990e-16 6.227795e-14
297 192.973251 -12.604387 0.468 16.27 0.000064 1.294926e-17 4.627079e-16
... ... ... ... ... ... ... ...
132 180.701767 21.145529 0.871 16.00 0.000090 0.000000e+00 0.000000e+00
131 110.026970 53.010151 0.757 14.28 0.000875 0.000000e+00 0.000000e+00
129 140.950027 2.112878 0.630 15.07 0.000301 0.000000e+00 0.000000e+00
128 21.381145 32.136272 1.582 13.63 0.002167 0.000000e+00 0.000000e+00
389 171.504303 1.984145 0.775 14.04 0.001219 0.000000e+00 0.000000e+00

390 rows × 7 columns

F’ing A

Multi order

https://emfollow.docs.ligo.org/userguide/tutorial/multiorder_skymaps.html

[28]:

mofile = os.path.join(resource_filename('astropath', 'data'), 'gw_examples',
                     'GW190814_PublicationSamples.multiorder.fits')