entrouvert
/
debian-skyfield
17
0
Fork 0
Code Issues Pull Requests Releases Activity

Give EarthSatellite builtin timescale by default 

Earth satellites used to use a timescale without leap seconds to
represent the epoch, which caused problems when users tried using that
date in calculations, so a `ts` parameter was added to the constructor.
But now that Skyfield has the concept of builtin leap second files,
let’s switch Earth satellites to those if no `ts` is passed.
This commit is contained in:
Brandon Rhodes 2020-06-12 11:21:31 -04:00
parent 0c86b6d568
commit 14356ac986
4 changed files with 138 additions and 137 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
skyfield/documentation/_static/goce-reentry.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 42 KiB

After

Width:  |  Height:  |  Size: 42 KiB

124
skyfield/io_timescale.py Normal file
View File

@ -0,0 +1,124 @@
import locale
import numpy as np
from datetime import date, datetime, timedelta
from pkgutil import get_data
from threading import Lock
from .timelib import Timescale, julian_date
try:
from io import BytesIO
except:
from StringIO import StringIO as BytesIO
_lock = Lock()
def _build_builtin_timescale():
b = get_data('skyfield', 'data/deltat.data')
expiration_date, data = parse_deltat_data(BytesIO(b))
b = get_data('skyfield', 'data/deltat.preds')
expiration_date, preds = parse_deltat_preds(BytesIO(b))
data_end_time = data[0, -1]
i = np.searchsorted(preds[0], data_end_time, side='right')
delta_t_recent = np.concatenate([data, preds[:,i:]], axis=1)
b = get_data('skyfield', 'data/Leap_Second.dat')
expiration_date, arrays = parse_leap_seconds(BytesIO(b))
leap_dates, leap_offsets = arrays
return Timescale(delta_t_recent, leap_dates, leap_offsets)
def parse_deltat_data(fileobj):
"""Parse the United States Naval Observatory ``deltat.data`` file.
Each line file gives the date and the value of Delta T::
2016 2 1 68.1577
This function returns a 2xN array of raw Julian dates and matching
Delta T values.
"""
array = np.loadtxt(fileobj)
year, month, day = array[-1,:3].astype(int)
expiration_date = date(year + 1, month, day)
year, month, day, delta_t = array.T
data = np.array((julian_date(year, month, day), delta_t))
return expiration_date, data
def parse_deltat_preds(fileobj):
"""Parse the United States Naval Observatory ``deltat.preds`` file.
The old format supplies a floating point year, the value of Delta T,
and one or two other fields::
2015.75 67.97 0.210 0.02
The new format adds a modified Julian day as the first field:
58484.000 2019.00 69.34 -0.152 0.117
This function returns a 2xN array of raw Julian dates and matching
Delta T values.
"""
lines = iter(fileobj)
header = next(lines)
if header.startswith(b'YEAR'):
# Format in use until 2019 February
next(lines) # discard blank line
year_float, delta_t = np.loadtxt(lines, usecols=[0, 1]).T
else:
# Format in use since 2019 February
year_float, delta_t = np.loadtxt(lines, usecols=[1, 2]).T
year = year_float.astype(int)
month = 1 + (year_float * 12.0).astype(int) % 12
expiration_date = date(year[0] + 2, month[0], 1)
data = np.array((julian_date(year, month, 1), delta_t))
return expiration_date, data
def parse_leap_seconds(fileobj):
"""Parse the IERS file ``Leap_Second.dat``.
The leap dates array can be searched with::
index = np.searchsorted(leap_dates, jd, 'right')
The resulting index allows (TAI - UTC) to be fetched with::
offset = leap_offsets[index]
"""
lines = iter(fileobj)
for line in lines:
if line.startswith(b'# File expires on'):
break
else:
raise ValueError('Leap_Second.dat is missing its expiration date')
line = line.decode('ascii')
with _lock:
original_locale = locale.setlocale(locale.LC_ALL)
locale.setlocale(locale.LC_ALL, 'C')
try:
dt = datetime.strptime(line, '# File expires on %d %B %Y\n')
finally:
locale.setlocale(locale.LC_ALL, original_locale)
# The file went out of date at the beginning of July 2016, and kept
# downloading every time a user ran a Skyfield program. So we now
# build in a grace period:
grace_period = timedelta(days=30)
expiration_date = dt.date() + grace_period
mjd, day, month, year, offsets = np.loadtxt(lines).T
leap_dates = np.ndarray(len(mjd) + 2)
leap_dates[0] = '-inf'
leap_dates[1:-1] = mjd + 2400000.5
leap_dates[-1] = 'inf'
leap_offsets = np.ndarray(len(mjd) + 2)
leap_offsets[0] = leap_offsets[1] = offsets[0]
leap_offsets[2:] = offsets
return expiration_date, (leap_dates, leap_offsets)

