Implementing custom satellites+data in fisheye animations when working with custom satellite and data

[ksd3]

Two-line element data of scientific satellites are often published online and users can derive LLA positions as functions of time of the satellites based on this information. If I want to plot the trajectory of a satellite and the custom data it records on an ASI image, how would I go about doing this? In the example given,

from datetime import datetime

import numpy as np
import matplotlib.pyplot as plt

import asilib
import asilib.asi

# ASI parameters
location_code = 'RANK'
alt=110  # km
time_range = (datetime(2017, 9, 15, 2, 32, 0), datetime(2017, 9, 15, 2, 35, 0))

fig, ax = plt.subplots(
    3, 1, figsize=(7, 10), gridspec_kw={'height_ratios': [4, 1, 1]}, constrained_layout=True
)

asi = asilib.asi.themis(location_code, time_range=time_range, alt=alt)

# Create the fake satellite track coordinates: latitude, longitude, altitude (LLA).
# This is a north-south satellite track oriented to the east of the THEMIS/RANK
# imager.
n = int((time_range[1] - time_range[0]).total_seconds() / 3)  # 3 second cadence.
lats = np.linspace(asi.meta["lat"] + 5, asi.meta["lat"] - 5, n)
lons = (asi.meta["lon"] - 0.5) * np.ones(n)
alts = alt * np.ones(n)  # Altitude needs to be the same as the skymap.
sat_lla = np.array([lats, lons, alts]).T
# Normally the satellite time stamps are not the same as the ASI.
# You may need to call Conjunction.interp_sat() to find the LLA coordinates
# at the ASI timestamps.
sat_time = asi.data.time

conjunction_obj = asilib.Conjunction(asi, (sat_time, sat_lla))

# Map the satellite track to the imager's azimuth and elevation coordinates and
# image pixels. NOTE: the mapping is not along the magnetic field lines! You need
# to install IRBEM and then use conjunction.lla_footprint() before
# calling conjunction_obj.map_azel.
sat_azel, sat_azel_pixels = conjunction_obj.map_azel()

# Calculate the auroral intensity near the satellite and mean intensity within a 10x10 km area.
nearest_pixel_intensity = conjunction_obj.intensity(box=None)
area_intensity = conjunction_obj.intensity(box=(10, 10))
area_mask = conjunction_obj.equal_area(box=(10,10))

# Need to change masked NaNs to 0s so we can plot the rectangular area contours.
area_mask[np.where(np.isnan(area_mask))] = 0

# Initiate the animation generator function.
gen = asi.animate_fisheye_gen(
    ax=ax[0], azel_contours=True, overwrite=True, cardinal_directions='NE'
)

for i, (time, image, _, im) in enumerate(gen):
    # Plot the entire satellite track, its current location, and a 20x20 km box
    # around its location.
    ax[0].plot(sat_azel_pixels[:, 0], sat_azel_pixels[:, 1], 'red')
    ax[0].scatter(sat_azel_pixels[i, 0], sat_azel_pixels[i, 1], c='red', marker='o', s=50)
    ax[0].contour(area_mask[i, :, :], levels=[0.99], colors=['yellow'])

    if 'vline1' in locals():
        vline1.remove()  # noqa: F821
        vline2.remove()  # noqa: F821
        text_obj.remove()  # noqa: F821
    else:
        # Plot the ASI intensity along the satellite path
        ax[1].plot(sat_time, nearest_pixel_intensity)
        ax[2].plot(sat_time, area_intensity)
    vline1 = ax[1].axvline(time, c='b')
    vline2 = ax[2].axvline(time, c='b')

    # Annotate the location_code and satellite info in the top-left corner.
    location_code_str = (
        f'THEMIS/{location_code} '
        f'LLA=({asi.meta["lat"]:.2f}, '
        f'{asi.meta["lon"]:.2f}, {asi.meta["alt"]:.2f})'
    )
    satellite_str = f'Satellite LLA=({sat_lla[i, 0]:.2f}, {sat_lla[i, 1]:.2f}, {sat_lla[i, 2]:.2f})'
    text_obj = ax[0].text(
        0,
        1,
        location_code_str + '\n' + satellite_str,
        va='top',
        transform=ax[0].transAxes,
        color='red',
    )
    ax[1].set(ylabel='ASI intensity\nnearest pixel [counts]')
    ax[2].set(xlabel='Time', ylabel='ASI intensity\n10x10 km area [counts]')

print(f'Animation saved in {asilib.config["ASI_DATA_DIR"] / "animations" / asi.animation_name}')

there does not seem to be a convenient way to include an extra subplot. Trying to arbitrarily insert ax[3] leads to indexing errors and the code for the fisheye animation generator shows that ax.plots is set to None by default. Is there a way for a user to insert custom data? I've been doing it by replacing the plotted area_intensity with my own data, but that requires using functions like np.tile in order to ensure that my data is in the correct form.

Is there a better way to do it?

For custom satellite data: I have been downsampling my satellite trajectory (LLA coords) and passing it to asilib.Conjunction, but I would like to know what the intended usage of interp_sat() is. Is it intended for use on position data where the original data is sampled at an extremely high rate? Of course, there are inaccuracies, but I want to know if asilib has an inbuilt method to handle this. However this is not a huge problem.

[mshumko]

Hi @ksd3,

asilib works with propagated geodetic ephemeris, so you'll have to run SGP4 on the TLEs, and pass that into the asilib.Conjunction() class. The skyfield's TLE propagator is probably the easiest approach, with another one being the python implementation of the SGP4 algorithm.

If you want to add another panel to the plot, you can change the plt.subplots() line to:

fig, ax = plt.subplots(
    4, 1, figsize=(7, 10), gridspec_kw={'height_ratios': [4, 1, 1, 1]}, constrained_layout=True
)

And you are correct, the asilib.Conjunction().interp_sat() interpolates the satellite's coordinates to the imager time stamps (at a three second cadence for the imager arrays currently supported by asilib).