Add excerpt to the command line

This commit is contained in:
Brandon Rhodes 2018-02-11 13:39:02 -05:00
parent b5f23ddbe2
commit dd51993296
4 changed files with 143 additions and 62 deletions

View File

@ -37,16 +37,33 @@ JPL domain.
Command Line Tool
-----------------
If you have downloaded a ``.bsp`` file and want to learn what ephemeris
segments are stored inside of it, you can have ``jplephem`` print them
out by invoking the module directly from the command line::
If you have downloaded a ``.bsp`` file, you can run ``jplephem`` from
the command line to display the data inside of it::
python -m jplephem de430.bsp
python -m jplephem comment de430.bsp
python -m jplephem dap de430.bsp
python -m jplephem spk de430.bsp
This will print out a summary identical to the one shown in the
following section, but without requiring that you type and run any
Python code.
You can also take a large ephemeris and produce a smaller excerpt by
limiting the range of dates that it covers:
python -m jplephem excerpt 2018/1/1 2018/4/1 de421.bsp outjup.bsp
If the input ephemeris is a URL, then `jplephem` will try to save
bandwidth by fetching only the blocks of the remote file that are
necessary to cover the dates you have specified. For example, the
Jupiter satellite ephemeris `jup310.bsp` is famously large, weighing in
a nearly a gigabyte in length. But if all you need are Jupiter's
satellites for one month, you can download considerably less data:
$ python -m jplephem excerpt 2018/1/1 2018/4/1 \
https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/jup310.bsp \
excerpt.bsp
$ ls -lh excerpt.bsp
-rw-r----- 1 brandon brandon 1.2M Feb 11 13:36 excerpt.bsp
In this case only about one-thousandth of the ephemeris's data needed to
be downloaded, which took less than one minute.
Getting Started With DE430
--------------------------

View File