126
skyfield/iokit.py
View File

@ -3,28 +3,21 @@ from __future__ import print_function
import itertools
import os
import errno
import locale
import numpy as np
import sys
from datetime import date, datetime, timedelta
from datetime import date
from fnmatch import fnmatch
from pkgutil import get_data
from threading import Lock
from time import time
import certifi
from .io_timescale import (_build_builtin_timescale, parse_deltat_data,
parse_deltat_preds, parse_leap_seconds)
from .jpllib import SpiceKernel
from .sgp4lib import EarthSatellite
from .timelib import Timescale, julian_date
try:
from io import BytesIO
except:
from StringIO import StringIO as BytesIO
from .timelib import Timescale
today = date.today
_lock = Lock()
try:
from fcntl import LOCK_EX, LOCK_UN, lockf
@ -322,19 +315,7 @@ class Loader(object):
"""
if builtin:
b = get_data('skyfield', 'data/deltat.data')
expiration_date, data = parse_deltat_data(BytesIO(b))
b = get_data('skyfield', 'data/deltat.preds')
expiration_date, preds = parse_deltat_preds(BytesIO(b))
data_end_time = data[0, -1]
i = np.searchsorted(preds[0], data_end_time, side='right')
delta_t_recent = np.concatenate([data, preds[:,i:]], axis=1)
b = get_data('skyfield', 'data/Leap_Second.dat')
expiration_date, arrays = parse_leap_seconds(BytesIO(b))
leap_dates, leap_offsets = arrays
return Timescale(delta_t_recent, leap_dates, leap_offsets)
return _build_builtin_timescale()
if delta_t is not None:
delta_t_recent = np.array(((-1e99, 1e99), (delta_t, delta_t)))
@ -396,103 +377,6 @@ def load_file(path):
return SpiceKernel(path)
raise ValueError('unrecognized file extension: {}'.format(path))
def parse_deltat_data(fileobj):
"""Parse the United States Naval Observatory ``deltat.data`` file.
Each line file gives the date and the value of Delta T::
2016 2 1 68.1577
This function returns a 2xN array of raw Julian dates and matching
Delta T values.
"""
array = np.loadtxt(fileobj)
year, month, day = array[-1,:3].astype(int)
expiration_date = date(year + 1, month, day)
year, month, day, delta_t = array.T
data = np.array((julian_date(year, month, day), delta_t))
return expiration_date, data
def parse_deltat_preds(fileobj):
"""Parse the United States Naval Observatory ``deltat.preds`` file.
The old format supplies a floating point year, the value of Delta T,
and one or two other fields::
2015.75 67.97 0.210 0.02
The new format adds a modified Julian day as the first field:
58484.000 2019.00 69.34 -0.152 0.117
This function returns a 2xN array of raw Julian dates and matching
Delta T values.
"""
lines = iter(fileobj)
header = next(lines)
if header.startswith(b'YEAR'):
# Format in use until 2019 February
next(lines) # discard blank line
year_float, delta_t = np.loadtxt(lines, usecols=[0, 1]).T
else:
# Format in use since 2019 February
year_float, delta_t = np.loadtxt(lines, usecols=[1, 2]).T
year = year_float.astype(int)
month = 1 + (year_float * 12.0).astype(int) % 12
expiration_date = date(year[0] + 2, month[0], 1)
data = np.array((julian_date(year, month, 1), delta_t))
return expiration_date, data
def parse_leap_seconds(fileobj):
"""Parse the IERS file ``Leap_Second.dat``.
The leap dates array can be searched with::
index = np.searchsorted(leap_dates, jd, 'right')
The resulting index allows (TAI - UTC) to be fetched with::
offset = leap_offsets[index]
"""
lines = iter(fileobj)
for line in lines:
if line.startswith(b'# File expires on'):
break
else:
raise ValueError('Leap_Second.dat is missing its expiration date')
line = line.decode('ascii')
with _lock: # won't help if anyone user threads are doing parsing, alas
original_locale = locale.setlocale(locale.LC_ALL)
locale.setlocale(locale.LC_ALL, 'C')
try:
dt = datetime.strptime(line, '# File expires on %d %B %Y\n')
finally:
locale.setlocale(locale.LC_ALL, original_locale)
# The file went out of date at the beginning of July 2016, and kept
# downloading every time a user ran a Skyfield program. So we now
# build in a grace period:
grace_period = timedelta(days=30)
expiration_date = dt.date() + grace_period
mjd, day, month, year, offsets = np.loadtxt(lines).T
leap_dates = np.ndarray(len(mjd) + 2)
leap_dates[0] = '-inf'
leap_dates[1:-1] = mjd + 2400000.5
leap_dates[-1] = 'inf'
leap_offsets = np.ndarray(len(mjd) + 2)
leap_offsets[0] = leap_offsets[1] = offsets[0]
leap_offsets[2:] = offsets
return expiration_date, (leap_dates, leap_offsets)
def parse_tle(fileobj):
"""Parse a file of TLE satellite element sets.

