Commit d71b0b66 authored by James Clark's avatar James Clark

Merge branch 'megasky-python3-fix' into 'master'

Update megasky with ligo.skymap

See merge request !59
parents 6c68b87a a4901745
Pipeline #47193 passed with stages
in 3 minutes and 39 seconds
......@@ -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)
......@@ -45,9 +41,10 @@ def parser():
parser.add_option("-N", "--nside", type=int, default=128)
parser.add_option("-I", "--inj", type=str, default=None, help='Injection xml')
parser.add_option("-e", "--eventnum", type=int, default=None, help='Injection event id')
parser.add_option("-r", "--ra", type=float, default=None)
parser.add_option("-d", "--dec", type=float, default=None)
parser.add_option("--geo", default=False, action="store_true")
parser.add_option("-r", "--ra", type=float, default=None, help="Right Ascension of injection")
parser.add_option("-d", "--dec", type=float, default=None, help="Declination of injection")
parser.add_option("-n", "--npost", type=int, default=5000, help="Number of points to use in making skymap")
parser.add_option("--geo", default=False, action="store_true", help="Whether to project onto world map")
(opts,args) = parser.parse_args()
......@@ -59,7 +56,7 @@ def parser():
return opts, args,
return opts, args
def gps_to_mjd(gps_time):
......@@ -143,7 +140,7 @@ def gps_to_iso8601(gps_time):
# Works with single output directory from
# BayesWaveBurst
# -------------------------------------------
def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000, geo=False):
# -- Close any open figures
......@@ -170,7 +167,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 +204,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 +217,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,72 +249,78 @@ 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')
obstime = Time(trigtime, format='gps').utc.isot
ax = plt.subplot(111, projection='geo degrees mollweide', obstime=obstime)
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])
fmt = r'%g\%%' if matplotlib.rcParams['text.usetex'] else '%g%%'
plt.clabel(cs, fmt=fmt, fontsize=6, inline=True)
import json
geojson_filename = os.path.join(
os.path.dirname(plot.__file__), 'ne_simplified_coastline.json')
with open(geojson_filename, 'r') as geojson_file:
geoms = json.load(geojson_file)['geometries']
verts = [coord for geom in geoms
for coord in zip(*geom['coordinates'])]
plt.plot(*verts, color='0.5', linewidth=0.5,
# 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.rad2deg(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.rad2deg(injra), np.rad2deg(injdec), unit='deg'), '*', markerfacecolor='white', markeredgecolor='black', markersize=10)
rhmap.meta['origin'] = 'LIGO/Virgo'
rhmap.meta['gps_creation_time'] =
rhmap.meta['gps_time'] = trigtime
# 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))
io.write_sky_map('skymap_{0}.fits'.format(gps), rhmap, nest=True)
# -- 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))
true_ra = injra
true_dec = injdec
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))
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))
true_ra = inj['ra']
true_dec = inj['dec']
injcontour = 0
searcharea = 0
pixarea = hp.nside2pixarea(NSIDE, degrees=True)
true_ra = None
true_dec = None
( searched_area, searched_prob, offset, searched_modes, contour_areas,
area_probs, contour_modes, searched_prob_dist, contour_dists,
searched_vol, searched_prob_vol, contour_vols
) = postprocess.find_injection_moc(skymap,true_ra=true_ra, true_dec=true_dec, contours=[0.5, 0.9])
# -- 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))
outfile.write("# area50 area90 searcharea p_value \n")
outfile.write("{0} {1} {2} {3}".format(contour_areas[0], contour_areas[1], searched_area, searched_prob))
# -- Write main script for command line running
if __name__ == "__main__":
......@@ -336,9 +338,10 @@ if __name__ == "__main__":
ra = opts.ra
dec = opts.dec
geo = opts.dec
geo = opts.geo
mdc = opts.mdc
NSIDE = opts.nside
npost = opts.npost
injpos = None
# Checl if injection file is present
......@@ -381,9 +384,7 @@ if __name__ == "__main__":
print('Warning! Failed to read mdc')
mdc = None
make_skyview(args[0], mdc, NSIDE, injpos)
make_skyview(args[0], mdc, NSIDE, injpos, npost, geo)
# 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