Add support for Barycentric Julian Dates

.gitmodules

@@ -0,0 +1,4 @@
+[submodule "barycorr"]
+	path = barycorr
+	url = https://github.com/vterron/barycorr.git
+	branch = init

Doc/user/install.rst

@@ -13,7 +13,7 @@ These are the steps to install LEMON on a fresh `Debian 7`_ machine:
 #. ``apt-get build-dep python-matplotlib python-scipy``
 #. ``apt-get install openmpi-dev``
 #. ``easy_install -U distribute``
-#. ``git clone --branch v0.3 git://github.com/vterron/lemon.git ~/lemon``
+#. ``git clone --recurse-submodules git://github.com/vterron/lemon.git ~/lemon``
 #. ``cd ~/lemon``
 #. ``pip install "numpy>=1.7.1"``
 #. ``pip install -r pre-requirements.txt`` # :download:`[View] <../../pre-requirements.txt>`

README.rst

@@ -21,6 +21,13 @@ Homepage
 
 http://lemon.readthedocs.org/
 
+Acknowledgments
+---------------
+
+LEMON uses `Eastman et al. (2010), PASP 122, pp. 935–946`__ for  Barycentric Julian Date calculations.
+
+__ http://adsabs.harvard.edu/abs/2010PASP..122..935E
+
 Acknowledging us
 ----------------
 

fitsimage.py

@@ -28,8 +28,10 @@
 # tasks which attempt to display graphics will fail, of course, but we are not
 # going to make use of any of them, anyway.
 
-import astropy.io
+
+import astropy.time
 import astropy.wcs
+import barycorr
 import calendar
 import collections
 import datetime
@@ -187,7 +189,9 @@ def read_keyword(self, keyword):
         if not keyword:
             raise ValueError("keyword cannot be empty")
         try:
-            return self._header[keyword.upper()]
+            value = self._header[keyword.upper()]
+            logging.debug("%s: %s = %s", self.path, keyword, value)
+            return value
         except KeyError:
             msg = "%s: keyword '%s' not found" % (self.path, keyword)
             raise KeyError(msg)
@@ -326,15 +330,66 @@ def add_history(self, history):
         finally:
             handler.close(output_verify = 'ignore')
 
