Fix #493 by reworking GOCE reentry plot

brandon-rhodes · brandon-rhodes · commit 2d64f92623f7 · 2020-12-09T11:16:10.000-05:00
* Move the plot into its own Python file for faster iteration. * Switch the x-axis from ordinal dates to datetime objects, to dodge a breaking change in matplotlib: its default epoch for floating point days moved from the year AD 1 to the year AD 1970. * Upgrade Pandas to 1.0.0 to avoid: matplotlib/matplotlib#12310 pandas-dev/pandas#22859 * Switch to relying on `message` to determine if SGP4 failed. * Spruce up the axis labelling, just for fun. * Use the same plot in both places it’s discussed in the docs.
diff --git a/examples/goce_reentry_chart.py b/examples/goce_reentry_chart.py
@@ -0,0 +1,71 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.dates import HourLocator, DateFormatter
+
+from skyfield.api import load, EarthSatellite
+
+# Labels for both date and hour on the x axis, and km on y.
+
+def label_dates_and_hours(axes):
+    axes.xaxis.set_major_locator(HourLocator([0]))
+    axes.xaxis.set_minor_locator(HourLocator([0, 6, 12, 18]))
+    axes.xaxis.set_major_formatter(DateFormatter('0h\n%Y %b %d\n%A'))
+    axes.xaxis.set_minor_formatter(DateFormatter('%Hh'))
+    for label in ax.xaxis.get_ticklabels(which='both'):
+        label.set_horizontalalignment('left')
+    axes.yaxis.set_major_formatter('{x:.0f} km')
+    axes.tick_params(which='both', length=0)
+
+# Load the satellite's final TLE entry.
+
+sat = EarthSatellite(
+    '1 34602U 09013A   13314.96046236  .14220718  20669-5  50412-4 0   930',
+    '2 34602 096.5717 344.5256 0009826 296.2811 064.0942 16.58673376272979',
+    'GOCE',
+)
+
+# Build the time range `t` over which to plot, plus other values.
+
+ts = load.timescale()
+t = ts.tt_jd(np.arange(sat.epoch.tt - 2.0, sat.epoch.tt + 2.0, 0.005))
+reentry = ts.utc(2013, 11, 11, 0, 16)
+earth_radius_km = 6371.0
+
+# Compute geocentric positions for the satellite.
+
+g = sat.at(t)
+valid = [m is None for m in g.message]
+
+# Start a new figure.
+
+fig, ax = plt.subplots()
+
+# Draw the blue curve.
+
+x = t.utc_datetime()
+y = np.where(valid, g.distance().km - earth_radius_km, np.nan)
+ax.plot(x, y)
+
+# Label the TLE epoch.
+
+x = sat.epoch.utc_datetime()
+y = sat.at(sat.epoch).distance().km - earth_radius_km
+ax.plot(x, y, 'k.')
+ax.text(x, y - 9, 'Epoch of TLE data ', ha='right')
+
+# Label the official moment of reentry.
+
+x = reentry.utc_datetime()
+y = sat.at(reentry).distance().km - earth_radius_km
+ax.plot(x, y, 'r.')
+ax.text(x, y + 6, ' Moment of re-entry', c='r')
+
+# Grid lines and labels.
+
+ax.grid(which='both')
+ax.set(title='GOCE satellite: altitude above sea level', xlabel='UTC')
+label_dates_and_hours(ax)
+
+# Render the plot to a PNG file.
+
+fig.savefig('goce-reentry.png', bbox_inches='tight')
diff --git a/requirements.txt b/requirements.txt
@@ -5,7 +5,7 @@ lxml==4.4.1
 mock==2.0.0
 numpy==1.15.4
 matplotlib==3.3.0
-pandas==0.23.3
+pandas==1.0.0
 pyflakes==2.1.1
 python-dateutil>=2.5.0
 pytz
