Merge pull request #65 from fact-project/mjd_axis

Mjd axis

examples/plot_qla.py

@@ -11,9 +11,8 @@
 qla_results = bin_runs(
     runs,
     bin_width_minutes=20,
-    discard_ontime_fraction=0.9,
 )
-ax1, ax2 = plot_excess_rate(qla_results)
+ax1, ax2, ax1_mjd = plot_excess_rate(qla_results)
 ax1.grid()
 
 plt.show()

fact/VERSION

@@ -1 +1 @@
-0.14.0
+0.14.0

fact/plotting/analysis.py

@@ -1,9 +1,66 @@
 import matplotlib.pyplot as plt
 from matplotlib.dates import DateFormatter
+import matplotlib.transforms as mtransforms
+from mpl_toolkits.axes_grid1.parasite_axes import SubplotHost
 import numpy as np
+from ..time import MJD_EPOCH
 
 
-def plot_excess_rate(binned_runs, outputfile=None):
+# Matplotlib uses ordinal for internal date represantion
+# to get from ordinal to mjd, we shift by the ordinal value
+# of the MJD_EPOCH
+MJD_AXES_TRANSFORM = (
+    mtransforms.Affine2D().translate(MJD_EPOCH.toordinal(), 0)
+)
+
+
+def create_datetime_mjd_axes(fig=None, *args, **kwargs):
+    '''
+    Create a plot with two x-axis, bottom axis using
+    dates, top axis using mjd.
+
+    Parameters
+    ----------
+    fig: matplotlib.Figure or None
+        the figure to use, if None use plt.gcf()
+
+    Returns
+    -------
+    ax: mpl_toolkits.axes_grid1.parasite_axes.SubplotHost
+        The ax for the dates
+    mjd_ax: mpl_toolkits.axes_grid1.parasite_axes.ParasiteAxis
+        The axis with the mjd axis
+
+    '''
+    if fig is None:
+        fig = plt.gcf()
+
+    if args == []:
+        ax = SubplotHost(fig, 1, 1, 1, **kwargs)
+    else:
+        ax = SubplotHost(fig, *args, **kwargs)
+
+    # The second axis shows MJD if the first axis uses dates
+    mjd_ax = ax.twin(MJD_AXES_TRANSFORM)
+    mjd_ax.set_viewlim_mode('transform')
+
+    # disable unwanted axes
+    mjd_ax.axis['right'].toggle(ticklabels=False, ticks=False)
+    mjd_ax.axis['bottom'].toggle(ticklabels=False, ticks=False)
+    mjd_ax.axis['bottom'].toggle(label=False)
+
+    # add/remove label
+    mjd_ax.axis['top'].set_label('MJD')
+
+    # Deactivate offset
+    mjd_ax.ticklabel_format(useOffset=False)
+
+    fig.add_subplot(ax)
+
+    return ax, mjd_ax
+
+
+def plot_excess_rate(binned_runs, outputfile=None, mjd=True):
     '''
     Create an excess rate plot from given data
 
@@ -23,14 +80,30 @@ def plot_excess_rate(binned_runs, outputfile=None):
     '''
     fig = plt.figure()
 
-    ax_sig = plt.subplot2grid((8, 1), (6, 0), rowspan=2)
-    ax_rate = plt.subplot2grid((8, 1), (1, 0), rowspan=5, sharex=ax_sig)
+    gridspec = plt.GridSpec(8, 1)
+    spec_sig = gridspec.new_subplotspec((6, 0), rowspan=2)
+    spec_rate = gridspec.new_subplotspec((1, 0), rowspan=5)
+
+    ax_sig = plt.subplot(spec_sig)
+
+    if mjd is True:
+        ax_rate, ax_rate_mjd = create_datetime_mjd_axes(
+            fig, spec_rate, sharex=ax_sig
+        )
+    else:
+        ax_rate = plt.subplot(spec_rate, sharex=ax_sig)
 
     colors = [e['color'] for e in plt.rcParams['axes.prop_cycle']]
 
     labels = []
     plots = []
-    for (name, group), color in zip(binned_runs.groupby('source'), colors):
+
+    groups = list(sorted(
+        binned_runs.groupby('source'),
+        key=lambda g: g[1].time_mean.min()
+    ))
+
+    for (name, group), color in zip(groups, colors):
         if len(group.index) == 0:
             continue
 
