Merge pull request cbassa#22 from kerel-fs/idea/ground_track

WIP: Idea/ground track

examples/irvine01.txt

@@ -0,0 +1,3 @@
+IRVINE01
+1 70002U 18599A   18315.16151858  .00000000  00000-0  00000-0 0    07
+2 70002  85.1205  90.1568 0012705 292.5520 107.9249 15.20792276    04

examples/iss.txt

@@ -0,0 +1,3 @@
+ISS (ZARYA)
+1 25544U 98067A   19004.59354167  .00000715  00000-0  18267-4 0  9995
+2 25544  51.6416  95.0104 0002419 236.2184 323.8248 15.53730729149833

examples/sathyabamasat.txt

@@ -0,0 +1,3 @@
+SATHYABAMASAT
+1 41600U 16040B   19310.17193730  .00001455  00000-0  63216-4 0  9994
+2 41600  97.3537  10.1403 0012192 255.4467 104.5418 15.23822105187350

python/README.md

@@ -0,0 +1,34 @@
+## Notes
+
+
+## Setup
+- Dependency installation via pip:
+  ```
+  pip install -r requirements.txt
+  ```
+
+- In order to be able to use `ground_track.py`,
+  download [earth.png](https://raw.githubusercontent.com/galactics/beyond/master/doc/source/_static/earth.png) into this directory.
+
+## Usage
+
+- Plot the ground track of satellite based on it's TLE around the epoch:
+  ```
+  ./ground_track.py ../examples/sathyabamasat.txt
+  ```
+
+- Calculate and plot satellite passes (requires a site.txt)
+  ```
+  ./pass.py ../examples/sathyabamasat.txt $ST_DATADIR/sites.txt -s 7300 --starttime 2019-11-06T19:30:00
+  ```
+
+- Adjust a TLE to a new launch time
+  ```
+  # Data from https://community.libre.space/t/electron-its-business-time-launch-this-weekend-irvine-01-amateur-payload/2819/4
+  $ ./launchtle.py ../examples/irvine01.txt 2018-11-11T03:00:00 2018-11-11T04:05:00
+  New TLE for a launch delay of 1:05:00h :
+
+  IRVINE01-delayed
+  1 70002U 18599A   18315.20665747  .00000000  00000-0  00000-0 0    09
+  2 70002  85.1205 106.4513 0012705 292.5520 107.9249 15.20792276    05
+  ```

python/ground_track.py

@@ -0,0 +1,114 @@
+#!/usr/bin/env python3
+
+"""
+Script showing the position of the ISS at the time of the TLE
+and the ground track for the previous and the next orbit
+
+(original) source: https://github.com/galactics/beyond/blob/master/doc/source/_static/ground-track.py
+"""
+
+import sys
+import argparse
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pathlib import Path
+
+from beyond.dates import Date, timedelta
+from utils import tle_from_file
+
+
+if __name__ == "__main__":
+    # Parse arguments
+    parser = argparse.ArgumentParser(
+        description="Plot the ground track of a TLE")
+    parser.add_argument('TLEFILE', help="the tle file", type=str)
+
+    args = parser.parse_args()
+
+    # Parsing of TLE
+    tle, name = tle_from_file(args.TLEFILE)
+
+    # Conversion into `Orbit` object
+    orb = tle.orbit()
+
+    # Tables containing the positions of the ground track
+    latitudes, longitudes = [], []
+    prev_lon, prev_lat = None, None
+
+    period = orb.infos.period
+    # start = orb.date - period
+    start = Date.now()
+    stop = period
+    step = period / 100
+
+    for point in orb.ephemeris(start=start, stop=stop, step=step):
+
+        # Conversion to earth rotating frame
+        point.frame = 'ITRF'
+
+        # Conversion from cartesian to spherical coordinates (range, latitude, longitude)
+        point.form = 'spherical'
+
+        # Conversion from radians to degrees
+        lon, lat = np.degrees(point[1:3])
+
+        # Creation of multiple segments in order to not have a ground track
+        # doing impossible paths
+        if prev_lon is None:
+            lons = []
+            lats = []
+            longitudes.append(lons)
+            latitudes.append(lats)
+        elif orb.i < np.pi /2 and (np.sign(prev_lon) == 1 and np.sign(lon) == -1):
+            lons.append(lon + 360)
+            lats.append(lat)
+            lons = [prev_lon - 360]
+            lats = [prev_lat]
+            longitudes.append(lons)
+            latitudes.append(lats)
+        elif orb.i > np.pi/2 and (np.sign(prev_lon) == -1 and np.sign(lon) == 1):
+            lons.append(lon - 360)
+            lats.append(lat)
+            lons = [prev_lon + 360]
+            lats = [prev_lat]
+            longitudes.append(lons)
+            latitudes.append(lats)
+
+        lons.append(lon)
+        lats.append(lat)
+        prev_lon = lon
+        prev_lat = lat
+
+
+    fig = plt.figure(figsize=(15.2, 8.2))
+    ax = fig.add_subplot(111)
+
+    ax.set_title(start.datetime.strftime('%Y-%m-%dT%H:%M:%SZ'))
+
+    # Plot earth
+    img = Path(__file__).parent / "earth.png"
+    im = plt.imread(str(img))
+    ax.imshow(im, extent=[-180, 180, -90, 90])
+
+    # Plot ground track 
+    for lons, lats in zip(longitudes, latitudes):
+        ax.plot(lons, lats, 'r-')
+
+    # Plot current location
+        ax.plot(lons[0], lats[0], 'bo')
+
+    # Plot location at epoch
+    # lon, lat = np.degrees(orb.copy(frame='ITRF', form='spherical')[1:3])
+    # ax.plot([lon], [lat], 'bo')
+
+    ax.set_xlim([-180, 180])
+    ax.set_ylim([-90, 90])
+    ax.grid(True, color='w', linestyle=":", alpha=0.4)
+    ax.set_xticks(range(-180, 181, 30))
+    ax.set_yticks(range(-90, 91, 30))
+    fig.tight_layout()
+
+    if "no-display" not in sys.argv:
+        plt.show()

python/launchtle.py

@@ -0,0 +1,84 @@
+#!/usr/bin/env python3
+
+import argparse
+from datetime import datetime, timedelta
+
+from sgp4.io import twoline2rv
+from sgp4.earth_gravity import wgs84
+
+from astropy.coordinates import Longitude, Angle
+from astropy import units as u
+from astropy.time import Time
+
+
+
+def gmst(t):
+    return Time(t).sidereal_time('mean', 'greenwich')
+
+def fractional_days(t):
+    d = t - datetime(t.year, 1, 1)
+    return t.timetuple().tm_yday + (d.seconds + d.microseconds/1e6)/(60*60*24)
+
+def print_tle(tle, new_epoch, new_nodeo):
+    def checksum(proto_tle):
+        s = 0
+        for c in proto_tle:
+            if c.isdigit():
+                s += int(c)
+            if c == '-':
+                s += 1
+        return s%10
+
+    tle0,old1,old2 = tle
+    tle1_proto = '{} {:2d}{:.8f} {}'.format(old1[:17],
+                              abs(new_epoch.year)%100,
+                              fractional_days(new_epoch),
+                              old1[33:-1])
+    tle2_proto = '{} {:>8.4f} {}'.format(old2[:16],
+                                   Angle(new_nodeo*u.radian).degree,
+                                   old2[26:68])
+
+    return (tle0,
+            tle1_proto + str(checksum(tle1_proto)),
+            tle2_proto + str(checksum(tle2_proto)))
+
+def launch_tle(tle, launch_date, new_launch_date):
+    sat = twoline2rv(tle[1], tle[2], whichconst=wgs84)
+
+    # New Epoch
+    new_epoch = sat.epoch - launch_date + new_launch_date
+
+    # New Right ascension of ascending node in radians
+    new_nodeo = (gmst(new_launch_date) - gmst(launch_date) + Longitude(sat.nodeo * u.radian)).rad
+
+    return print_tle(tle, new_epoch, new_nodeo)
+
+
+if __name__ == "__main__":
+    # Parse arguments
+    parser = argparse.ArgumentParser(
+        description="Adjust a TLE to a new launch time")
+    parser.add_argument('TLEFILE', help="reference tle file", type=str)
+    parser.add_argument("OLD",
+                        help="Reference (old) launch time (YYYY-MM-DDTHH:MM:SS)",
+                        type=str)
+    parser.add_argument("NEW",
+                        help="New launch time (YYYY-MM-DDTHH:MM:SS) [default: now]",
+                        default=datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%S"),
+                        type=str)
+
+    args = parser.parse_args()
+
+    with open(args.TLEFILE, 'r') as f:
+        tle = f.readlines()
+    tle = [l.strip() for l in tle]
+
+    launch_date_ref = datetime.strptime(args.OLD, '%Y-%m-%dT%H:%M:%S')
+    new_launch_date = datetime.strptime(args.NEW, '%Y-%m-%dT%H:%M:%S')
+
+    new_tle = launch_tle(tle, launch_date_ref, new_launch_date)
+
+    print('New TLE for a launch delay of {}h :\n'.format(new_launch_date - launch_date_ref))
+    print(new_tle[0].strip()+'-delayed')
+    print(new_tle[1])
+    print(new_tle[2])