Merge pull request #106 from pysat/enh/skyfield

ENH: Use skyfield to propagate and convert

CHANGELOG.md

@@ -3,6 +3,7 @@ All notable changes to this project will be documented in this file.
 This project adheres to [Semantic Versioning](https://semver.org/).
 
 ## [0.3.4] - 2023-XX-XX
+* Add support for skyfield propagation
 * Maintenance
   * Update pytest syntax
   * Update Github Actions versions

README.md

@@ -36,6 +36,7 @@ Python 3.8+.
 | numpy          | pysat>=3.0.4      | aacgmv2          |
 | pandas         | pyEphem           | apexpy           |
 |                | sgp4>=2.7         | OMMBV            |
+|                | skyfield          |                  |
 
 
 One way to install is through pip.  Just type

docs/installation.rst

@@ -23,6 +23,7 @@ Python 3.8+ and pysat 3.0.4+.
   numpy            pysat>=3.0.4        aacgmv2
   pandas           pyEphem             apexpy
                    sgp4>=2.7           OMMBV>=1.0
+                   skyfield
  ================ =================== ==================
 
 

docs/supported_instruments.rst

@@ -16,3 +16,9 @@ Missions SGP4
 
 .. automodule:: pysatMissions.instruments.missions_sgp4
    :members:
+
+Missions Skyfield
+-----------------
+
+.. automodule:: pysatMissions.instruments.missions_skyfield
+  :members:

pyproject.toml

@@ -43,7 +43,8 @@ dependencies = [
   "pandas",
   "pysat >= 3.0.4",
   "pyEphem",
-  "sgp4 >= 2.7"
+  "sgp4 >= 2.7",
+  "skyfield"
 ]
 
 [project.optional-dependencies]

pysatMissions/instruments/__init__.py

@@ -3,5 +3,6 @@
 from pysatMissions.instruments import methods  # noqa: F401
 from pysatMissions.instruments import missions_ephem
 from pysatMissions.instruments import missions_sgp4
+from pysatMissions.instruments import missions_skyfield
 
-__all__ = ['missions_ephem', 'missions_sgp4']
+__all__ = ['missions_ephem', 'missions_sgp4', 'missions_skyfield']

pysatMissions/instruments/missions_skyfield.py

@@ -0,0 +1,316 @@
+# -*- coding: utf-8 -*-
+"""Simulate satellite orbit data using Two Line Elements (TLEs) and SGP4.
+
+Properties
+----------
+platform
+    'missions'
+name
+    'skyfield'
+tag
+    None supported
+inst_id
+    None supported
+
+"""
+
+import datetime as dt
+import functools
+import numpy as np
+import pandas as pds
+
+import pysat
+from pysat.instruments.methods import testing as ps_meth
+from pysatMissions.instruments import _core as mcore
+from pysatMissions.instruments.methods import orbits
+
+from sgp4 import api as sapi
+from skyfield import api as skyapi
+from skyfield import framelib
+
+# -------------------------------
+# Required Instrument attributes
+platform = 'missions'
+name = 'skyfield'
+tags = {'': 'Satellite simulation data set'}
+inst_ids = {'': ['']}
+
+# ------------------
+# Testing attributes
+_test_dates = {'': {'': dt.datetime(2018, 1, 1)}}
+
+
+def init(self):
+    """Initialize the Instrument object with required values.
+
+    Runs once upon instantiation.
+
+    """
+
+    orbits._check_orbital_params(self.kwargs)
+
+    self.acknowledgements = ' '.join((
+        'The project uses the skyfield library available at',
+        'https://github.com/skyfielders/python-skyfield'))
+    self.references = ' '.join((
+        'Rhodes, B., Skyfield: High precision research-grade positions for',
+        'planets and Earth satellites generator, Astrophysics Source Code',
+        'Library, 2019. ascl:1907.024.'))
+    pysat.logger.info(self.acknowledgements)
+
+    if 'epoch' not in self.kwargs['load'].keys():
+        self.kwargs['load']['epoch'] = self.files.files.index[0]
+
+    return
+
+
+# Clean method
+clean = mcore._clean
+
+
+def load(fnames, tag=None, inst_id=None, tle1=None, tle2=None,
+         alt_periapsis=None, alt_apoapsis=None,
+         inclination=None, raan=0., arg_periapsis=0., mean_anomaly=0.,
+         epoch=None, bstar=0., one_orbit=False, num_samples=None, cadence='1S'):
+    """Generate position of satellite in ECI co-ordinates.
+
+    Parameters
+    ----------
+    fnames : list-like collection
+        File name that contains date in its name.
+    tag : string
+        Identifies a particular subset of satellite data
+    inst_id : string
+        Instrument satellite ID (accepts '' or a number (i.e., '10'), which
+        specifies the number of seconds to simulate the satellite)
+        (default='')
+    tle1 : str or NoneType
+        First string for Two Line Element. Must be in TLE format.  tle1 and tle2
+        both required if instantiating instrument by TLEs. (defalt=None)
+    tle2 : str or NoneType
+        Second string for Two Line Element. Must be in TLE format.  tle1 and tle2
+        both required if instantiating instrument by TLEs. (default=None)
+    alt_periapsis : float or NoneType
+        The lowest altitude from the mean planet surface along the orbit (km).
+        Required along with inclination if instantiating via orbital elements.
+        (default=None)
+    alt_apoapsis : float or NoneType
+        The highest altitude from the mean planet surface along the orbit (km)
+        If None, assumed to be equal to periapsis (ie, circular orbit). Optional
+        when instantiating via orbital elements. (default=None)
+    inclination : float or NoneType
+        Orbital Inclination in degrees.  Required along with alt_periapsis if
+        instantiating via orbital elements. (default=None)
+    raan : float
+        Right Ascension of the Ascending Node (RAAN) in degrees. This defines
+        the orientation of the orbital plane to the generalized reference frame.
+        The Ascending Node is the point in the orbit where the spacecraft passes
+        through the plane of reference moving northward.  For Earth orbits, the
+        location of the RAAN is defined as the angle eastward of the First Point
+        of Aries. Optional when instantiating via orbital elements.
+        (default=0.)
+    arg_periapsis : float
+        Argument of Periapsis in degrees.  This defines the orientation of the
+        ellipse in the orbital plane, as an angle measured from the ascending
+        node to the periapsis.  Optional when instantiating via orbital elements.
+        (default=0.)
+    mean_anomaly : float
+        The fraction of an elliptical orbit's period that has elapsed since the
+        orbiting body passed periapsis.  Note that this is a "fictitious angle"
+        (input in degrees) which defines the location of the spacecraft in the
+        orbital plane based on the orbital period.  Optional when instantiating
+        via orbital elements.
+        (default=0.)
+    epoch : dt.datetime or NoneType
+        The epoch used for calculating orbital propagation from Keplerians.
+        If None, then use the first date in the file list for consistency across
+        multiple days. Note that this will be set in `init`. (default=None)
+    bstar : float
+        Inverse of the ballistic coefficient. Used to model satellite drag.
+        Measured in inverse distance (1 / earth radius).  Optional when
+        instantiating via orbital elements. (default=0.)
+    one_orbit : bool
+        Flag to override num_samples and only provide a single orbit.
+        (default=False)
+    num_samples : int
+        Number of samples per day.  (default=None)
+    cadence : str
+        Uses pandas.frequency string formatting ('1S', etc)
+        (default='1S')
+
+    Returns
+    -------
+    data : pandas.DataFrame
+        Object containing satellite data
+    meta : pysat.Meta
+        Object containing metadata such as column names and units
+
+    Note
+    ----
+    * Routine is directly called by pysat and not the user.
+    * Altitude accuracy expected to be on the order of 10 km in Low Earth Orbit.
+      Efforts to improve accuracy documented under issue #79.
+
+
+    Example
+    -------
+    ::
+
+          tle1='1 25544U 98067A   18135.61844383  .00002728  00000-0  48567-4 0  9998'
+          tle2='2 25544  51.6402 181.0633 0004018  88.8954  22.2246 15.54059185113452'
+          inst = pysat.Instrument('pysat', 'sgp4', tle1=tle1, tle2=tle2)
+          inst.load(2018, 1)
+
+    """
+
+    # TLEs (Two Line Elements for ISS)
+    # format of TLEs is fixed and available from wikipedia...
+    # lines encode list of orbital elements of an Earth-orbiting object
+    # for a given point in time
+    line1 = '1 25544U 98067A   18135.61844383  .00002728  00000-0  48567-4 0  9998'
+    line2 = '2 25544  51.6402 181.0633 0004018  88.8954  22.2246 15.54059185113452'
+
+    # If provided, use user-specified TLEs.  Otherwise use ISS defaults above.
+    if tle1 is not None:
+        line1 = tle1
+    if tle2 is not None:
+        line2 = tle2
+
+    if (num_samples is None) or one_orbit:
+        # Calculate one day of samples for default
+        num_samples = len(pds.date_range('2018/1/1', '2018/1/2',
+                                         freq=cadence)) - 1
+
+    # Extract list of times from filenames and inst_id
+    times, index, dates = ps_meth.generate_times(fnames, num_samples,
+                                                 freq=cadence)
+    # Calculate epoch for orbital propagator
+    epoch_days = (epoch - dt.datetime(1949, 12, 31)).days
+    jd, _ = sapi.jday(dates[0].year, dates[0].month, dates[0].day, 0, 0, 0)
+
+    if inclination is not None:
+        # If an inclination is provided, specify by Keplerian elements
+        eccentricity, mean_motion = orbits.convert_to_keplerian(alt_periapsis,
+                                                                alt_apoapsis)
+        satrec = sapi.Satrec()
+        satrec.sgp4init(sapi.WGS84, 'i', 0, epoch_days, bstar, 0, 0,
+                        eccentricity, np.radians(arg_periapsis),
+                        np.radians(inclination), np.radians(mean_anomaly),
+                        mean_motion, np.radians(raan))
+    else:
+        # Otherwise, use TLEs
+        satrec = sapi.Satrec.twoline2rv(line1, line2, sapi.WGS84)
+        mean_motion = satrec.mm
+
+    if one_orbit:
+        ind = times <= (2 * np.pi / mean_motion * 60)
+        times = times[ind]
+        index = index[ind]
+
+    skytimes = skyapi.load.timescale().from_datetimes(index.tz_localize('UTC'))
+
+    esat = skyapi.EarthSatellite.from_satrec(satrec, skyapi.load.timescale())
+    # Generate object with info for all times
+    sat_obj = esat.at(skytimes)
+    # ECI coords are True Equator and Equinox of date
+    position = sat_obj.position.km
+    velocity = sat_obj.velocity.km_per_s
+
+    # Add ECEF values to instrument.
+    # ECEF coords are International Terrestrial Reference System (ITRS)
+    ecef = sat_obj.frame_xyz_and_velocity(framelib.itrs)
+
+    # Convert to geocentric latitude, longitude, altitude.
+    lat, lon = skyapi.wgs84.latlon_of(sat_obj)
+    alt = skyapi.wgs84.height_of(sat_obj).km
+
+    # Put data into DataFrame
+    data = pds.DataFrame({'position_eci_x': position[0, :],
+                          'position_eci_y': position[1, :],
+                          'position_eci_z': position[2, :],
+                          'velocity_eci_x': velocity[0, :],
+                          'velocity_eci_y': velocity[1, :],
+                          'velocity_eci_z': velocity[2, :],
+                          'position_ecef_x': ecef[0].km[0],
+                          'position_ecef_y': ecef[0].km[1],
+                          'position_ecef_z': ecef[0].km[2],
+                          'velocity_ecef_x': ecef[1].km_per_s[0],
+                          'velocity_ecef_y': ecef[1].km_per_s[1],
+                          'velocity_ecef_z': ecef[1].km_per_s[2],
+                          'geod_latitude': lat.degrees,
+                          'geod_longitude': lon.degrees,
+                          'geod_altitude': alt},
+                         index=index)
+    data.index.name = 'Epoch'
+
+    # Create metadata corresponding to variables in load routine
+    meta = pysat.Meta()
+    meta['Epoch'] = {
+        meta.labels.units: 'Milliseconds since 1970-1-1',
+        meta.labels.notes: 'UTC time at middle of geophysical measurement.',
+        meta.labels.desc: 'UTC seconds',
+        meta.labels.name: 'Time index in milliseconds'}
+    for v in ['x', 'y', 'z']:
+        meta['position_eci_{:}'.format(v)] = {
+            meta.labels.units: 'km',
+            meta.labels.name: 'ECI {:}-position'.format(v),
+            meta.labels.notes: 'Calculated using skyfield with TEME',
+            meta.labels.desc: 'Earth Centered Inertial {:}-position'.format(v),
+            meta.labels.min_val: -np.inf,
+            meta.labels.max_val: np.inf,
+            meta.labels.fill_val: np.nan}
+        meta['velocity_eci_{:}'.format(v)] = {
+            meta.labels.units: 'km/s',
+            meta.labels.name: 'Satellite velocity ECI-{:}'.format(v),
+            meta.labels.notes: 'Calculated using skyfield with TEME',
+            meta.labels.desc: 'Satellite velocity along ECI-{:}'.format(v),
+            meta.labels.min_val: -np.inf,
+            meta.labels.max_val: np.inf,
+            meta.labels.fill_val: np.nan}
+        meta['position_ecef_{:}'.format(v)] = {
+            meta.labels.units: 'km',
+            meta.labels.notes: 'Calculated using skyfield with ITRS',
+            meta.labels.name: 'ECEF {:}-position'.format(v),
+            meta.labels.desc: 'Earth Centered Earth Fixed {:}-position'.format(v),
+            meta.labels.min_val: -np.inf,
+            meta.labels.max_val: np.inf,
+            meta.labels.fill_val: np.nan}
+        meta['velocity_ecef_{:}'.format(v)] = {
+            meta.labels.units: 'km',
+            meta.labels.notes: 'Calculated using skyfield with ITRS',
+            meta.labels.name: 'ECEF {:}-velocity'.format(v),
+            meta.labels.desc: 'Earth Centered Earth Fixed {:}-velocity'.format(v),
+            meta.labels.min_val: -np.inf,
+            meta.labels.max_val: np.inf,
+            meta.labels.fill_val: np.nan}
+    meta['geod_latitude'] = {
+        meta.labels.units: 'degrees',
+        meta.labels.notes: 'Calculated using skyfield with WGS84',
+        meta.labels.name: 'Geodetic Latitude',
+        meta.labels.desc: 'Geodetic Latitude of satellite',
+        meta.labels.min_val: -90.,
+        meta.labels.max_val: 90.,
+        meta.labels.fill_val: np.nan}
+    meta['geod_longitude'] = {
+        meta.labels.units: 'degrees',
+        meta.labels.notes: 'Calculated using skyfield with WGS84',
+        meta.labels.name: 'Geodetic Longitude',
+        meta.labels.desc: 'Geodetic Longitude of satellite',
+        meta.labels.min_val: -180.,
+        meta.labels.max_val: 180.,
+        meta.labels.fill_val: np.nan}
+    meta['geod_altitude'] = {
+        meta.labels.units: 'km',
+        meta.labels.notes: 'Calculated using skyfield with WGS84',
+        meta.labels.name: 'Altitude',
+        meta.labels.desc: 'Altitude of satellite above an ellipsoidal Earth',
+        meta.labels.min_val: -np.inf,
+        meta.labels.max_val: np.inf,
+        meta.labels.fill_val: np.nan}
+
+    return data, meta
+
+
+list_files = functools.partial(ps_meth.list_files, test_dates=_test_dates)
+download = functools.partial(ps_meth.download)
+clean = functools.partial(mcore._clean)

pysatMissions/tests/test_instruments.py

@@ -33,6 +33,8 @@
 for inst in instruments['download']:
     if 'sgp4' in inst['inst_module'].name:
         instruments['sgp4'].append(inst)
+    if 'skyfield' in inst['inst_module'].name:
+        instruments['sgp4'].append(inst)
 
 
 class TestInstruments(clslib.InstLibTests):

requirements.txt

@@ -4,3 +4,4 @@ pandas
 pysat>=3.0.4
 pyEphem
 sgp4>=2.7
+skyfield