A bare minimum script to query the NASA HORIZONS system for positional data

horizons.py

"""
The URL format is documented at http://ssd.jpl.nasa.gov/horizons_batch.cgi
I figured out the template by using the webinterface to make a query I wanted and then
clicking 'batch-file' to see how the batch file commands would look like.

Also see https://github.com/mommermi/callhorizons by Michael Mommert

"""
import urllib.request
import time

from matplotlib import pyplot as plt
import numpy as np


def query(target, start, stop, step, center="500@0", retry_count=50):
  query.url = '&'.join([
    "https://ssd.jpl.nasa.gov/horizons_batch.cgi?batch=1",
    "COMMAND="+target,
    "CENTER="+center,
    "MAKE_EPHEM='YES'",
    "TABLE_TYPE='VECTORS'",
    "START_TIME="+start,
    "STOP_TIME="+stop,
    "STEP_SIZE="+step,
    "OUT_UNITS='KM-S'",
    "REF_PLANE='ECLIPTIC'",
    "REF_SYSTEM='J2000'",
    "VECT_CORR='NONE'",
    "VEC_LABELS='NO'",
    "VEC_DELTA_T='NO'",
    "CSV_FORMAT='YES'",
    "OBJ_DATA='YES'",
    "VEC_TABLE='1'"
  ])

  for n in range(retry_count):
    try:
      return parse(urllib.request.urlopen(query.url).readlines())
    except urllib.request.URLError:
      time.sleep(0.1)
      # in case the HORIZONS website is blocked (due to another query)
      # wait 0.1 second and try again
  else:
    raise TimeoutError


def parse(lines):
  return {
    'raw_text': lines,
    'ephemerides': parse_ephemerides(lines)
  }


def parse_ephemerides(lines):
  se = [n for n, l in enumerate(lines) if l.startswith(b'$$SOE')] + \
       [n for n, l in enumerate(lines) if l.startswith(b'$$EOE')]
  return np.array(
    [parse_eph_line(l) for l in lines[se[0] + 1:se[1]]],
    dtype=[('jdtdb', float), ('x', float), ('y', float), ('z', float)]
  )


def parse_eph_line(line):
  sp = line.split(b',')
  return (sp[0],) + tuple(sp[2:5])


def plot_query(target, start, stop, step, center="0", retry_count=50):
  res = query(target=target, center=center, start=start, stop=stop, step=step, retry_count=retry_count)

  eph = res['ephemerides']

  fig = plt.figure(figsize=(4, 4))
  plt.plot(eph['x'], eph['y'])
  plt.gca().set_aspect('equal')

  fig = plt.figure(figsize=(4, 4))
  plt.plot(eph['x'], eph['z'])


# https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/naif_ids.html
oi = object_ids = {
  'solar-bary': '0',
  'em-bary': '3',
  'earth': '399',
  'moon': '301',
  'mercury': '199',
  'jupiter': '599',
  'saturn': '699'
}

plot_query(target=oi['earth'], center=oi['solar-bary'],
           start='2100-01-01', stop='2100-12-31', step='1d')