another sky map fix

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/tags/O2_online@691 c56465c9-8126-4a4f-9d7d-ac845eff4865

postprocess/skymap/megasky.py

 #!/usr/bin/env python
+import lal
 import matplotlib
 matplotlib.use('Agg')
 import numpy as np
@@ -18,6 +19,7 @@ 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.io.fits as lf
@@ -25,6 +27,84 @@ import lalinference.plot as lp
 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.
 # Works with single output directory from 
@@ -194,7 +274,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
     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. The fits files work on gracedb though so ¯\_(ツ)_/¯	
+    # 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:
@@ -202,10 +282,16 @@ 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)
+
+
     # 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)