@ -5,6 +5,7 @@ from __future__ import print_function
import argparse
import sys
from .daf import DAF
from .excerpter import RemoteFile, write_excerpt
from .spk import SPK
def main(args):
@ -33,9 +34,10 @@ def main(args):
help="Create an SPK covering a narrower range of dates",
)
p.set_defaults(func=excerpt)
p.add_argument('start_date', help='Start date yyyy/mm/dd')
p.add_argument('end_date', help='End date yyyy/mm/dd')
p.add_argument('file', help='Local filename or remote URL')
p.add_argument('start_date', help='Start date yyyy/mm/dd', type=parse_date)
p.add_argument('end_date', help='End date yyyy/mm/dd', type=parse_date)
p.add_argument('path_or_url', help='Local filename or remote URL')
p.add_argument('output_path', help='Output file to create')
p = subparsers.add_parser(
'spk',
@ -68,8 +70,39 @@ def daf_segments(args):
' '.join(repr(v) for v in values))
def excerpt(args):
pass
if args.path_or_url.startswith(('http://', 'https://')):
url = args.path_or_url
f = RemoteFile(url)
spk = SPK(DAF(f))
with open(args.output_path, 'w+b') as output_file:
write_excerpt(spk, output_file, args.start_date, args.end_date)
else:
path = args.path_or_url
with open(path, 'rb') as f:
spk = SPK(DAF(f))
with open(args.output_path, 'w+b') as output_file:
write_excerpt(spk, output_file, args.start_date, args.end_date)
return ()
def spk_segments(args):
with open(args.path, 'rb') as f:
yield str(SPK(DAF(f)))
def parse_date(s):
try:
fields = [int(f) for f in s.split('/')]
except ValueError:
fields = []
if len(fields) < 1 or len(fields) > 3:
E = argparse.ArgumentTypeError
raise E('specify each date as YYYY or YYYY/MM or YYYY/MM/DD')
return julian_day(*fields)
def julian_day(year, month=1, day=1):
"""Given a proleptic Gregorian calendar date, return a Julian day int."""
janfeb = month < 3
return (day
+ 1461 * (year + 4800 - janfeb) // 4
+ 367 * (month - 2 + janfeb * 12) // 12
- 3 * ((year + 4900 - janfeb) // 100) // 4
- 32075)

80
jplephem/excerpter.py Normal file
View File

@ -0,0 +1,80 @@
"""Extract data for a specific date range from an SPK file."""
try:
from urllib.request import URLopener
except:
from urllib import URLopener
from numpy import copy
from .daf import DAF
from .spk import S_PER_DAY, T0
clip_lower = max
clip_upper = min
def _seconds(jd):
"""Convert a Julian Date to a number of seconds since J2000."""
return (jd - T0) * S_PER_DAY
def write_excerpt(input_spk, output_file, start_jd, end_jd):
start_seconds = _seconds(start_jd)
end_seconds = _seconds(end_jd)
old = input_spk.daf
# Copy the file record and the comments verbatim.
f = output_file
f.seek(0)
f.truncate()
for n in range(1, old.fward):
data = old.read_record(n)
f.write(data)
# Start an initial summary and name block.
summary_data = b'\0' * 1024
name_data = b' ' * 1024
f.write(summary_data)
f.write(name_data)
d = DAF(f)
d.fward = d.bward = old.fward
d.free = (d.fward + 1) * (1024 // 8) + 1
d.write_file_record()
# Copy over an excerpt of each array.
for name, values in old.summaries():
start, end = values[-2], values[-1]
init, intlen, rsize, n = old.read_array(end - 3, end)
rsize = int(rsize)
i = int(clip_lower(0, (start_seconds - init) // intlen))
j = int(clip_upper(n, (end_seconds - init) // intlen + 1))
init = init + i * intlen
n = j - i
extra = 4 # enough room to rebuild [init intlen rsize n]
excerpt = copy(old.read_array(
start + rsize * i,
start + rsize * j + extra - 1,
))
excerpt[-4:] = (init, intlen, rsize, n)
values = (init, init + n * intlen) + values[2:]
d.add_array(b'X' + name[1:], values, excerpt)
class RemoteFile(object):
def __init__(self, url):
self.opener = URLopener()
self.url = url
self.offset = 0
def seek(self, offset, whence=0):
assert whence == 0
self.offset = offset
def read(self, size):
start = self.offset
end = start + size - 1
assert end > start
h = 'Range', 'bytes={}-{}'.format(start, end)
self.opener.addheaders.append(h)
data = self.opener.open(self.url).read()
return data

View File

@ -3,23 +3,18 @@
http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/spk.html
"""
from numpy import array, copy, empty, empty_like, rollaxis
from numpy import array, empty, empty_like, rollaxis
from .daf import DAF
from .names import target_names
T0 = 2451545.0
S_PER_DAY = 86400.0
clip_lower = max
clip_upper = min
def jd(seconds):
"""Convert a number of seconds since J2000 to a Julian Date."""
return T0 + seconds / S_PER_DAY
def _seconds(jd):
"""Convert a Julian Date to a number of seconds since J2000."""
return (jd - T0) * S_PER_DAY
class SPK(object):
"""A JPL SPK ephemeris kernel for computing positions and velocities.
@ -70,50 +65,6 @@ class SPK(object):
"""Return the file comments, as a string."""
return self.daf.comments()
def excerpt(self, new_file, start_jd, end_jd):
start_seconds = _seconds(start_jd)
end_seconds = _seconds(end_jd)
old = self.daf
# Copy the file record and the comments verbatim.
f = new_file
f.seek(0)
f.truncate()
for n in range(1, old.fward):
data = old.read_record(n)
f.write(data)
# Start an initial summary and name block.
summary_data = b'\0' * 1024
name_data = b' ' * 1024
f.write(summary_data)
f.write(name_data)
d = DAF(f)
d.fward = d.bward = old.fward
d.free = (d.fward + 1) * (1024 // 8) + 1
d.write_file_record()
# Copy over an excerpt of each array.
for name, values in old.summaries():
start, end = values[-2], values[-1]
init, intlen, rsize, n = old.read_array(end - 3, end)
rsize = int(rsize)
i = int(clip_lower(0, (start_seconds - init) // intlen))
j = int(clip_upper(n, (end_seconds - init) // intlen + 1))
init = init + i * intlen
n = j - i
extra = 4 # enough room to rebuild [init intlen rsize n]
excerpt = copy(old.read_array(
start + rsize * i,
start + rsize * j + extra - 1,
))
excerpt[-4:] = (init, intlen, rsize, n)
values = (init, init + n * intlen) + values[2:]
d.add_array(b'X' + name[1:], values, excerpt)
class Segment(object):
"""A single segment of an SPK file.