Commit 0011902a authored by Tiffany Summerscales's avatar Tiffany Summerscales
Browse files

Replaced broken trunk version of megasky.py with working version from tags/O2_online

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@737 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent e11fbf1d
#!/usr/bin/env python
import lal
import matplotlib
matplotlib.use('Agg')
import numpy as np
......@@ -7,7 +8,8 @@ import glob
import os
import sys
from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS
from pylal.xlal.date import XLALGreenwichSiderealTime
#from pylal.xlal.date import XLALGreenwichSiderealTime
from lal import GreenwichSiderealTime
import find_trig as ft
import healpy as hp
from subprocess import call
......@@ -17,10 +19,91 @@ from pylal import bayespputils
import sky_area.sky_area_clustering as sky
import getopt
from readbwb import BwbParams
import math
import lalinference.fits as lf
#import lalinference.fits as lf
import lalinference.io.fits as lf
import lalinference.plot as lp
import lalinference.cmap
import lalinference.plot.cmap
#import lalinference.cmap
def gps_to_mjd(gps_time):
"""
Convert a floating-point GPS time in seconds to a modified Julian day.
Parameters
----------
gps_time : float
Time in seconds since GPS epoch
Returns
-------
mjd : float
Modified Julian day
Example
-------
>>> '%.9f' % round(gps_to_mjd(1129501781.2), 9)
'57316.937085648'
"""
gps_seconds_fraction, gps_seconds = math.modf(gps_time)
mjd = lal.ConvertCivilTimeToMJD(lal.GPSToUTC(int(gps_seconds)))
return mjd + gps_seconds_fraction / 86400.
def gps_to_iso8601(gps_time):
"""
Convert a floating-point GPS time in seconds to an ISO 8601 date string.
Parameters
----------
gps : float
Time in seconds since GPS epoch
Returns
-------
iso8601 : str
ISO 8601 date string (with fractional seconds)
Example
-------
>>> gps_to_iso8601(1000000000.01)
'2011-09-14T01:46:25.010000'
>>> gps_to_iso8601(1000000000)
'2011-09-14T01:46:25.000000'
>>> gps_to_iso8601(1000000000.999999)
'2011-09-14T01:46:25.999999'
>>> gps_to_iso8601(1000000000.9999999)
'2011-09-14T01:46:26.000000'
>>> gps_to_iso8601(1000000814.999999)
'2011-09-14T01:59:59.999999'
>>> gps_to_iso8601(1000000814.9999999)
'2011-09-14T02:00:00.000000'
"""
gps_seconds_fraction, gps_seconds = math.modf(gps_time)
gps_seconds_fraction = '{0:6f}'.format(gps_seconds_fraction)
if gps_seconds_fraction[0] == '1':
gps_seconds += 1
else:
assert gps_seconds_fraction[0] == '0'
assert gps_seconds_fraction[1] == '.'
gps_seconds_fraction = gps_seconds_fraction[1:]
year, month, day, hour, minute, second, _, _, _ = lal.GPSToUTC(int(gps_seconds))
return '{0:04d}-{1:02d}-{2:02d}T{3:02d}:{4:02d}:{5:02d}{6:s}'.format(year, month, day, hour, minute, second, gps_seconds_fraction)
# -------------------------------------------
# Module to make skymaps, skyview webpage, etc.
......@@ -45,12 +128,18 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
os.chdir(directory)
print "Entered", os.getcwd()
jobname = glob.glob('*bayeswave.run')[0]
try:
jobname = 'bayeswave.run'
bayeswave = open(jobname,'r')
except:
sys.exit("\n *bayeswave.run not found! \n")
# ---------------------------------------
# Extract information from bayeswave.run
# Note: We may only need the trigger time from this
# ---------------------------------------
bayeswave = open(jobname, 'r')
#bayeswave = open(jobname, 'r')
ifoNames = []
for line in bayeswave:
spl = line.split()
......@@ -116,7 +205,8 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
print "Using GPS {0}".format(gps)
trigtime = float(gps)
lalgps = LIGOTimeGPS(trigtime)
sidtime = XLALGreenwichSiderealTime(lalgps, 0)
#sidtime = XLALGreenwichSiderealTime(lalgps, 0)
sidtime = GreenwichSiderealTime(lalgps, 0)
sidtime = sidtime % (np.pi*2)
# -- Get the injection location
......@@ -134,6 +224,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
print "GOT INJECTION PARAMTERS"
print injtheta, injra
# -- Special handling for injections drawn from prior
if params.ra is not None:
injdec = params.dec
......@@ -169,9 +260,11 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
rarad = rarad + 2*np.pi
# Shift the declination to match hp.projscatter conventions
decrad = np.pi*0.5 - decrad
# -- Plot the skymap and injection location
skymap = kde.as_healpix(NSIDE, nest=True)
#fig = plt.figure(frameon=False)
fig = plt.figure(figsize=(8,6), frameon=False)
......@@ -179,7 +272,9 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
ax.cla()
ax.grid()
lp.healpix_heatmap(skymap, nest=True, vmin=0.0, vmax=np.max(skymap), cmap=plt.get_cmap('cylon'))
# TODO: our .png skymaps don't work anymore.
if (ra is not None and dec is not None):
plt.plot(rarad, dec, 'kx', ms=30, mew=1)
if mdc is not None:
......@@ -187,7 +282,18 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
plt.savefig(plotsDir+'/skymap.png')
plt.close()
lf.write_sky_map('skymap_{0}.fits'.format(gps), skymap, nest=True, gps_time=trigtime)
# lf.write_sky_map('skymap_{0}.fits'.format(gps), skymap, nest=True, 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
......@@ -207,9 +313,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
else:
injcontour = 0
searcharea = 0
pixarea = hp.nside2pixarea(NSIDE, degrees=True)
print("here")
if results is not None:
for name in ['injcontour', 'area50', 'area90', 'searcharea']:
if name not in results: results[name] = []
......@@ -226,16 +330,25 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
# -- Write main script for command line running
if __name__ == "__main__":
# Example command line function calls
# megasky.py /path/to/working/dir --mdc=/path/to/mdc/mdclog.txt
# megasky.py --mdc=/path/to/mdc/mdclog.txt
# Allow navigation into specified working directory
topdir=os.getcwd()
try:
workdir=sys.argv[1]
except IndexError:
# No work directory specified, workdir=./
workdir=os.getcwd()
os.chdir(workdir)
opts, args = getopt.getopt(sys.argv[1:], "", ['directory=', 'mdc=', 'NSIDE=', 'ra=', 'dec=', 'geo'])
# Working directory is current directory unless it is specified by the
# first argument of the function call
workdir=os.getcwd()
inargs = sys.argv[1:]
if len(sys.argv) >= 2:
if os.path.isdir(sys.argv[1]):
workdir = sys.argv[1]
inargs = sys.argv[2:]
os.chdir(workdir)
opts, args = getopt.getopt(inargs, "", ['directory=', 'mdc=', 'NSIDE=', 'ra=', 'dec=', 'geo'])
# -- Set default argument values
directory='.'
......@@ -275,7 +388,7 @@ if __name__ == "__main__":
dec = float(arg)
if opt == '--geo':
geo = 'True'
make_skyview(directory, mdc, NSIDE, ra, dec, results)
# 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