Interactive-3D-Earth-Satellites.ipynb

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