-    def date(self, date_keyword = 'DATE-OBS', time_keyword = 'TIME-OBS',
-             exp_keyword = 'EXPTIME'):
+    def date(self,
+             date_keyword=keywords.datek,
+             time_keyword=keywords.timek,
+             exp_keyword=keywords.exptimek,
+             barycentric=False,
+             bjd_keyword=keywords.bjdk,
+             ra_keyword=keywords.rak,
+             dec_keyword=keywords.deck,
+        ):
+        """ Return the date at mid-exposure (UTC), in Unix time.
+
+        This function always returns a UTC Unix timestamp. By default, the date
+        at mid-exposure is derived using the traditional, standard FITS keywords
+        (e.g. DATE-OBS, EXPTIME). However, if 'barycentric' is set to True, the
+        value used is the Barycentric Julian Date in Barycentric Dynamical Time
+        (BJD_TDB)). This is a non-standard FITS keyword, but commonly used in
+        the exoplanet community for high-accuracy observations. For more
+        information, see http://astroutils.astronomy.ohio-state.edu/time/.
+        """
+
+        if barycentric:
+            return self._read_barycentric_date(
+                bjd_keyword, ra_keyword, dec_keyword)
+        return self._read_uncorrected_date(
+            date_keyword, time_keyword, exp_keyword=exp_keyword)
+
+    def _read_barycentric_date(self, bjd_keyword, ra_keyword, dec_keyword):
+        """Returns the Barycentric Julian Date as a UTC timestamp."""
+
+        bjd_tdb = self.read_keyword(bjd_keyword)
+        ra = self.ra(ra_keyword=ra_keyword)
+        dec = self.dec(dec_keyword=dec_keyword)
+
+        logging.debug(
+            "Converting %f (%s, α=%f, δ=%f) to JD UTC",
+            bjd_tdb, bjd_keyword, ra, dec)
+
+        jd_utc = barycorr.bjd2utc(bjd_tdb, ra, dec)
+
+        logging.debug(
+            "%f (%s, α=%f, δ=%f) = %f (JD UTC)",
+            bjd_tdb, bjd_keyword, ra, dec, jd_utc)
+
+        t = astropy.time.Time(jd_utc, format='jd', scale='utc')
+        logging.debug("%f (JD UTC) = %f (UTC timestamp)", jd_utc, t.unix)
+        return t.unix
+
+    def _read_uncorrected_date(self, date_keyword, time_keyword, exp_keyword):
         """ Return the date of observation (UTC), in Unix time.
 
         This method returns, in seconds after the Unix epoch (the time 00:00:00
         UTC on 1 January 1970), the date at the 'midpoint' of the observation.
         This is defined as the time of the start of the exposure plus one-half
-        the exposure duration. It must be pointed out that these dates are
-        *always* interpreted as Coordinated Universal Time (UTC).
+        the exposure duration.
+
+        This function makes no assumptions or corrections for the time reference
+        frame (e.g. geocentric, heliocentric or barycentric; i.e. the different
+        geometric locations from which one could measure time, differing by the
+        light-travel time between them). It simply just returns the date found
+        in the header as a UTC timestamp.
 
         The KeyError exception is raised if any of the specified keywords
         cannot be found in the FITS header. NonStandardFITS is thrown if the
@@ -541,7 +596,7 @@ def ra(self, ra_keyword = keywords.rak):
             logging.debug(msg3)
 
             # HH:MM:SS[.sss]
-            regexp = '^(?P<hh>\d{2}):(?P<mm>\d{2}):(?P<ss>\d{2}(\.\d{0,3})?)$'
+            regexp = '^(?P<hh>\d{2})[:\s](?P<mm>\d{2})[:\s](?P<ss>\d{2}(\.\d{0,3})?)$'
             match = re.match(regexp, ra_str.strip())
             if match:
                 hh =   int(match.group('hh'))
@@ -556,8 +611,8 @@ def ra(self, ra_keyword = keywords.rak):
                 logging.debug("{0}: α = {1:.5f}".format(self.path, ra))
                 return ra
             else:
-                msg = "{0}: '{1}' not in decimal degrees or 'HH:MM:SS.sss' format"
-                raise ValueError(msg.format(self.path, ra_keyword))
+                msg = "{0}: '{1}' not in decimal degrees or 'HH:MM:SS.sss' format (got = {2!r})"
+                raise ValueError(msg.format(self.path, ra_keyword, ra_str))
 
     def dec(self, dec_keyword = keywords.deck):
         """ Return the declination, in decimal degrees.
@@ -590,7 +645,7 @@ def dec(self, dec_keyword = keywords.deck):
             logging.debug(msg3)
 
             # DD:MM:SS[.sss]
-            regexp = '^(?P<dd>([-+])?\d{2}):(?P<mm>\d{2}):(?P<ss>\d{2}(\.\d{0,3})?)$'
+            regexp = '^(?P<dd>([-+])?\d{2})[:\s](?P<mm>\d{2})[:\s](?P<ss>\d{2}(\.\d{0,3})?)$'
             match = re.match(regexp, dec_str.strip())
             if match:
                 dd =   int(match.group('dd'))
@@ -605,8 +660,8 @@ def dec(self, dec_keyword = keywords.deck):
                 logging.debug("{0}: δ = {1:.5f}".format(self.path, dec))
                 return dec
             else:
-                msg = "{0}: '{1}' not in decimal degrees or 'DD:MM:SS.sss' format"
-                raise ValueError(msg.format(self.path, dec_keyword))
+                msg = "{0}: '{1}' not in decimal degrees or 'DD:MM:SS.sss' format (got = {2!r}"
+                raise ValueError(msg.format(self.path, dec_keyword, dec_str))
 
     @property
     def prefix(self):

keywords.py

@@ -104,3 +104,11 @@
 "determine the number of counts at which saturation arises in coadded " \
 "observations. If the keyword is missing, we assume a value of one (that " \
 "is, that the observation consisted of a single exposure) [default: %default]"
