{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Working with satellites" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## License\n", "\n", "```\n", "Working with the satellites module of pycraf.\n", "Copyright (C) 2015+ Benjamin Winkel (bwinkel@mpifr.de)\n", "\n", "This program is free software; you can redistribute it and/or\n", "modify it under the terms of the GNU General Public License\n", "as published by the Free Software Foundation; either version 2\n", "of the License, or (at your option) any later version.\n", "\n", "This program is distributed in the hope that it will be useful,\n", "but WITHOUT ANY WARRANTY; without even the implied warranty of\n", "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n", "GNU General Public License for more details.\n", "\n", "You should have received a copy of the GNU General Public License\n", "along with this program; if not, write to the Free Software\n", "Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note: GitHub's integrated notebook viewer won't show the embedded videos. You can use the [nbviewer](http://nbviewer.jupyter.org/) service to open the notebook. Or follow the direct [link](http://nbviewer.jupyter.org/github/bwinkel/pycraf/blob/master/notebooks/04a_satellites_plot.ipynb)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import requests\n", "import numpy as np\n", "\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D\n", "from matplotlib import animation, rc\n", "\n", "from pycraf import satellite\n", "from astropy import time\n", "from astropy import units as u\n", "from astropy.coordinates import EarthLocation\n", "\n", "rc('animation', html='html5')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Two-line elements sets (TLE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To compute positions of a satellite, two-line element sets (TLE) have to be used. These look like the following\n", "\n", " ISS (ZARYA)\n", " 1 25544U 98067A 13165.59097222 .00004759 00000-0 88814-4 0 47\n", " 2 25544 51.6478 121.2152 0011003 68.5125 263.9959 15.50783143834295\n", "\n", "which is the latest TLE for the international space station. The first line has only the name (or identifier) of a satellite. The following two lines contain all the necessary information about the satellite orbit to calculate its position for a certain time. Note that the TLEs are usually published once a day, because the contained parameters quickly change; drag forces cause rapid changes in the orbits of almost all satellites." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course, the position of a satellite depends on the choosen coordinate system. With the help of the [sgp4 package](https://pypi.python.org/pypi/sgp4/) pycraf calculates the ECI position of the satellite (in Cartesian coordinates). ECI is the [Earth-centered inertial frame](https://en.wikipedia.org/wiki/Earth-centered_inertial). But often we are interested in the apparent sky position (horizontal frame) of a satellite with respect to an observer. The latter is not provided by the `sgp4` library, but `pycraf` adds this functionality." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fetching latest satellites from Celestrak" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an example how to use the `pycraf.satellites` module, we will download the latest TLE data of all satellites launched in the last 30 days." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "celestrak_latest = requests.get('http://celestrak.com/NORAD/elements/tle-new.txt')" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "'FENGYUN 4B \\r\\n1 48808U 21047A 21179.92803759 -.00000367 00000-0 00000-0 0 9992\\r\\n2 48808 0.0541 226.6478 0000694 252.0694 256.3446 1.00272877 387\\r\\nCZ-3B R/B \\r\\n1 48809U 21047B 21179.71501365 .00048070 97412-6 15215-2 0 9995\\r\\n2 48809 28.6217 122.6465 7215152 195.7071 116.8720 2.40687722 638\\r\\nDRAGON CRS-22 \\r\\n1 48831U 21048A 21179.59233343 .00001944 00000-0 43787-4 0 9992\\r\\n2 48831 51.6454 286.3205 0002279 118.9439 5.7591 15.48743536290365\\r\\nSXM-8 \\r\\n1 48838U 21049A 21180.57409069 .00000001 00000-0 00000+0 0 9999\\r\\n2 48838 0.0430 263.7940 0002147 166.6985 293.5692 1.00273981 441\\r\\nFALCON 9 R/B \\r\\n1 48839U 21049B 21180.25559672 .00006015 00000-0 42660-3 0 9994\\r\\n2 48839 27.0231 135.4209 5912344 205.6044 106.4870 4.22186156 979\\r\\n2021-050A \\r\\n1 48840U 21050A 21180.43817509 -.00000073 00000-0 00000+0 0 9991\\r\\n2 48840 97.5118 250.1699 0011777 203.6227 211.1975 15.23623308'" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "celestrak_latest.text[:1000]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 'tle-new.txt' file apparently simply contains a list of TLEs. Let's fix the line endings (`\\r\\n` to `\\n`) and split into a list of TLEs:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "all_lines = celestrak_latest.text.split('\\r\\n')\n", "tle_list = []\n", "for idx in range(len(all_lines) // 3):\n", " tle_list.append('\\n'.join(all_lines[3 * idx: 3 * (idx + 1)]))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "FENGYUN 4B \n", "1 48808U 21047A 21179.92803759 -.00000367 00000-0 00000-0 0 9992\n", "2 48808 0.0541 226.6478 0000694 252.0694 256.3446 1.00272877 387\n", "CZ-3B R/B \n", "1 48809U 21047B 21179.71501365 .00048070 97412-6 15215-2 0 9995\n", "2 48809 28.6217 122.6465 7215152 195.7071 116.8720 2.40687722 638\n" ] } ], "source": [ "print(tle_list[0])\n", "print(tle_list[1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting satellite position" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we start plotting, we'll use have `sgp4` parse all the satellite strings and construct `sgp4.io.Satellite` instances. These merely contain the orbital parameters, but it will be faster to work the the already parsed data instead of letting the TLEs being parsed in each time step (yes, we are going to do a nice animation)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "('FENGYUN 4B ', )" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sat_list = [satellite.get_sat(tle_string) for tle_string in tle_list]\n", "# contains a list of tuples: (sat_name, sat_instance)\n", "sat_list[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The second ingredient, we need is an array of \"observing\" times. pycraf wants an `astropy.time.Time` object for this, because it has (accurate!) built-in conversion between UTC and sidereal time." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Start time: 2021-06-29 18:32:03.306425\n" ] } ], "source": [ "start_time = time.Time.now()\n", "print('Start time:', start_time)\n", "td = time.TimeDelta(np.arange(0, 3600 * 24, 60 * 1), format='sec')\n", "times = start_time + td # one day in steps of 1 min" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### ECI position\n", "For this we make use of `sgp4` built-in `propagate` method. This is easy to use, but it doesn't work with numpy arrays. Let's vectorize it:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def _propagate(sat, dt):\n", " '''\n", " True equator mean equinox (TEME) position from `sgp4` at given time.\n", "\n", " Parameters\n", " ----------\n", " sat : `sgp4.io.Satellite` instance\n", " Satellite object filled from TLE\n", " dt : `~datetime.datetime`\n", " Time\n", "\n", " Returns\n", " -------\n", " xs, ys, zs : float\n", " TEME (=True equator mean equinox) position of satellite [km]\n", " '''\n", "\n", " # pos [km], vel [km/s]\n", "\n", " try:\n", " from sgp4.api import jday, SGP4_ERRORS\n", " except ImportError:\n", " raise ImportError(\n", " 'The \"sgp4\" package (version 2+) is necessary to use this '\n", " 'function.'\n", " )\n", "\n", " jd, fr = jday(\n", " dt.year, dt.month, dt.day,\n", " dt.hour, dt.minute, dt.second + dt.microsecond / 1e6\n", " )\n", " err_code, position, velocity = sat.sgp4(jd, fr)\n", "\n", " if err_code:\n", " raise ValueError(\n", " 'Satellite propagation error', err_code,\n", " '({})'.format(SGP4_ERRORS[err_code])\n", " )\n", "\n", " return position\n", "\n", "\n", "vec_propagate = np.vectorize(\n", " _propagate, excluded=['sat'], otypes=[np.float64] * 3\n", " )" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "pos_eci = np.array([vec_propagate(sat[1], times.datetime) for sat in sat_list])" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plim = 6000\n", "\n", "fig = plt.figure(figsize=(12, 9))\n", "ax = Axes3D(fig)\n", "ax.view_init(azim=60, elev=30)\n", "# Aspect ratio is not implemented;\n", "# see https://github.com/matplotlib/matplotlib/issues/1077/\n", "# ax.set_aspect('equal')\n", "# need to manually stretch to make it approx. right\n", "ax.set_xlim((-2 * plim, 2 * plim))\n", "ax.set_ylim((-2 * plim, 2 * plim))\n", "ax.set_zlim((-plim, plim))\n", "ax.set_xlabel('x [km]')\n", "ax.set_ylabel('y [km]')\n", "ax.set_zlabel('z [km]')\n", "\n", "points = ax.scatter(pos_eci[:, 0, 0], pos_eci[:, 1, 0], pos_eci[:, 2, 0], marker='o')\n", "title = ax.set_title('{:%y/%m/%d %H:%M}'.format(times[0].datetime), loc='center', fontsize=20)\n", "\n", "def init():\n", " points._offsets3d = (pos_eci[:, 0, 0], pos_eci[:, 1, 0], pos_eci[:, 2, 0])\n", " title.set_text('{:%y/%m/%d %H:%M}'.format(times[0].datetime))\n", " return points, title\n", "\n", "def animate(i):\n", " points._offsets3d = (pos_eci[:, 0, i], pos_eci[:, 1, i], pos_eci[:, 2, i])\n", " title.set_text('{:%y/%m/%d %H:%M}'.format(times[i].datetime))\n", " return points, title\n", "\n", "# call the animator. blit=True means only re-draw the parts that have changed.\n", "anim = animation.FuncAnimation(\n", " fig, animate, init_func=init, frames=pos_eci.shape[2], interval=20, blit=True\n", " )\n", "\n", "# this takes a while!\n", "plt.close(anim._fig)\n", "anim" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Horizontal position" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the next plot, we need to define an observer position, say the famous Parkes 64-m dish and the Effelsberg 100-m radio telescope." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "obs_loc_parkes = EarthLocation(148.25738, -32.9933, 414.8)\n", "obs_loc_effbg = EarthLocation(6.88375, 50.525, 366.)\n", "\n", "# create a SatelliteObserver instance\n", "sat_obs_parkes = satellite.SatelliteObserver(obs_loc_parkes)\n", "sat_obs_effbg = satellite.SatelliteObserver(obs_loc_effbg)\n", "pos_horiz_parkes = np.array([\n", " sat_obs_parkes.azel_from_sat(sat[1], times)\n", " for sat in sat_list\n", " ])\n", "pos_horiz_effbg = np.array([\n", " sat_obs_effbg.azel_from_sat(sat[1], times)\n", " for sat in sat_list\n", " ])\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vmin, vmax = np.log10(100), np.log10(100000)\n", "\n", "fig = plt.figure(figsize=(12, 7))\n", "ax1 = fig.add_axes((0.1, 0.5, 0.8, 0.35))\n", "ax2 = fig.add_axes((0.1, 0.1, 0.8, 0.35))\n", "cax = fig.add_axes((0.91, 0.2, 0.02, 0.5))\n", "ax2.set_xlabel('Azimuth [deg]')\n", "ax1.set_ylabel('Elevation [deg]')\n", "for ax in [ax1, ax2]:\n", " ax.set_xlim((-180, 180))\n", " ax.set_ylim((0, 90))\n", " ax.set_xticks(range(-150, 180, 30))\n", " ax.set_yticks(range(0, 91, 30))\n", " ax.set_aspect('equal')\n", " ax.grid()\n", "\n", "points1 = ax1.scatter(\n", " pos_horiz_parkes[:, 0, 0], pos_horiz_parkes[:, 1, 0],\n", " c=np.log10(pos_horiz_parkes[:, 2, 0]),\n", " cmap='viridis', vmin=vmin, vmax=vmax,\n", " )\n", "points2 = ax2.scatter(\n", " pos_horiz_effbg[:, 0, 0], pos_horiz_effbg[:, 1, 0],\n", " c=np.log10(pos_horiz_effbg[:, 2, 0]),\n", " cmap='viridis', vmin=vmin, vmax=vmax,\n", " )\n", "cbar = fig.colorbar(points1, cax=cax)\n", "cbar.set_label('Distance (km)')\n", "cbar.set_ticks([2, 3, 4, 5])\n", "cbar.set_ticklabels([100, 1000, 10000, 100000])\n", "\n", "ax1.text(-170, 75, 'Parkes 64-m', fontsize=16)\n", "ax2.text(-170, 75, 'Effelsberg 100-m', fontsize=16)\n", "title = ax1.text(\n", " 174, 75, '{:%y/%m/%d %H:%M}'.format(times[0].datetime),\n", " fontsize=15, ha='right'\n", " )\n", "\n", "def init():\n", " points1.set_offsets(pos_horiz_parkes[:, 0:2, 0])\n", " points1.set_array(np.log10(pos_horiz_parkes[:, 2, 0]))\n", " points2.set_offsets(pos_horiz_effbg[:, 0:2, 0])\n", " points2.set_array(np.log10(pos_horiz_effbg[:, 2, 0]))\n", " title.set_text('{:%y/%m/%d %H:%M}'.format(times[0].datetime))\n", " return points, title\n", "\n", "def animate(i):\n", " points1.set_offsets(pos_horiz_parkes[:, 0:2, i])\n", " points1.set_array(np.log10(pos_horiz_parkes[:, 2, i]))\n", " points2.set_offsets(pos_horiz_effbg[:, 0:2, i])\n", " points2.set_array(np.log10(pos_horiz_effbg[:, 2, i]))\n", " title.set_text('{:%y/%m/%d %H:%M}'.format(times[i].datetime))\n", " return points, title\n", "\n", "anim = animation.FuncAnimation(\n", " fig, animate, init_func=init, frames=pos_horiz_parkes.shape[2], interval=20, blit=True\n", " )\n", "\n", "# this takes a while!\n", "plt.close(anim._fig)\n", "anim" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python [conda root]", "language": "python", "name": "conda-root-py" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.2" } }, "nbformat": 4, "nbformat_minor": 4 }