megasky.py 10.2 KB
Newer Older
1
#!/usr/bin/env python
Duncan Macleod's avatar
Duncan Macleod committed
2
3
4

from __future__ import print_function

5
import lal
Jonah Kanner's avatar
Jonah Kanner committed
6
7
8
9
10
11
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
12
from lal import GreenwichSiderealTime
Jonah Kanner's avatar
Jonah Kanner committed
13
import healpy as hp
14
15
16
from ligo.skymap import kde, plot, bayestar, postprocess, io
from astropy.coordinates import SkyCoord
from astropy.time import Time
Jonah Kanner's avatar
Jonah Kanner committed
17
import getopt
18
import math
19
20
21
22
23
24
25
26
27
28
29
30
31
from optparse import OptionParser

from glue.ligolw import ligolw
from glue.ligolw import lsctables
from glue.ligolw import table
from glue.ligolw import utils

class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
    pass


lsctables.use_in(LIGOLWContentHandler)

Jonah Kanner's avatar
Jonah Kanner committed
32

33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
def parser():
    """
    Parser for input (command line and ini file)
    """

    # --- cmd line
    parser = OptionParser()
    parser.add_option("-M", "--mdc", type=str, default=None)
    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")
    
    (opts,args) = parser.parse_args()

    if len(args)==0:
Duncan Macleod's avatar
Duncan Macleod committed
51
        print("ERROR: require BW directory", file=sys.stderr)
52
53
        sys.exit()
    if not os.path.isdir(args[0]):
Duncan Macleod's avatar
Duncan Macleod committed
54
        print("ERROR: BW dir %s does not exist"%args[0], file=sys.stderr)
55
        sys.exit()
56
57


58
    return opts, args
59

60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
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)




Tyson Littenberg's avatar
Tyson Littenberg committed
136

Jonah Kanner's avatar
Jonah Kanner committed
137
138
139
140
141
# -------------------------------------------
# Module to make skymaps, skyview webpage, etc.
# Works with single output directory from 
# BayesWaveBurst
# -------------------------------------------
142
def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
Jonah Kanner's avatar
Jonah Kanner committed
143
144
145
146
147
148
149

    # -- Close any open figures
    plt.close('all')
    # ----------------
    # Read S6 MDC log file
    # -----------------
    if mdc is None:
Duncan Macleod's avatar
Duncan Macleod committed
150
        print("No MDC log provided.")
Jonah Kanner's avatar
Jonah Kanner committed
151
152
153
154
155
156

    # -- Save PWD for future reference
    topdir = os.getcwd()

    # -- Enter requested directory
    os.chdir(directory)
Duncan Macleod's avatar
Duncan Macleod committed
157
    print("Entered", os.getcwd())
158
159
160
161
162
163

    try:
        jobname = 'bayeswave.run'
        bayeswave = open(jobname,'r')
    except:
        sys.exit("\n *bayeswave.run not found! \n")
Jonah Kanner's avatar
Jonah Kanner committed
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
    
    # ---------------------------------------
    # Extract information from bayeswave.run
    # Note: We may only need the trigger time from this
    # ---------------------------------------
    ifoNames = []
    for line in bayeswave:
        spl = line.split()
        if len(spl) < 1: continue
        
        # Determine names of IFOs and MDC scale factor
        if spl[0] == 'Command' and spl[1] == 'line:':
            for index, splElement in enumerate(spl):
                if splElement == '--ifo':
                    ifoNames.append(spl[index+1]) 
                    
        del spl[-1]
        # Determine number of IFOs
        if ' '.join(spl) == 'number of detectors':
            spl = line.split()
            ifoNum = spl[3]
            continue         
        # Determine gps trigger time
        if ' '.join(spl) == 'trigger time (s)':
            spl = line.split()
            gps = spl[3]
            break
    bayeswave.close()

    # --------------------------------
    # Create skymap summary statistics
    # --------------------------------
Jonah Kanner's avatar
Jonah Kanner committed
196

Jonah Kanner's avatar
Jonah Kanner committed
197
    # -- Input skymap data
Duncan Macleod's avatar
Duncan Macleod committed
198
    print("Extracting RA/DEC samples")
199
    filename = './chains/' + 'signal_params_h0.dat.0'
Tyson Littenberg's avatar
Tyson Littenberg committed
200
    data = np.loadtxt(filename, unpack=True,usecols=(0,1,2))
201
202
    ralist = data[1]
    sin_dec = data[2]
Duncan Macleod's avatar
Duncan Macleod committed
203
    print("Total samples are {0}".format(ralist.size))
Jonah Kanner's avatar
Jonah Kanner committed
204
205

    # -- Remove burn in samples
206
    burnin = int(ralist.size/4)
Jonah Kanner's avatar
Jonah Kanner committed
207
208
    ralist = ralist[burnin:]
    sin_dec = sin_dec[burnin:]
Duncan Macleod's avatar
Duncan Macleod committed
209
    print("After removing burn-in samples are {0}".format(ralist.size))
Jonah Kanner's avatar
Jonah Kanner committed
210

Jonah Kanner's avatar
Jonah Kanner committed
211
212
213
214
215
    declist = np.arcsin(sin_dec)
    thetalist = np.pi/2.0 - declist

    radec = np.column_stack((ralist, declist))

Tyson Littenberg's avatar
Tyson Littenberg committed
216
217
218
    if radec.shape[0] > npost:
        radec = np.random.permutation(radec)[:npost,:]

219
    skypost = kde.Clustered2DSkyKDE(radec)
Jonah Kanner's avatar
Jonah Kanner committed
220
221

    # -- Get the sidereal time