+
+bjdk = 'BJD_TDB'
+desc['bjdk'] = (
+"keyword for the barycentric JD (TDB, Barycentric Dynamical Time) at "
+"mid-exposure. Internally, LEMON stores all times as UTC Unix timestamps, so "
+"BJD_TDB is converted to JD_UTC using Jason Eastman's web applet for time "
+"correction (htttp://astroutils.astronomy.ohio-state.edu/time/bjd2utc.html) "
+"[default: %default]")

photometry.py

@@ -190,9 +190,15 @@ def parallel_photometry(args):
     pfilter = image.pfilter(options.filterk)
     logging.debug("%s: filter = %s" % (image.path, pfilter))
 
-    kwargs = dict(date_keyword = options.datek,
-                  time_keyword = options.timek,
-                  exp_keyword = options.exptimek)
+    kwargs = dict(
+        date_keyword=options.datek,
+        time_keyword=options.timek,
+        exp_keyword=options.exptimek,
+        barycentric=options.barycentric,
+        bjd_keyword=options.bjdk,
+        ra_keyword=options.rak,
+        dec_keyword=options.deck,
+    )
     unix_time = image.date(**kwargs)
     msg = "%s: observation date: %.2f (%s)"
     args = (image.path, unix_time, util.utctime(unix_time))
@@ -460,6 +466,36 @@ def parallel_photometry(args):
                      dest = 'uncimgk', default = keywords.uncimgk,
                      help = keywords.desc['uncimgk'])
 
+key_group.add_option('--rak', action = 'store', type = 'str',
+                     dest = 'rak', default = keywords.rak,
+                     help = keywords.desc['rak'])
+
+key_group.add_option('--deck', action = 'store', type = 'str',
+                     dest = 'deck', default = keywords.deck,
+                     help = keywords.desc['deck'])
+
+bjd_group = optparse.OptionGroup(parser,
+    "Barycentric Julian Dates [experimental]" ,
+    "By default, LEMON determines the date of mid-exposure via the --datek, "
+    "--timek and --expk keywords in the FITS header. However, high-accuracy "
+    "observations (e.g. transits of binaries or exoplanets) commonly use "
+    "Barycentric Julian Dates in Barycentric Dynamical Time (BJD_TDB). This "
+    "corrects for the amount of time it took light to travel to the changing "
+    "location of the observatory. For a brief, non-technical explanation, see "
+    "http://astroutils.astronomy.ohio-state.edu/time/bjd_explanation.html")
+
+bjd_group.add_option('--use_barycentric', action = 'store_true',
+                     dest = 'barycentric',
+                     help = "Use Barycentric Julian Dates in Barycentric "
+                     "Dynamical Time (BJD_TDB). These dates are read from "
+                     "the --bjdk FITS keyword.")
+
+bjd_group.add_option('--bjdk', action = 'store', type = 'str',
+                     dest = 'bjdk', default = keywords.bjdk,
+                     help = keywords.desc['bjdk'])
+
+parser.add_option_group(bjd_group)
+
 parser.add_option_group(key_group)
 customparser.clear_metavars(parser)
 
@@ -675,9 +711,15 @@ def main(arguments = None):
 
     def get_date(img):
         """ Return the date() of a FITSImage object """
-        return img.date(date_keyword = options.datek,
-                        time_keyword = options.timek,
-                        exp_keyword = options.exptimek)
+        return img.date(
+            date_keyword=options.datek,
+            time_keyword=options.timek,
+            exp_keyword=options.exptimek,
+            barycentric=options.barycentric,
+            bjd_keyword=options.bjdk,
+            ra_keyword=options.rak,
+            dec_keyword=options.deck,
+        )
 
     files = fitsimage.InputFITSFiles()
     img_dates = {}

requirements.txt

@@ -6,6 +6,7 @@
 # setup.py files, not pip's fault. Serial installation is the only
 # answer here [https://github.com/pypa/pip/issues/25]
 
+absl-py>=0.9.0
 APLpy>=0.9.9
 scipy>=0.12.0
 matplotlib>=1.2.1
@@ -17,3 +18,5 @@ unittest2>=0.5.1
 montage-wrapper>=0.9.7
 requests>=2.0.1
 subprocess32>=3.2.6
+requests-cache>=0.5.2
+requests-mock>=1.7.0