25
skyfield/sgp4lib.py
View File

@ -8,24 +8,13 @@ from sgp4.api import SGP4_ERRORS, Satrec
from .constants import AU_KM, DAY_S, T0, tau
from .functions import _mxv, rot_x, rot_y, rot_z
from .io_timescale import _build_builtin_timescale
from .positionlib import ITRF_to_GCRS2
from .searchlib import _find_discrete, find_maxima
from .timelib import Timescale, calendar_date
from .vectorlib import VectorFunction
_minutes_per_day = 1440.
# Since satellite calculations are done entirely in UTC, we can display
# and return their epoch using a null Timescale. Let's just hope that
# no one ever pulls the .epoch off of an EarthSatellite and tries to
# generate high precision positions with it. If that becomes a problem,
# we can make a more severe change to the satellite API to require a
# Timescale when the user asks for the epoch as a Time. (It was
# probably a mistake to have an .epoch convenience attribute.)
_identity = identity(3)
_infs = array(('-inf', 'inf'), float)
_ts = Timescale(array((_infs, (0.0, 0.0))), _infs, array((37.0, 37.0)))
class EarthSatellite(VectorFunction):
"""An Earth satellite loaded from a TLE file and propagated with SGP4.
@ -84,9 +73,13 @@ class EarthSatellite(VectorFunction):
"""
center = 399
ts = None # see __init__()
def __init__(self, line1, line2, name=None, ts=None):
ts = ts or _ts
if ts is None:
ts = self.ts
if ts is None:
ts = EarthSatellite.ts = _build_builtin_timescale()
self.name = None if name is None else name.strip()
satrec = Satrec.twoline2rv(line1, line2)
@ -100,9 +93,9 @@ class EarthSatellite(VectorFunction):
self.epoch = ts.utc(year, 1, satrec.epochdays)
self._setup(satrec, ts)
self._setup(satrec)
def _setup(self, satrec, ts):
def _setup(self, satrec):
# If only I had not made __init__() specific to TLE lines, but
# had put them in an alternate construtor instead, this would
# simply have lived in __init__(). Alas! I was so young then.
@ -127,7 +120,7 @@ class EarthSatellite(VectorFunction):
year, month, day = calendar_date(satrec.jdsatepoch)
self.epoch = ts.utc(year, month, day + satrec.jdsatepochF)
self._setup(satrec, ts)
self._setup(satrec)
return self
def __str__(self):