Duncan Macleod's avatar
Duncan Macleod committed
222
223
    print("Finding sidereal time")
    print("Using GPS {0}".format(gps))
Jonah Kanner's avatar
Jonah Kanner committed
224
    trigtime = float(gps)
225
226
    lalgps = lal.LIGOTimeGPS(trigtime)
    sidtime = GreenwichSiderealTime(lalgps, 0)
Jonah Kanner's avatar
Jonah Kanner committed
227
228
229
    sidtime = sidtime % (np.pi*2)

    # -- Get the injection location
Duncan Macleod's avatar
Duncan Macleod committed
230
231
    print("GOT MDC?")
    print(mdc)
Jonah Kanner's avatar
Jonah Kanner committed
232
233
234
235
236
237
238
239
240
    if mdc is None:
        injtheta = 0
        injphi   = 0
        injdec   = 0
        injra    = 0
    else:
        injtheta, injphi = mdc.get_theta_phi(trigtime)
        injra = injphi + sidtime
        injdec = np.pi/2 - injtheta
Duncan Macleod's avatar
Duncan Macleod committed
241
242
        print("GOT INJECTION PARAMTERS")
        print(injtheta, injra)
Jonah Kanner's avatar
Jonah Kanner committed
243
244

    # -- Make plots directory, if needed
245
246
247
248
    plotsDir = './plots'
    if not os.path.exists(plotsDir):
        os.makedirs(plotsDir)

249
 
Jonah Kanner's avatar
Jonah Kanner committed
250
    # -- Plot the skymap and injection location
251
252
253
254
255
256
    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()
257

Tyson Littenberg's avatar
Tyson Littenberg committed
258
259
   
    fig = plt.figure(figsize=(8,6), frameon=False)
260
    ax = plt.subplot(111, projection='astro hours mollweide')
Tyson Littenberg's avatar
Tyson Littenberg committed
261
    ax.grid()
262
263
264
265
266
267
    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])


268
269
270


    # TODO: our .png skymaps don't work anymore. 
271
    if (inj is not None):
272
        ax.plot_coord(Skycoord(np.rad2drg(inj['ra']), np.rad2deg(inj['dec']), unit='deg'), '*',markerfacecolor='white', markeredgecolor='black', markersize=10)
Jonah Kanner's avatar
Jonah Kanner committed
273
    if mdc is not None:
274
        ax.plot_coord(Skycoord(np.rad2drg(injra), np.rad2deg(injdec), unit='deg'), '*', markerfacecolor='white', markeredgecolor='black', markersize=10)
275
    plt.savefig(plotsDir+'/skymap.png')
Jonah Kanner's avatar
Jonah Kanner committed
276
277
    plt.close()
    
278
279
280
    rhmap.meta['origin'] = 'LIGO/Virgo'
    rhmap.meta['gps_creation_time'] = Time.now().gps
    rhmap.meta['gps_time'] = trigtime
281

282
283
    io.write_sky_map('skymap_{0}.fits'.format(gps), rhmap, nest=True)
    print("Sky Area calculation currently not supported")
Jonah Kanner's avatar
Jonah Kanner committed
284
285
286
287
288


# -- Write main script for command line running
if __name__ == "__main__":

289
290
291
292
    # 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

293
294
295
    # Allow navigation into specified working directory
    topdir=os.getcwd()

296
297
    # Working directory is current directory unless it is specified by the 
    # first argument of the function call
298
    opts, args = parser()
299
         
300
301
302
303
304
305
306
    ra = opts.ra
    dec = opts.dec
    geo = opts.dec
    mdc = opts.mdc
    NSIDE = opts.nside

    injpos = None
307
    # Checl if injection file is present
308
    if opts.inj is not None:
309
310
      # If present, check if event id is provided
      # If yes, proceed. If no, throw error message and exit
311
      if(opts.eventnum is None):
Duncan Macleod's avatar
Duncan Macleod committed
312
        print("Provide event num if giving injfile", file=sys.stderr)
313
        sys.exit()
Duncan Macleod's avatar
Duncan Macleod committed
314
      print("Loading xml")    
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
      xmldoc = utils.load_filename(
           opts.inj, contenthandler=LIGOLWContentHandler)
      try:
        print('Checking if using a sim_inspiral table...')
        injs = table.get_table(
            xmldoc, lsctables.SimInspiralTable.tableName)
        inj = injs[opts.eventnum]
        injpos = {
            'ra': inj.longitude,
            'dec': inj.latitude,
            'id': inj.simulation_id}
        print(' yes')
      except:
        print('Checking if using a sim_burst table...')
        injs = table.get_table(xmldoc, lsctables.SimBurstTable.tableName)
        inj = injs[opts.eventnum]
        injpos = {'ra': inj.ra, 'dec': inj.dec, 'id': inj.simulation_id}
        print(' yes')
Jonah Kanner's avatar
Jonah Kanner committed
333
    
334
    # If ra and dec provided manually, use those
335
336
    if ra is not None and dec is not None:
      injpos  ={'ra': ra, 'dec':dec, 'id':0}
337
338
    # If md log provided, use that. Currently,
    # there is no scheme to assert exclusivity of injection and mdcs
339
    if mdc is not None:
Duncan Macleod's avatar
Duncan Macleod committed
340
      print('Reading mdc log {0}'.format(mdc))
341
342
343
      try:
        mdc = ft.Mdc(arg)
      except:
Duncan Macleod's avatar
Duncan Macleod committed
344
        print('Warning! Failed to read mdc')
345
346
        mdc = None

347
    make_skyview(args[0], mdc, NSIDE, injpos)
Jonah Kanner's avatar
Jonah Kanner committed
348

349
350
351
352
    # Move back to original dir
    os.chdir(topdir)