Maintenance will be performed on and on the 25th April 2019 between 9-10pm CDT, towards the end of the maintenance window there will be a short (5-10 minute) downtime.

Commit f2833be0 authored by Sudarshan Ghonge's avatar Sudarshan Ghonge

Add sky area calculation and option of geocentric skymap plot

parent 930cbbbf
......@@ -41,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()
......@@ -139,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
......@@ -257,7 +258,11 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
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')
ax.imshow_hpx((probperdeg2, 'ICRS'), nested=True, vmin=0., vmax=vmax, cmap='cylon')
cls = 100 * postprocess.find_greedy_credible_levels(rhmap['PROB'].data)
......@@ -266,11 +271,18 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
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):
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:
......@@ -283,8 +295,32 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
rhmap.meta['gps_time'] = trigtime
io.write_sky_map('skymap_{0}.fits'.format(gps), rhmap, nest=True)
print("Sky Area calculation currently not supported")
# -- Get the found contour and searched areas
if mdc is not None:
print("Calculating p value of injection")
true_ra = 180/np.pi*injra
true_dec = 180/np.pi*injdec
elif (inj is not None):
print("Calculating p value of command line injection")
true_ra = 180/np.pi*inj['ra']
true_dec = 180/np.pi*inj['dec']
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])
outfile = open('skystats.txt', 'w')
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__":
......@@ -302,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
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