Commit 41f98216 authored by Sudarshan Ghonge's avatar Sudarshan Ghonge
Browse files

Update megsky with ligo.skymap

parent dae3f100
......@@ -7,15 +7,14 @@ import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import glob
import os
import sys
from lal import GreenwichSiderealTime
from bayeswave_plot import find_trig as ft
import healpy as hp
import sky_area.sky_area_clustering as sky
from ligo.skymap import kde, plot, bayestar, postprocess, io
from astropy.coordinates import SkyCoord
from astropy.time import Time
import getopt
from bayeswave_plot.readbwb import BwbParams
import math
from optparse import OptionParser
......@@ -31,9 +30,6 @@ class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
import as lf
import lalinference.plot as lp
import lalinference.plot.cmap
def parser():
Parser for input (command line and ini file)
......@@ -59,7 +55,7 @@ def parser():
return opts, args,
return opts, args
def gps_to_mjd(gps_time):
......@@ -170,7 +166,6 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
# Extract information from
# Note: We may only need the trigger time from this
# ---------------------------------------
#bayeswave = open(jobname, 'r')
ifoNames = []
for line in bayeswave:
spl = line.split()
......@@ -208,7 +203,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
print("Total samples are {0}".format(ralist.size))
# -- Remove burn in samples
burnin = ralist.size/4
burnin = int(ralist.size/4)
ralist = ralist[burnin:]
sin_dec = sin_dec[burnin:]
print("After removing burn-in samples are {0}".format(ralist.size))
......@@ -221,7 +216,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
if radec.shape[0] > npost:
radec = np.random.permutation(radec)[:npost,:]
kde = sky.ClusteredSkyKDEPosterior(radec)
skypost = kde.Clustered2DSkyKDE(radec)
# -- Get the sidereal time
print("Finding sidereal time")
......@@ -253,71 +248,39 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
# -- Plot the skymap and injection location
skymap = kde.as_healpix(NSIDE, nest=True)
skymap = skypost.as_healpix()
rhmap = bayestar.rasterize(skymap)
nside = hp.npix2nside(len(rhmap))
deg2perpix = hp.nside2pixarea(nside, degrees=True)
probperdeg2 = rhmap['PROB']/deg2perpix
vmax = probperdeg2.max()
fig = plt.figure(figsize=(8,6), frameon=False)
ax = plt.subplot(111, projection='astro hours mollweide')
lp.healpix_heatmap(skymap, nest=True, vmin=0.0, vmax=np.max(skymap), cmap=plt.get_cmap('cylon'))
ax.imshow_hpx((probperdeg2, 'ICRS'), nested=True, vmin=0., vmax=vmax, cmap='cylon')
cls = 100 * postprocess.find_greedy_credible_levels(rhmap['PROB'].data)
cs = ax.contour_hpx((cls, 'ICRS'), nested=True,
colors='k', linewidths=0.5, levels=[50,90])
# TODO: our .png skymaps don't work anymore.
if (inj is not None):
plt.plot(inj['ra'], inj['dec'], 'kx', ms=30, mew=1)
ax.plot_coord(Skycoord(np.rad2drg(inj['ra']), np.rad2deg(inj['dec']), unit='deg'), '*',markerfacecolor='white', markeredgecolor='black', markersize=10)
if mdc is not None:
plt.plot(injra, injdec, 'kx', ms=30, mew=1)
ax.plot_coord(Skycoord(np.rad2drg(injra), np.rad2deg(injdec), unit='deg'), '*', markerfacecolor='white', markeredgecolor='black', markersize=10)
# Meg's hacky way to get sky maps to work, because this: lf.write_sky_map('skymap_{0}.fits'.format(gps), skymap, nest=True, gps_time=trigtime) doesn't work anymore
skymap = np.atleast_2d(skymap)
extra_header = []
extra_header.append(('DATE-OBS', gps_to_iso8601(trigtime), 'UTC date of the observation'))
extra_header.append(('MJD-OBS', gps_to_mjd(trigtime), 'modified Julian date of the observation'))
hp.write_map('skymap_{0}.fits'.format(gps), skymap, nest=True, fits_IDL=True, coord='C', column_names=('PROB', 'DISTMU', 'DISTSIGMA', 'DISTNORM')[:len(skymap)], column_units=('pix-1', 'Mpc', 'Mpc', 'Mpc-2')[:len(skymap)], extra_header=extra_header)
# -- Calculate the 50 and 90% credible intervals
sq_deg = (180.0/np.pi)**2
print("Calculating sky area ...")
(area50, area90) = kde.sky_area( [0.5, 0.9] )*sq_deg
print("Got area50 of {0}".format(area50))
# -- Get the found contour and searched areas
if mdc is not None:
print("Calculating p value of injection")
injcontour = kde.p_values( np.array([[injra, injdec]]) )[0]
print("Got contour of {0}".format(injcontour))
print("Calculating searched area...")
searcharea = (kde.searched_area( np.array([[injra, injdec]]) )*sq_deg)[0]
print("Got searched area of {0}".format(searcharea))
elif (inj is not None):
print("Calculating p value of command line injection")
injcontour = kde.p_values( np.array([inj['ra'], inj['dec']]) )[0]
print("Got contour of {0}".format(injcontour))
rhmap.meta['origin'] = 'LIGO/Virgo'
rhmap.meta['gps_creation_time'] =
rhmap.meta['gps_time'] = trigtime
print("Calculating searched area...")
searcharea = (kde.searched_area( np.array([inj['ra'], inj['dec']]) )*sq_deg)[0]
print("Got searched area of {0}".format(searcharea))
injcontour = 0
searcharea = 0
pixarea = hp.nside2pixarea(NSIDE, degrees=True)
# -- Make an output file with key statistics
outfile = open('skystats.txt', 'w')
outfile.write("# area50 area90 searcharea injcontour \n")
outfile.write("{0} {1} {2} {3}".format(area50, area90, searcharea, injcontour))
io.write_sky_map('skymap_{0}.fits'.format(gps), rhmap, nest=True)
print("Sky Area calculation currently not supported")
# -- Write main script for command line running
......@@ -381,7 +344,6 @@ if __name__ == "__main__":
print('Warning! Failed to read mdc')
mdc = None
make_skyview(args[0], mdc, NSIDE, injpos)
# Move back to original dir
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment