take_screenshots = True
%pylab wx
if take_screenshots:
%pylab inline
rcParams['figure.figsize'] = (8.0, 8.0)
Populating the interactive namespace from numpy and matplotlib Populating the interactive namespace from numpy and matplotlib
from mayavi import mlab
ocean_blue = (0.4, 0.5, 1.0)
yellow = (1.0, 1.0, 0.0)
start_date = (2013, 1, 11, 15, 45, 0)
We use the sgp4
Python library
to interpret CelesTrak two-line orbital element entries from a text file.
from skyfield.api import JulianDate, earth, now
def find_satellites(lines, earth=earth):
"""Yield up satellites found in a sequence of text lines.
Each item yielded is a tuple ``('name', satellite)`` or
``(None, satellite)`` if the satellite TLE entry does not
appear to have been preceded by a name on the line above.
"""
line0 = line1 = None
for line2 in lines:
line2 = line2.rstrip('\r\n')
if (line1 and line1.startswith('1 ')
and line2.startswith('2 ')
and len(line1) == len(line2) == 69):
if line0 and len(line0) == 24:
name = line0.strip()
else:
name = None
sat = earth.satellite(line1 + '\n' + line2)
yield name, sat
line0 = line1
line1 = line2
with open('data/visual.txt') as lines:
names, satellites = zip(*find_satellites(lines))
print 'Successfully loaded', len(satellites), 'satellites'
Successfully loaded 145 satellites
from sgp4.earth_gravity import wgs72
r = wgs72.radiusearthkm
class Geocenter(object):
"""Coordinates of the geocenter."""
rGCRS = array([0, 0, 0])
vGCRS = array([0, 0, 0])
geocenter = Geocenter()
geocenter.jd = now()
def positions_at(jd):
geocenter.jd = jd
positions = []
for sat in satellites:
try:
position = sat.observe_from(geocenter).position.km
except TypeError:
pass
else:
positions.append(position)
return array(positions).T
tt0 = now().tt
minute = 1.0 / 24.0 / 60.0
second = minute / 60.0
jd = JulianDate(tt=arange(tt0, tt0 + 2 * minute, second * 5))
print jd.shape
(24,)
position_list = []
for tt in jd.tt:
position_list.append(positions_at(JulianDate(tt=tt)))
The simplest possible satellite display
uses a single mlab.points3d
call
to set up a graphics pipeline.
from tvtk.api import tvtk
style = tvtk.InteractorStyleTerrain()
fig = mlab.gcf()
fig.scene.interactor.interactor_style = style
mlab.clf()
def recenter():
mlab.view(distance=3e4, focalpoint=(0.0, 0.0, 0.0))
if take_screenshots:
imshow(mlab.screenshot(antialiased=True))
axis('off')
from mayavi.sources.builtin_surface import BuiltinSurface
def draw_globe():
sphere = mlab.points3d(0, 0, 0, name='Globe',
scale_mode='none', scale_factor=r * 2.0,
color=ocean_blue, resolution=50)
sphere.actor.property.specular = 0.20
sphere.actor.property.specular_power = 10
continents_src = BuiltinSurface(
source='earth', name='Continents')
continents_src.data_source.on_ratio = 1 # detail level
continents_src.data_source.radius = r
continents = mlab.pipeline.surface(
continents_src, color=(0, 0, 0))
draw_globe()
recenter()
x, y, z = position_list[0]
points_glyph = mlab.points3d(
x, y, z, x * 0.0 + 0.01,
mode='sphere', color=(1.0, 1.0, 0.8),
scale_mode='scalar', scale_factor=30000.0)
recenter()
One can interact with (that is, navigate in) the Mayavi window with mouse and keyboard as documented on the Using the Mayavi application web page.
for position in position_list:
x, y, z = position
radius = x * 0.0 + 0.01
points_glyph.mlab_source.set(x=x, y=y, z=z, s=radius)
This makes it easier to see the altitude of each satellite.
# Compute `distance` of each satellite
# from the center of the earth.
distance = np.sqrt((position * position).sum(axis=0))
xb, yb, zb = position / distance * r
print 'Satellite #0 is', distance[0], 'km from the geocenter'
print 'The base of its stalk will be at:'
print xb[0], yb[0], zb[0]
Satellite #0 is 6836.24046832 km from the geocenter The base of its stalk will be at: -2098.28789652 -5835.13265204 1492.9905943
x, y, z = position
src = mlab.pipeline.scalar_scatter(
np.hstack([xb, x]),
np.hstack([yb, y]),
np.hstack([zb, z]),
)
src.mlab_source.dataset.lines = np.array([
[i, i + len(x)] for i in range(len(x))
])
# The stripper filter cleans up connected lines
lines = mlab.pipeline.stripper(src)
# Finally, display the set of lines
mlab.pipeline.surface(lines, color=yellow, line_width=3,
opacity=0.4)
recenter()