diff --git a/skyfield/documentation/_static/goce-decay.png b/skyfield/documentation/_static/goce-decay.png
diff --git a/skyfield/documentation/_static/goce-reentry.png b/skyfield/documentation/_static/goce-reentry.png
diff --git a/skyfield/documentation/earth-satellites.rst b/skyfield/documentation/earth-satellites.rst
@@ -687,18 +687,18 @@ At that point, the satellite will return a position and velocity
 are the special floating-point value ``nan`` which means *not-a-number*.
 
 When a propagation error occurs and you get ``nan`` values,
-you can examine the ``sgp4_error`` attribute of the returned position
+you can examine the ``message`` attribute of the returned position
 to learn the error that the SGP4 propagator encountered.
 
 We can take as an example the satellite elements above.
 They are the last elements ever issued for GOCE,
-about one day before the satellite re-entered the atmosphere
+just before the satellite re-entered the atmosphere
 after an extended and successful mission.
 Because of the steep decay of its orbit,
 the elements are valid over an unusually short period —
 from just before noon on Saturday to just past noon on Tuesday:
 
-.. image:: _static/goce-decay.png
+.. image:: _static/goce-reentry.png
 
 By asking for GOCE’s position just before or after this window,
 we can learn about the propagation errors
@@ -727,7 +727,7 @@ that are limiting this TLE set’s predictions:
     mrt is less than 1.0 which indicates the satellite has decayed
 
 If you use a ``Time`` array to ask about an entire range of dates,
-then ``sgp4_error`` will be a sequence filled in with ``None``
+then ``message`` will be a sequence filled in with ``None``
 whenever the SGP4 propagator was successful
 and otherwise recording the propagator error:
 
diff --git a/skyfield/documentation/example-plots.rst b/skyfield/documentation/example-plots.rst
@@ -75,62 +75,12 @@ The code to produce the diagram using `matplotlib`_,
 including custom tick marks that are based on the date,
 is:
 
-.. testcode::
-
-    from matplotlib import pyplot as plt
-    from matplotlib.dates import HourLocator, DateFormatter
-
-    from numpy import arange
-
-    def label_dates_and_hours(axes):
-        axes.xaxis.set_major_locator(HourLocator([0]))
-        axes.xaxis.set_minor_locator(HourLocator([0, 12]))
-        axes.xaxis.set_major_formatter(DateFormatter('\n%a %d'))
-        axes.xaxis.set_minor_formatter(DateFormatter('%Hh'))
-
-    from skyfield.api import load, EarthSatellite
-
-    # Load the satellite's final TLE entry.
-
-    sat = EarthSatellite(
-        '1 34602U 09013A   13314.96046236  .14220718  20669-5  50412-4 0   930',
-        '2 34602 096.5717 344.5256 0009826 296.2811 064.0942 16.58673376272979',
-        'GOCE',
-    )
-
-    # Build the time range `t` over which to plot, plus other values.
-
-    ts = load.timescale()
-    t = ts.tt_jd(arange(sat.epoch.tt - 1.0, sat.epoch.tt + 3.0, 0.01))
-    reentry = ts.utc(2013, 11, 11, 0, 16)
-    earth_radius_km = 6371.
-
-    # Start a new figure.
-
-    fig, ax = plt.subplots()
-
-    # Draw the blue curve.
-
-    x = t.toordinal()
-    y = sat.at(t).distance().km - earth_radius_km
-    ax.plot(x, y)
-
-    # Label the official moment of reentry.
-
-    x = reentry.toordinal()
-    y = sat.at(reentry).distance().km - earth_radius_km
-    ax.plot(x, y, 'ro')
-    ax.text(x, y + 10, 'Moment of re-entry')
-
-    # Grid lines and labels.
-
-    label_dates_and_hours(ax)
-    ax.grid()
-    ax.set(title='GOCE satellite altitude', ylabel='km above sea level')
+.. testsetup::
 
-    # Render the plot to a PNG file.
+    import goce_reentry_chart
 
-    fig.savefig('goce-reentry.png', bbox_inches='tight')
+.. include:: ../../examples/goce_reentry_chart.py
+   :literal:
 
 .. testcleanup::
 
