<jupyter_start><jupyter_text>LVC GW TargetsGenerate a target list of host galaxies given the 50% and 90% credible localization regions of a GW alert.<jupyter_code>from astropy.io import fits
from astropy.table import Table

import healpy as hp

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

from gw_skymap import compute_quantiles, compute_contours, plot_gwmap

from desitarget.io import read_targets_in_hp, read_targets_in_box, read_targets_in_cap, check_hp_target_dir
from desitarget.targetmask import desi_mask, bgs_mask<jupyter_output><empty_output><jupyter_text>Read GW LALInference MapExtract localization PDF for a GW event and plot it using the DESI footprint.Data can be downloaded from the [LVC GraceDB public alert page](https://gracedb.ligo.org/).<jupyter_code>fits_file = 'GW190412_skymap.fits.gz'

# Read metadata.
hdus = fits.open(fits_file)
header = hdus[1].header

# instruments = header['INSTRUME']
distmean = header['DISTMEAN']
diststd = header['DISTSTD']
origin = header['ORIGIN']
date = header['DATE']

gwmap = hp.read_map(fits_file)
npix = len(gwmap)
nside = hp.npix2nside(npix)
# Compute GW quantiles.
levels = [0.50, 0.90]
lev50, lev90 = compute_quantiles(levels, gwmap)
print(levels, lev50, lev90)

# # Compute GW quantiles for degraded map.
# gwmap64 = hp.pixelfunc.ud_grade(gwmap, 64) #reduce nside to make it faster
# gwmap64 = gwmap64/np.sum(gwmap64)
# lev50, lev90 = compute_quantiles(levels, gwmap64)
# print(levels, lev50, lev90)
fig = plot_gwmap('GW190412_skymap.fits.gz', levels=[0.5, 0.9])
fig.savefig('gw_skymap.png', dpi=100)
maxpix = np.argmax(gwmap)
ra_c, dec_c = hp.pix2ang(nside, maxpix, lonlat=True)
print(ra_c, dec_c)<jupyter_output>218.49609375 36.469945487867896<jupyter_text>Load DR8 TargetsUsing the code from the [DESI target example notebook](https://github.com/qmxp55/desi_uk_meeting_2019/blob/master/desitarget_local_example.ipynb) created by <NAME>, extract the targets within a box surrounding the 90% C.I.<jupyter_code># pixids = hp.ring2nest(64, np.arange(len(gwmap64)))
# pix50 = pixids[gwmap64 > lev50]
# pix90 = pixids[gwmap64 > lev90]

pixids = hp.ring2nest(nside, np.arange(npix))
pix50 = pixids[gwmap > lev50]
pix90 = pixids[gwmap > lev90]
# Load sample of resolved DR8 target catalog using HEALPixels above threshold in the GW map.
hpdirname = '/global/project/projectdirs/desi/target/catalogs/dr9m/0.42.0/targets/main/resolve/dark'
targ90 = Table(read_targets_in_hp(hpdirname, nside=nside, pixlist=pix90))
len(targ90)
bgs_mask_bits = '|'.join([_ for _ in bgs_mask.names() if 'BGS' in _])
bgs_mask_bits<jupyter_output><empty_output><jupyter_text>Select TargetsIdentify targets within the 50% and 90% contours, and apply a basic cut on bright BGS objects.<jupyter_code>bgs_targets = (targ90['BGS_TARGET'] & bgs_mask.mask(bgs_mask_bits)) != 0
ra90, dec90 = [targ90[_][bgs_targets] for _ in ['RA', 'DEC']]

rmag = 22.5 - 2.5*np.log10(targ90['FLUX_R'][bgs_targets])
rthr = np.max(rmag)
print('{:8d} targets in 90% C.I. with rmag > {:g}'.format(len(ra90), rthr))<jupyter_output>6845 targets in 90% C.I. with rmag > 20.0333<jupyter_text>Plot the Targets<jupyter_code># fig = plt.figure(1, figsize=(4,6))

hp.cartview(gwmap, fig=1, cmap='OrRd',
            lonra=[ra_c-360-3,ra_c-360+3],
            latra=[dec_c-3,dec_c+3],
            min=0, max=np.max(gwmap),
            title='Targets: GW190412', unit=r'Pr$(\alpha,\delta)$', flip='astro')
# hp.projplot(ra50, dec50, 'k,', alpha=0.1, lonlat=True);
hp.projplot(ra90, dec90, 'k,', alpha=0.3, lonlat=True);

# Compute GW contours for degraded map.
# gwmap64 = hp.pixelfunc.ud_grade(gwmap, 64) #reduce nside to make it faster
# gwmap64 = gwmap64/np.sum(gwmap64)
# ra_contour, dec_contour = compute_contours(levels, gwmap64)
gwmap256 = hp.pixelfunc.ud_grade(gwmap, 256) #reduce nside to make it faster
gwmap256 = gwmap256/np.sum(gwmap256)
ra_contour, dec_contour = compute_contours(levels, gwmap256)
for i, (rc, dc) in enumerate(zip(ra_contour, dec_contour)):
    hp.projplot(rc, dc, lonlat=True, linewidth=1, c='g')
    
# fig = plt.gcf()
# fig.savefig('targlist.png', dpi=150)
targ90[bgs_targets]
targ90[bgs_targets].write('GW190412_BGS_targets_dr9m_0.42.0_dark.fits', format='fits', overwrite=True)
def angle_to_180(angle):
    """Normalize longitude angle between -180 and 180."""
    a = angle
    while a < -180:
        a += 360
    while a > 180:
        a -= 360
    return a

xmin = np.round(ra_c - 3)
xmax = np.round(ra_c + 3)
# xmin = angle_to_180(np.round(ra_c - 5))
# xmax = angle_to_180(np.round(ra_c + 5))
if xmax < xmin:
    xmin, xmax = xmax, xmin
cxmin, cxmax = xmin, xmax
frot = 0.
if xmax > 90 and xmin < -90:
    frot, cxmin, cmax = 180., xmax-180., xmax+180.
ymin = np.round(dec_c - 3)
ymax = np.round(dec_c + 3)

cxmin, cxmax, ymin, ymax

faspect = np.abs(cxmax - cxmin)/np.abs(ymax-ymin)
fysize = 4
figsize = (fysize*faspect+2, fysize+2.75)

tfig   = plt.figure(num=2,figsize=figsize)
rotimg = hp.cartview(gwmap, fig=2,coord='C', title="", cbar=False, flip='astro',
                     lonra=[cxmin,cxmax], latra=[ymin,ymax], rot=frot,
                     notext=True, xsize=1000,
                     return_projected_map=True)
plt.close(tfig)
fig, ax = plt.subplots(1,1, num=1, figsize=figsize)
img = ax.imshow(rotimg, extent=[cxmax, cxmin, ymin, ymax],
                origin='lower', cmap='OrRd')

for i, (rc, dc, lstyle, clev) in enumerate(zip(ra_contour, dec_contour, ['--', '-'], ['50', '90'])):
    p = ax.plot(rc, dc, 'g-', ls=lstyle, lw=1, label='{}% CI'.format(clev))

ratarg, dectarg = ra90, dec90
ratarg[ratarg < 0] += 360

ax.plot(ratarg, dectarg, 'k,', alpha=0.3)
ax.set(xlim=(cxmax, cxmin),
       xlabel='RA [deg]',
       ylabel='Dec [deg]')
ax.grid(ls=':')

circ = plt.Circle((ra_c, dec_c), radius=1.6, fc='None', ec='k', ls='--')
ax.add_artist(circ)

_h, _l = ax.get_legend_handles_labels()
_h.append(circ)
_l.append('DESI FOV')

ax.legend(handles=_h, labels=_l)

cb = fig.colorbar(img, orientation='horizontal', shrink=0.85, fraction=0.1,
                  pad=0.1, ax=ax)
cb.set_label(r'GW angular posterior Pr$(\alpha,\delta)$')
fig.savefig('targlist.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14938431626701232, 'fraction_numerical': 0.05524951393389501, 'mean_word_length': 4.8071495766698025, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 28, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '444581', 'license': 'MIT', 'n_tokens_mistral': 2796, 'n_tokens_neox': 2533, 'n_words': 537, 'path': 'gwtarget/targeting_gw.ipynb', 'repo_name': 'divijsharma18/timedomain', 'size': '6177'}