@@ -75,16 +148,25 @@ def plot_excess_rate(binned_runs, outputfile=None):
     ax_sig.set_yticks(np.arange(0, ymax + 0.1, ymax // 4 + 1))
 
     plt.setp(ax_rate.get_xticklabels(), visible=False)
+    plt.setp(ax_rate_mjd.get_xticklabels(), visible=False)
     plt.setp(ax_sig.get_xticklabels(), rotation=30, va='top', ha='right')
 
+    ax_rate_mjd.set_xlabel('')
+
     if binned_runs.night.nunique() <= 7:
         ax_sig.xaxis.set_major_formatter(DateFormatter('%Y-%m-%d %H:%M'))
     else:
         ax_sig.xaxis.set_major_formatter(DateFormatter('%Y-%m-%d'))
 
     fig.tight_layout()
 
+    if mjd is True:
+        plt.setp(ax_rate_mjd.get_xticklabels(), visible=True)
+
     if outputfile is not None:
         fig.savefig(outputfile)
 
-    return ax_rate, ax_sig
+    if mjd is True:
+        return ax_rate, ax_sig, ax_rate_mjd
+    else:
+        return ax_rate, ax_sig

fact/time.py

@@ -18,16 +18,62 @@
 from __future__ import print_function, division
 import time
 import calendar
-from datetime import datetime, timedelta
+from datetime import datetime, timedelta, timezone
 import logging
 
 import dateutil
 import dateutil.parser
+import numpy as np
 
 import pandas as pd
 
 OFFSET = (datetime(1970, 1, 1) - datetime(1, 1, 1)).days
 
+UNIX_EPOCH = datetime(1970, 1, 1, 0, 0, tzinfo=timezone.utc)
+MJD_EPOCH = datetime(1858, 11, 17, 0, tzinfo=timezone.utc)
+MJD_EPOCH_NUMPY = np.array(MJD_EPOCH.timestamp()).astype('datetime64[s]')
+
+
+def unixtime_to_mjd(unixtime):
+    return (unixtime - (MJD_EPOCH - UNIX_EPOCH).total_seconds()) / 3600 / 24
+
+
+def mjd_to_unixtime(mjd):
+    return (mjd + (MJD_EPOCH - UNIX_EPOCH).total_seconds()) * 3600 * 24
+
+
+def datetime_to_mjd(dt):
+    # handle numpy arrays
+    if isinstance(dt, np.ndarray):
+        mjd_ns = (dt - MJD_EPOCH_NUMPY).astype('timedelta64[ns]').astype(float)
+        return mjd_ns / 1e9 / 3600 / 24
+
+    # assume datetimes without timezone are utc
+    if isinstance(dt, datetime):
+        if dt.tzinfo is None:
+            dt = dt.replace(tzinfo=timezone.utc)
+
+    elif isinstance(dt, (pd.Series, pd.DatetimeIndex)):
+        if dt.tz is None:
+            dt.tz = timezone.utc
+
+    return (dt - MJD_EPOCH).total_seconds() / 24 / 3600
+
+
+def mjd_to_datetime(mjd):
+    if isinstance(mjd, (int, float)):
+        delta = timedelta(microseconds=mjd * 24 * 3600 * 1e6)
+        return MJD_EPOCH + delta
+
+    if isinstance(mjd, pd.Series):
+        delta = (mjd * 24 * 3600 * 1e9).astype('timedelta64[ns]')
+        return MJD_EPOCH + delta
+
+    # other types will be returned as numpy array if possible
+    mjd = np.asanyarray(mjd)
+    delta = (mjd * 24 * 3600 * 1e9).astype('timedelta64[ns]')
+    return MJD_EPOCH_NUMPY + delta
+
 
 def fjd(datetime_inst):
     ''' convert datetime instance to FJD

tests/test_time.py

@@ -0,0 +1,64 @@
+from datetime import datetime, timezone
+import pandas as pd
+import numpy as np
+
+
+def test_datetime_to_mjd():
+    from fact.time import datetime_to_mjd, MJD_EPOCH
+
+    assert datetime_to_mjd(MJD_EPOCH) == 0.0
+
+    dt = datetime(2017, 1, 1, 12, 0)
+    assert datetime_to_mjd(dt) == 57754.5
+
+
+def test_datetime_to_mjd_pandas():
+    from fact.time import datetime_to_mjd
+
+    dates = pd.date_range('2017-05-17 12:00', freq='1d', periods=5)
+
+    mjd = pd.Series([57890.5, 57891.5, 57892.5, 57893.5, 57894.5])
+
+    assert all(datetime_to_mjd(dates) == mjd)
+
+
+def test_datetime_to_mjd_numpy():
+    from fact.time import datetime_to_mjd
+
+    dates = np.array(['2017-05-17 12:00', '2017-05-18 12:00'], dtype='datetime64')
+    mjd = np.array([57890.5, 57891.5])
+
+    assert all(datetime_to_mjd(dates) == mjd)
+
+
+def test_mjd_to_datetime_float():
+    from fact.time import mjd_to_datetime, MJD_EPOCH
+
+    assert mjd_to_datetime(0) == MJD_EPOCH
+
+    dt = datetime(2017, 5, 17, 18, 0, tzinfo=timezone.utc)
+    assert mjd_to_datetime(57890.75) == dt
+
+
+def test_mjd_to_datetime_numpy():
+    from fact.time import mjd_to_datetime
+
+    mjd = np.arange(57890.0, 57891, 0.25)
+
+    dt = np.array(
+        ['2017-05-17 00:00', '2017-05-17 06:00',
+         '2017-05-17 12:00', '2017-05-17 18:00'],
+        dtype='datetime64[us]'
+    )
+
+    assert all(mjd_to_datetime(mjd) == dt)
+
+
+def test_mjd_to_datetime_pandas():
+    from fact.time import mjd_to_datetime
+
+    mjd = pd.Series(np.arange(57890.0, 57891, 0.25))
+
+    dt = pd.Series(pd.date_range('2017-05-17 00:00', freq='6h', periods=4))
+
+    assert all(mjd_to_datetime(mjd) == dt)