Create a gist now

Instantly share code, notes, and snippets.

@kghose /horizons.py
Last active Mar 21, 2017

What would you like to do?
A bare minimum script to query the NASA HORIZONS system for positional data
"""
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')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment