Conjunctions#

In this advanced tutorial we will combine a low Earth orbiting satellite ephemeris with ASIs, to calculate the mean auroral intensity at the satellite footprint as a function of time.

The mean auroral intensity is calculated two ways: 1) the nearest pixel, and 2) a (20x20) km area at 110 km emission altitude.

[2]:

from datetime import datetime, timedelta
from IPython.display import Video
import numpy as np

import matplotlib.pyplot as plt
import asilib
import asilib.asi
import asilib.map

plt.style.use('dark_background')

print(f'asilib version: {asilib.__version__}')

asilib version: 0.20.3

[3]:

location_code = 'RANK'
alt=110
time_range = (datetime(2017, 9, 15, 2, 32, 0), datetime(2017, 9, 15, 2, 35, 0))
asi = asilib.asi.themis(location_code, time_range=time_range, alt=alt)

Downloaded 20170915_0232_rank_themis12_full.pgm.gz.
Downloaded 20170915_0233_rank_themis12_full.pgm.gz.
Downloaded 20170915_0234_rank_themis12_full.pgm.gz.
Downloaded themis_skymap_rank_20061101-20071203_vXX.sav.
Downloaded themis_skymap_rank_20071212-%2B_vXX.sav.
Downloaded themis_skymap_rank_20071218-%2B_vXX.sav.
Downloaded themis_skymap_rank_20071207-20090428_vXX.sav.
Downloaded themis_skymap_rank_20090911-%2B_vXX.sav.
Downloaded themis_skymap_rank_20100909-%2B_vXX.sav.
Downloaded themis_skymap_rank_20110202-%2B_vXX.sav.
Downloaded themis_skymap_rank_20120225-%2B_vXX.sav.
Downloaded themis_skymap_rank_20130107-%2B_vXX.sav.
Downloaded themis_skymap_rank_20150825-%2B_vXX.sav.
Downloaded themis_skymap_rank_20150825-20171003_vXX.sav.
Downloaded themis_skymap_rank_20170915-%2B_v02.sav.
Downloaded themis_skymap_rank_20170915-%2B_vXX.sav.
Downloaded themis_skymap_rank_20171003-%2B_v02.sav.
Downloaded themis_skymap_rank_20171003-%2B_vXX.sav.
Downloaded themis_skymap_rank_20190206-%2B_v02.sav.
Downloaded themis_skymap_rank_20190206-%2B_vXX.sav.
Downloaded themis_skymap_rank_20190827-%2B_v02.sav.
Downloaded themis_skymap_rank_20200227-%2B_v02.sav.
Downloaded themis_skymap_rank_20200919-%2B_v02.sav.
Downloaded themis_skymap_rank_20210309-%2B_v02.sav.
Downloaded themis_skymap_rank_20210829-%2B_v02.sav.
Downloaded themis_skymap_rank_20220304-%2B_v02.sav.

Satellite footprint#

Now we define an orbit path of a low Earth orbiting satellite (i.e. its footprint). This is a north-south satellite track oriented to the east of the THEMIS/RANK imager. In this context, lla stands for the (latitude, longitude, altitude) coordinates.

[4]:

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)
sat_lla = np.array([lats, lons, alts]).T
sat_time = asi.data.time

Create an asilib.Conjunction() object that handles mapping between the satellite and the imager. It takes in an Imager instance and arrays of the satellite times and LLA coordinates.

[5]:

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

Now, here are two steps that we are ignoring that you’ll likely need to implement:

  1. Map the satellite’s LLA coordinates along the magnetic field line from the satellite altitude down to 110 km (or whatever you chose for the alt kwarg.) This is done via Conjunction.lla_footprint() that requires the IRBEM library. IRBEM can be hard to install; in the future, I plan to change remove IRBEM in favor of geopack (or a similar package).

  2. Normally the satellite LLA time stamps are not the same as the ASI. In that case you will need to call Conjunction.interp_sat() to interpolate the LLA coordinates to the ASI timestamps. Note: this method does not handle interpolation well across the anti-meridian (-180/180 degree longitude). If it detects that you’re interpolating over it, it will issue a warning.

Nearest pixel intensity#

[6]:

sat_azel, sat_azel_pixels = conjunction_obj.map_azel()
nearest_pixel_intensity = conjunction_obj.intensity(box=None)

Mean pixel intensity in a 20x20 km area.#

The mean intensity is calculated with a masked array. It contains np.nan outside the 20x20 km area, and 1s inside.

[7]:

area_intensity = conjunction_obj.intensity(box=(10, 10))

# You don't need to calculate the area mask if you just need the intensity, but this is useful if you
# want to animate and visualize the area
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

The following code block contains many steps.

  1. We create three subplots and initialize the animation generator.

  2. We loop over each image.

    • In the first subplot plot the entire satellite footprint using the sat_azel_pixels array.

    • In the first subplot plot the instantaneous footprint.

    • In the first subplot plot the 20x20 km area contour.

    • In the second and third subplots plot the auroral intensity from a) the nearest pixel, and b) the 20x20 km area.

    • In the second and third subplots plot a vertical line at the current time.

    • Annotate the first subplot.

[8]:

fig, ax = plt.subplots(
    3, 1, figsize=(7, 10), gridspec_kw={'height_ratios': [4, 1, 1]}, constrained_layout=True
)
ax[1].set(ylabel='ASI intensity\nnearest pixel [counts]')
ax[2].set(xlabel='Time', ylabel='ASI intensity\n10x10 km area [counts]')

gen = asi.animate_fisheye_gen(
    ax=ax[0], azel_contours=True, overwrite=True, cardinal_directions='news'
)

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()
        vline2.remove()
        text_obj.remove()
    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',
    )
plt.close()

20170915_023200_023500_themis_rank_fisheye.mp4: |###################### | 97%
Animation saved to C:\Users\shumkms1\asilib-data\animations\20170915_023200_023500_themis_rank_fisheye.mp4

[9]:

# When you run this, you should see the video below in your asilib-data/movies directory.
Video('https://github.com/mshumko/asilib/blob/main/docs/_static/v1/20170915_023200_023500_themis_rank_fisheye.mp4?raw=true')

[9]: