Plot timeseries of equatorial crossing times using pyorbital.

ltan.py

"""Plot timeseries of equatorial crossing times using pyorbital."""

import pyorbital.orbital
from datetime import datetime, timedelta
from dateutil.rrule import rrule, MONTHLY
import matplotlib.pyplot as plt
import logging
import numpy as np
from matplotlib.dates import DateFormatter
from matplotlib.dates import HourLocator
import matplotlib.dates as mdates
import matplotlib.ticker as mticker


def tle2datetime64(times):
    """Convert TLE timestamps to numpy.datetime64

    Args:
       times (float): TLE timestamps as %y%j.1234, e.g. 18001.25
    """
    # Convert %y%j.12345 to %Y%j.12345 (valid for 1950-2049)
    times = np.where(times > 50000, times + 1900000, times + 2000000)

    # Convert float to datetime64
    doys = (times % 1000).astype('int') - 1
    years = (times // 1000).astype('int')
    msecs = np.rint(24 * 3600 * 1000 * (times % 1))
    times64 = (years - 1970).astype('datetime64[Y]').astype('datetime64[ms]')
    times64 += doys.astype('timedelta64[D]')
    times64 += msecs.astype('timedelta64[ms]')

    return times64


def get_tle_lines(tle_file, utc_time, tle_thresh=7):
    """Find closest two line elements (TLEs) for the current orbit
    
    Copied from pygac.

    Raises:
        IndexError, if the closest TLE is more than :meth:`pygac.GACReader.tle_thresh` days apart
    """
    with open(tle_file, 'r') as ifile:
        tle_data = ifile.readlines()

    sdate = np.datetime64(utc_time, 'ms')
    dates = tle2datetime64(
        np.array([float(line[18:32]) for line in tle_data[::2]]))

    # Find index "iindex" such that dates[iindex-1] < sdate <= dates[iindex]
    # Notes:
    #     1. If sdate < dates[0] then iindex = 0
    #     2. If sdate > dates[-1] then iindex = len(dates), beyond the right boundary!
    iindex = np.searchsorted(dates, sdate)

    if iindex in (0, len(dates)):
        if iindex == len(dates):
            # Reset index if beyond the right boundary (see note 2. above)
            iindex -= 1
    elif abs(sdate - dates[iindex - 1]) < abs(sdate - dates[iindex]):
        # Choose the closest of the two surrounding dates
        iindex -= 1

    # Make sure the TLE we found is within the threshold
    delta_days = abs(sdate - dates[iindex]) / np.timedelta64(1, 'D')
    if delta_days > tle_thresh:
        raise IndexError(
            "Can't find tle data within +/- %d days around %s" %
            (tle_thresh, sdate))

    # Select TLE data
    tle1 = tle_data[iindex * 2]
    tle2 = tle_data[iindex * 2 + 1]
    return tle1, tle2


if __name__ == '__main__':
    start_date = datetime(1998, 1, 1)
    end_date = datetime(2019, 1, 1)
    dates = list(rrule(dtstart=start_date, until=end_date, freq=MONTHLY))
    sats = ['NOAA-15', 'NOAA-16', 'NOAA-17', 'NOAA-18', 'NOAA-19', 'METOP-A', 'METOP-B']
    fig, ax = plt.subplots()

    for sat in sats:
        print(sat)
        tle_file = '/cmsaf/nfshome/sfinkens/software/devel/pygac/tle/TLE_{}.txt'.format(sat.replace('-', '').lower())

        ects = []
        ddates = []
        for date in dates:
            tstart = date
            tend = tstart + timedelta(hours=3)

            try:
                tle1, tle2 = get_tle_lines(tle_file, tstart)
            except IndexError:
                continue
    
            orb = pyorbital.orbital.Orbital(sat, line1=tle1, line2=tle2)
            ect = orb.get_equatorial_crossing_time(tstart, tend, rtol=1E-12, local_time=True)
            if ect is not None:
                ects.append(ect.time())
                ddates.append(date)

        ects_plot = [datetime.combine(date=datetime(2009, 7, 1).date(), time=ect)
                     for ect in ects]
        ax.plot(ddates, ects_plot, label=sat)

    ax.legend(loc='lower left')
    ax.set_xlabel('Year')
    ax.set_ylabel('Local Time of Ascending Node (Hour)')

    # ax.xaxis.set_major_locator(mdates.YearLocator())
    # ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))

    ax.yaxis.set_major_locator(mdates.HourLocator(byhour=range(0, 24, 1), interval=1))
    ax.yaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))

    plt.savefig('crossing_times.png')
    plt.show()