From d54977e7d2fb913393c82858b0b80088f7d7f57f Mon Sep 17 00:00:00 2001 From: Brandon Rhodes Date: Thu, 23 Jul 2020 23:40:51 -0400 Subject: [PATCH] Update test repr to match new ci/hip_main.dat.gz --- requirements.txt | 5 +- skyfield/documentation/example-plots.rst | 158 +++++++++++++++++++++++ skyfield/projections.py | 7 +- skyfield/tests/test_stars.py | 2 +- 4 files changed, 166 insertions(+), 6 deletions(-) diff --git a/requirements.txt b/requirements.txt index 2165489..8c51536 100755 --- a/requirements.txt +++ b/requirements.txt @@ -3,13 +3,12 @@ beautifulsoup4==4.6.0 html5lib==1.0.1 lxml==4.4.1 mock==2.0.0 -numpy==1.14.2 -matplotlib==2.2.2 +numpy==1.15.4 +matplotlib==3.3.0 pandas==0.23.3 pyflakes==2.1.1 python-dateutil>=2.5.0 pytz sphinx==1.7.2 https://github.com/brandon-rhodes/assay/archive/master.zip -https://github.com/whiskie14142/spktype21/archive/master.zip -e . diff --git a/skyfield/documentation/example-plots.rst b/skyfield/documentation/example-plots.rst index 8c76806..57fc0c4 100644 --- a/skyfield/documentation/example-plots.rst +++ b/skyfield/documentation/example-plots.rst @@ -9,12 +9,170 @@ for producing images from Skyfield computations. For the moment there’s only example so far, for plotting the elevation of a satellite over time: +Finder chart for comet NEOWISE +============================== + .. testsetup:: import matplotlib matplotlib.use('Agg') # to avoid “no display name” error on Travis CI del matplotlib +.. testcode:: + + import numpy as np + import pandas as pd + from matplotlib import pyplot as plt + + from skyfield import api + from skyfield.api import load + ts = load.timescale(builtin=True) + eph = load('de421.bsp') + sun = eph['sun'] + earth = eph['earth'] + + + # In[3]: + + + t_comet = ts.utc(2020, 7, range(17, 27)) + t = t_comet[len(t_comet) // 2] # middle date + + + # In[4]: + + + from skyfield.data import mpc + + with load.open(mpc.COMET_URL) as f: + comets = mpc.load_comets_dataframe(f) + + comets = comets.set_index('designation', drop=False) + row = comets.loc['C/2020 F3 (NEOWISE)'] + + + # In[5]: + + + from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN + + comet = sun + mpc.comet_orbit(row, ts, GM_SUN) + center = earth.at(t).observe(comet) + + + # In[6]: + + + from skyfield.api import Star + from skyfield.data import hipparcos + + + # In[53]: + + + from skyfield.projections import build_stereographic_projection + proj = build_stereographic_projection(center) + + + # In[54]: + + + with load.open(hipparcos.URL) as f: + stars = hipparcos.load_dataframe(f) + + + # In[55]: + + + star_positions = earth.at(t).observe(Star.from_dataframe(stars)) + stars['x'], stars['y'] = proj(star_positions) + + + # In[56]: + + + limiting_magnitude = 6.5 + bright_stars = (stars.magnitude <= limiting_magnitude) + + + # In[57]: + + + comet_x, comet_y = proj(earth.at(t_comet).observe(comet)) + + + # In[98]: + + + # = 'https://raw.githubusercontent.com/Stellarium/stellarium/master/skycultures/western_SnT/constellationship.fab' + from skyfield.data.stellarium import parse_constellations + + url = 'https://raw.githubusercontent.com/Stellarium/stellarium/master/skycultures/western_SnT/constellationship.fab' + with load.open(url) as f: + constellations = parse_constellations(f) + + edges_star1 = [star1 for name, edges in constellations for star1, star2 in edges] + edges_star2 = [star2 for name, edges in constellations for star1, star2 in edges] + + np.array([stars['x'].loc[edges_star1], stars['y'].loc[edges_star1]]) + xy1 = stars[['x', 'y']].loc[edges_star1].values + xy2 = stars[['x', 'y']].loc[edges_star2].values + lines_xy = np.rollaxis(np.array([xy1, xy2]), 1) + + + # In[94]: + + + fig, ax = plt.subplots(figsize=[9, 9]) + + from matplotlib.collections import LineCollection + from skyfield.api import tau + + field_of_view_degrees = 45.0 + + line_collection = LineCollection(lines_xy) #, color=['#0002'] * len(lines_xy)) + ax.add_collection(line_collection) + + marker_size = (0.5 + limiting_magnitude - stars.magnitude[bright_stars]) ** 2.0 + ax.scatter(stars['x'][bright_stars], stars['y'][bright_stars], s=marker_size, color='k') + #ax.scatter(comet_x, comet_y, s=100, color='b') + #comet_color = '#f00' + #ax.plot(comet_x, comet_y, '+', c=comet_color, zorder=3) + ax.plot(comet_x, comet_y, '+', zorder=3) + offset = 0.002 + for xi, yi, tstr in zip(comet_x, comet_y, t_comet.utc_strftime('%-m/%d')): + #text = ax.text(xi + offset, yi - offset, tstr, color=comet_color, ha='left', va='top', fontsize=12, + text = ax.text(xi + offset, yi - offset, tstr, ha='left', va='top', fontsize=12, + weight='bold') + text.set_alpha(0.3) + + ax.set_title('Comet NEOWISE {} through {}'.format( + t_comet[0].utc_strftime('%Y %B %d'), + t_comet[-1].utc_strftime('%Y %B %d'), + )) + ax.set_aspect(1.0) + ax.xaxis.set_visible(False) + ax.yaxis.set_visible(False) + + angle = (180.0 - field_of_view_degrees) / 720.0 * tau + limit = np.sin(angle) / (1.0 - np.cos(angle)) + ax.set_xlim(-limit, limit) + ax.set_ylim(-limit, limit) + + fig.savefig('neowise-finder-chart.png') + + + +.. image:: _static/neowise-finder-chart.png + +.. testcleanup:: + + import os + os.rename('neowise-finder-chart.png', '_static/neowise-finder-chart.png') + +Satellite altitude during re-entry +================================== + .. testcode:: from matplotlib import pyplot as plt diff --git a/skyfield/projections.py b/skyfield/projections.py index 8ac16bc..a652f54 100644 --- a/skyfield/projections.py +++ b/skyfield/projections.py @@ -43,8 +43,11 @@ def build_stereographic_projection(center): # https://math.stackexchange.com/questions/409217/ p = center.position.au u = p / length_of(p) - c = u.mean(axis=1) - c = c / length_of(c) + if len(u.shape) > 1: + c = u.mean(axis=1) + c = c / length_of(c) + else: + c = u x_c, y_c, z_c = c def project(position): diff --git a/skyfield/tests/test_stars.py b/skyfield/tests/test_stars.py index bff6a9e..a77e953 100644 --- a/skyfield/tests/test_stars.py +++ b/skyfield/tests/test_stars.py @@ -5,4 +5,4 @@ def test_dataframe(): with api.load.open('hip_main.dat.gz') as f: df = load_dataframe(f) star = api.Star.from_dataframe(df) - assert repr(star) == 'Star(ra shape=214, dec shape=214, ra_mas_per_year shape=214, dec_mas_per_year shape=214, parallax_mas shape=214, epoch shape=214)' + assert repr(star) == 'Star(ra shape=9933, dec shape=9933, ra_mas_per_year shape=9933, dec_mas_per_year shape=9933, parallax_mas shape=9933, epoch shape=9933)'