Add a command line script for ephemeris generation.

Author: mkelley

sbpy/data/ephem.py

@@ -1,15 +1,26 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """
 ======================
-sbpy data.Ephem Module
+sbpy.data.ephem Module
 ======================
 
-Class for storing and querying ephemerides
+Class for storing and querying ephemerides.
+
+Ephemerides may be generated by this module via command line:
+
+:code:
+    python3 -m sbpy.data.ephem
+
+Or with the ``sbpy-ephem`` command line script.
 
 created on June 04, 2017
+
 """
+
 import os
+import argparse
 from warnings import warn
+from typing import Union
 
 import numpy as np
 from numpy import ndarray, hstack, iterable
@@ -525,7 +536,7 @@ def from_miriade(cls, targetids, objtype='asteroid',
         * For detailed explanations of the queried fields, refer to
           `astroquery.imcce.MiriadeClass.get_ephemerides` and the
           `Miriade documentation
-          <http://vo.imcce.fr/webservices/miriade/?documentation>`_.
+          <https://vo.imcce.fr/webservices/miriade/?documentation>`_.
         * By default, all properties are provided in the J2000.0 reference
           system. Different settings can be chosen using
           additional keyword arguments as used by
@@ -863,3 +874,211 @@ def from_oo(cls, orbit, epochs=None, location='500', scope='full',
         ephem.table.remove_column('MJD')
 
         return ephem
+
+
+@requires("astroquery")
+def main():
+    """Command line script for ephemeris generation."""
+
+    from .names import Names, TargetNameParseError
+
+    today: Time = Time(Time.now().iso[:10])
+
+    parser: argparse.ArgumentParser = argparse.ArgumentParser()
+    parser.add_argument(
+        "service",
+        choices=("horizons", "miriade", "mpc"),
+        help="online service to use for ephemeris",
+    )
+    parser.add_argument(
+        "target",
+    )
+    parser.add_argument(
+        "--location",
+        "-l",
+        default="500",
+        help="observer location",
+    )
+    parser.add_argument(
+        "--start",
+        default=today,
+        type=Time,
+        help="start time (UTC)",
+    )
+    parser.add_argument(
+        "--stop", default=today + 30 * u.day, type=Time, help="stop time (UTC)"
+    )
+    step_mutex = parser.add_mutually_exclusive_group()
+    step_mutex.add_argument(
+        "--step",
+        default="1 d",
+        type=u.Quantity,
+        help="ephemeris step size with units or number of steps",
+    )
+    parser.add_argument(
+        "--radec",
+        "-r",
+        default="hmsdms",
+        choices=("hmsdms", "deg"),
+        help="RA, Dec coordinate format",
+    )
+
+    horizons: argparse._ArgumentGroup = parser.add_argument_group(
+        "JPL Horizons options"
+    )
+    horizons.add_argument(
+        "--id-type",
+        "-i",
+        default=None,
+        choices=(
+            "smallbody",
+            "designation",
+            "name",
+            "asteroid_name",
+            "comet_name",
+        ),
+        help=(
+            "target identifer type; automatically set to 'designation' for "
+            "target names that appear to be a cometary designation"
+        ),
+    )
+    horizons.add_argument(
+        "--quantities",
+        "-q",
+        default="1,3,9,19,20,23,24,27,33",
+        help="table quantities to return as a comma-separated list of integers",
+    )
+
+    miriade: argparse._ArgumentGroup = parser.add_argument_group(
+        "IMCCE Miriade options"
+    )
+    miriade.add_argument(
+        "--type",
+        "-t",
+        default="asteroid",
+        choices=(
+            "asteroid",
+            "comet",
+            "dwarf planet",
+            "planet",
+            "satellite",
+        ),
+        help="target identifer type",
+    )
+
+    args: argparse.Namespace = parser.parse_args()
+
+    eph: Ephem
+    epochs: dict = {
+        "start": args.start,
+        "stop": args.stop,
+    }
+    if args.service == "mpc":
+        if args.step.unit == u.dimensionless_unscaled:
+            epochs["number"] = int(args.step)
+        else:
+            epochs["step"] = args.step
+
+        eph = Ephem.from_mpc(
+            args.target,
+            epochs=epochs,
+            location=args.location,
+        )
+
+        # simplify output
+        for k in ["Targetname"]:
+            del eph.table[k]
+    elif args.service == "horizons":
+        epochs["step"] = (
+            int(args.step) if args.step.unit == u.dimensionless_unscaled else args.step
+        )
+
+        # comet specific options: always avoid multiple fragment matches, always
+        # return closest apparition
+        closest_apparition: bool = False
+        no_fragments: bool = False
+        id_type: Union[str, None] = args.id_type
+        try:
+            Names.parse_comet(args.target)
+            closest_apparition = True
+            no_fragments = True
+            id_type = "designation" if id_type is None else id_type
+        except TargetNameParseError:
+            pass
+
+        eph = Ephem.from_horizons(
+            args.target,
+            id_type=id_type,
+            closest_apparition=closest_apparition,
+            no_fragments=no_fragments,
+            epochs=epochs,
+            location=args.location,
+            quantities=args.quantities,
+        )
+
+        # simplify output
+        eph["epoch"].format = "iso"
+        eph["ra"] = eph["ra"].unmasked
+        eph["dec"] = eph["dec"].unmasked
+        for k in ["M1", "k1", "solar_presence", "lunar_presence", "targetname"]:
+            del eph.table[k]
+    elif args.service == "miriade":
+        if args.step.unit == u.dimensionless_unscaled:
+            epochs["number"] = int(args.step)
+        else:
+            epochs["step"] = args.step
+
+        eph = Ephem.from_miriade(
+            args.target,
+            objtype=args.type,
+            epochs=epochs,
+            location=args.location,
+            coordtype=1,
+        )
+
+        # simplify output
+        eph["epoch"].format = "iso"
+        for k in ["target"]:
+            del eph.table[k]
+        for k in eph.field_names:
+            if hasattr(eph[k], "unmasked"):
+                eph[k] = eph[k].unmasked
+
+    # convert RA and Dec to Angle in units of degree
+    eph["ra"] = Angle(eph["ra"], eph["ra"].unit).to("deg")
+    eph["dec"] = Angle(eph["dec"], eph["dec"].unit).to("deg")
+
+    if args.radec == "hmsdms":
+        eph["ra"].info.format = lambda x: x.to_string(
+            sep=":", precision=2, unit="hourangle"
+        )
+        eph["dec"].info.format = lambda x: x.to_string(sep=":", precision=1)
+    else:
+        eph["ra"].info.format = lambda x: x.to_string()
+        eph["dec"].info.format = lambda x: x.to_string()
+
+    # unified output order for most common columns
+    fields = eph.field_names
+    for k in eph._translate_columns(
+        [
+            "RA*cos(Dec)_rate",
+            "DEC_rate",
+            "elong",
+            "phase",
+            "delta",
+            "rh",
+            "dec",
+            "ra",
+            "epoch",
+        ],
+        ignore_missing=True,
+    ):
+        if k not in fields:
+            continue
+        fields.insert(0, fields.pop(fields.index(k)))
+
+    eph[fields].table.pprint_all()
+
+
+if __name__ == "__main__":
+    main()

setup.cfg

@@ -67,6 +67,9 @@ docs =
 ginga.rv.plugins =
     cometaryenhancements = sbpy.ginga_plugins:setup_cometaryenhancements
 
+console_scripts =
+    sbpy-ephem = sbpy.data.ephem:main
+
 [tool:pytest]
 minversion = 7.0
 testpaths = "sbpy" "docs"