AEBS - AEJxLPS_2F and AEJxPBS_2F products (using 2019 preliminary data)¶

#------------------------------------------------------------------------------
# User input:

SPACECRAFT = 'A'
TIME = '2015-06-02T00:00:00Z'

SERVER_URL = 'https://staging.viresdisc.vires.services/ows'
#SERVER_URL = 'https://testing.vires.services/ows'

#------------------------------------------------------------------------------

from viresclient import SwarmRequest
from viresclient._wps.time_util import parse_datetime

# point types
MASK_EJ_TYPE = PT_EEJ = 0x1
MASK_BOUNDARY_TYPE = PT_POLAR = 0x2
MASK_POINT_TYPE = PT_PEAK = 0x4
PT_WEJ = PT_BOUNDARY = PT_EQUATORIAL = 0x0

SwarmRequest.COLLECTIONS.update({
    'AEJ_LPS': ['SW_OPER_AEJ{}LPS_2F'.format(x) for x in 'ABC'],
    'AEJ_PBS': ['SW_OPER_AEJ{}PBS_2F'.format(x) for x in 'ABC'],
    'AEJ_PBS:GroundMagneticDisturbance': [
        'SW_OPER_AEJ{}PBS_2F:GroundMagneticDisturbance'.format(x) for x in 'ABC'
    ],
})

SwarmRequest.PRODUCT_VARIABLES.update({
    'AEJ_LPS': [
        'Latitude_QD', 'Longitude_QD', 'MLT_QD',
        'J_CF', 'J_DF', 'J_CF_SemiQD', 'J_DF_SemiQD', 'J_C'
    ],
    'AEJ_PBS': [
        'Latitude_QD', 'Longitude_QD', 'MLT_QD',
        'Flags', 'PointType',
    ],
    'AEJ_PBS:GroundMagneticDisturbance': ['B_NEC'],
})

try:
    request = SwarmRequest(SERVER_URL)

    orbit_number = request.get_orbit_number(SPACECRAFT, parse_datetime(TIME))
    start_time, end_time =  request.get_times_for_orbits(SPACECRAFT, orbit_number, orbit_number)
    print('Spacecraft:   %s' % SPACECRAFT)
    print('Orbit number: %s' % orbit_number)
    print('Start time:  ', start_time)
    print('End time:    ', end_time)
    
    # retrieve latitude profiles
    request.set_collection('SW_OPER_AEJ%sLPS_2F' % SPACECRAFT)
    request.set_products(
        measurements=['Latitude_QD', 'Longitude_QD', 'MLT_QD', 'J_CF', 'J_DF', 'J_CF_SemiQD', 'J_DF_SemiQD', 'J_C'],
        auxiliaries=['QDOrbitDirection', 'OrbitNumber'],
    )
    response = request.get_between(
        start_time=start_time,
        end_time=end_time,
        asynchronous=False,
    )
    print('Sources:\n\t%s' % '\n\t'.join(response.sources))
    profile = response.as_xarray()
    
    # retrieve peeks and boundaries
    request.set_collection('SW_OPER_AEJ%sPBS_2F' % SPACECRAFT)
    request.set_products(
        measurements=['Latitude_QD', 'Longitude_QD', 'MLT_QD', 'Flags', 'PointType'],
        auxiliaries=['QDOrbitDirection', 'OrbitNumber'],
    )
    response = request.get_between(
        start_time=start_time,
        end_time=end_time,
        asynchronous=False,
    )
    print('Sources:\n\t%s' % '\n\t'.join(response.sources))
    boundaries = response.as_xarray()

    # retrieve ground magnetic disturbances
    request.set_collection('SW_OPER_AEJ%sPBS_2F:GroundMagneticDisturbance' % SPACECRAFT)
    request.set_products(
        measurements=['B_NEC'],
    )
    response = request.get_between(
        start_time=start_time,
        end_time=end_time,
        asynchronous=False,
    )
    print('Sources:\n\t%s' % '\n\t'.join(response.sources))
    ground_mag_dist = response.as_xarray()
    
    # retrieve spacecraft orbit
    request.set_collection('SW_OPER_MAG%s_LR_1B' % SPACECRAFT)
    request.set_products(
        measurements=[],
        auxiliaries=['QDLat', 'QDLon', 'MLT'],
    )
    response = request.get_between(
        start_time=start_time,
        end_time=end_time,
        asynchronous=False,
    )
    print('Sources:\n\t%s' % '\n\t'.join(response.sources))
    orbit = response.as_xarray()
    
except Exception as error:
    print('ERROR: ', error)
else:
    print(profile)
    print(boundaries)
    print(ground_mag_dist)
    print(orbit)

from numpy import isnan, logical_and, abs as aabs, zeros
from matplotlib.pyplot import figure, subplot, show
from scipy.interpolate import interp1d
%matplotlib inline

def get_quadrant_mask(data, ascending, north):
    orbit_direction = data['QDOrbitDirection'].values
    latitude = data['Latitude_QD'].values
    return logical_and(
        (orbit_direction > 0 if ascending else orbit_direction < 0),
        (latitude > 0 if north else latitude < 0)
    )


def format_time(time):
    return time.replace(microsecond=0).isoformat()


def get_title(ascending, north):
    pass_ = 'Ascending' if ascending else 'Descending'
    hemisphere = 'North' if north else 'South'
    return '%s / %s (%s:%s:%s)' % (hemisphere, pass_, SPACECRAFT, orbit_number, format_time(start_time))


def plot(ax, ascending, north):
    p_mask = get_quadrant_mask(profile, ascending=ascending, north=north)
    b_mask = get_quadrant_mask(boundaries, ascending=ascending, north=north)
    
    b_mask_eb = b_mask & (point_type&(MASK_POINT_TYPE|MASK_BOUNDARY_TYPE) == PT_BOUNDARY|PT_EQUATORIAL)
    b_mask_pb = b_mask & (point_type&(MASK_POINT_TYPE|MASK_BOUNDARY_TYPE) == PT_BOUNDARY|PT_POLAR)
    b_mask_min = b_mask & (point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_WEJ|PT_PEAK)
    b_mask_max = b_mask & (point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_EEJ|PT_PEAK)

    b_mask_wej = b_mask & (point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_WEJ|PT_BOUNDARY)
    b_mask_eej = b_mask & (point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_EEJ|PT_BOUNDARY)
    
    l_wej, = ax.plot(b_qd_latitude[b_mask_wej], zeros(b_qd_latitude[b_mask_wej].shape), color='tab:red')
    l_eej, = ax.plot(b_qd_latitude[b_mask_eej], zeros(b_qd_latitude[b_mask_eej].shape), color='tab:purple')
    
    l_cf, = ax.plot(p_qd_latitude[p_mask], p_qd_j_cf[p_mask], label='J_CF')
    l_df, = ax.plot(p_qd_latitude[p_mask], p_qd_j_df[p_mask], label='J_DF')
    
    l_eb, = ax.plot(b_qd_latitude[b_mask_eb], b_qd_j_df[b_mask_eb], '<' if north==ascending else '>', color='tab:green')
    l_pb, = ax.plot(b_qd_latitude[b_mask_pb], b_qd_j_df[b_mask_pb], '>' if north==ascending else '<', color='tab:orange')
    l_min, = ax.plot(b_qd_latitude[b_mask_min], b_qd_j_df[b_mask_min], 'v', color='tab:red')
    l_max, = ax.plot(b_qd_latitude[b_mask_max], b_qd_j_df[b_mask_max], '^', color='tab:purple')
    
    if not ascending:
        xmin, xmax = ax.get_xlim(); ax.set_xlim((xmax, xmin)) # flip x axis    
    ymax = aabs(ax.get_ylim()).max(); ax.set_ylim((-ymax, ymax)) # center y axis around 0
        
    ax2 = ax.twinx()
    l_r, = ax2.plot(p_qd_latitude[p_mask], p_qd_j_r[p_mask], label='J_R', color='tab:cyan')
    ax2.set_ylabel('J_C / A/km2')
    ymax = aabs(ax2.get_ylim()).max(); ax2.set_ylim((-ymax, ymax))
    
    ax.grid()
    ax.set_title(get_title(ascending=ascending, north=north))
    ax.set_xlabel('QD Latitude / deg')
    ax.set_ylabel('J_SemiQD / A/km')
    
    ax.legend((l_wej, l_eej, l_df, l_cf, l_r, l_min, l_max, l_eb, l_pb), ('WEJ', 'EEJ', 'J_DF', 'J_CF', 'J_R', 'MIN', 'MAX', 'EB', 'PB' ))
 
    for x, y in zip(b_qd_latitude[b_mask_min], b_qd_j_df[b_mask_min]):
        ax.text(x, y, ' %.3g A/km' % y, color='tab:red', va='top', ha='left')
    for x, y in zip(b_qd_latitude[b_mask_max], b_qd_j_df[b_mask_max]):
        ax.text(x, y, ' %.3g A/km' % y, color='tab:purple', va='bottom', ha='left')


p_time = profile['Timestamp'].values
p_qd_latitude = profile['Latitude_QD'].values
p_qd_j_cf = profile['J_CF_SemiQD'].values
p_qd_j_df = profile['J_DF_SemiQD'].values
p_qd_j_r = profile['J_C'].values

b_time = boundaries['Timestamp'].values
b_qd_latitude = boundaries['Latitude_QD'].values
#b_qd_j_df = boundaries['J_DF_SemiQD'].values
b_qd_j_df = interp1d(p_time.astype('int64'), p_qd_j_df, kind="nearest")(b_time.astype('int64'))


p_qd_j_df
b_qd_j_df[isnan(b_qd_j_df)] = 0
point_type = boundaries['PointType'].values

fig = figure(figsize=(18, 15), dpi=100)

plot(subplot(221), ascending=True, north=True)
plot(subplot(222), ascending=False, north=True)
plot(subplot(223), ascending=False, north=False)
plot(subplot(224), ascending=True, north=False)

from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

p_time = profile['Timestamp'].values
p_qd_j_cf = profile['J_CF_SemiQD'].values
p_qd_j_df = profile['J_DF_SemiQD'].values
p_qd_j_r = profile['J_C'].values

b_time = boundaries['Timestamp'].values
#b_qd_j_df = boundaries['J_DF_SemiQD'].values
b_qd_j_df = interp1d(p_time.astype('int64'), p_qd_j_df, kind="nearest")(b_time.astype('int64'))
b_qd_j_df[isnan(b_qd_j_df)] = 0
point_type = boundaries['PointType'].values

def plot(ax):
    
    def _plot(ascending, north):
        p_mask = get_quadrant_mask(profile, ascending=ascending, north=north)
        l_cf, = ax.plot(p_time[p_mask], p_qd_j_cf[p_mask], '-', markersize=4, label='J_CF', color='tab:blue')
        l_df, = ax.plot(p_time[p_mask], p_qd_j_df[p_mask], '-', markersize=4, label='J_DF', color='tab:orange')
        return l_cf, l_df
    
    l_cf, l_df = _plot(ascending=True, north=True)
    l_cf, l_df = _plot(ascending=False, north=True)
    l_cf, l_df = _plot(ascending=False, north=False)
    l_cf, l_df = _plot(ascending=True, north=False)
    
    b_mask_eb = point_type&(MASK_POINT_TYPE|MASK_BOUNDARY_TYPE) == PT_BOUNDARY|PT_EQUATORIAL
    b_mask_pb = point_type&(MASK_POINT_TYPE|MASK_BOUNDARY_TYPE) == PT_BOUNDARY|PT_POLAR
    b_mask_min = point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_WEJ|PT_PEAK
    b_mask_max = point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_EEJ|PT_PEAK

    b_mask_wej = point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_WEJ|PT_BOUNDARY
    b_mask_eej = point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_EEJ|PT_BOUNDARY

    ymax = 1.2 * max(aabs(p_qd_j_cf).max(), aabs(p_qd_j_df).max())
    for x in b_time[b_mask_min]:
        ax.plot([x, x], [-ymax, ymax], '-', color='tab:red')
    for x in b_time[b_mask_max]:
        ax.plot([x, x], [-ymax, ymax], '-', color='tab:purple')
    
    for x0, x1 in zip(b_time[b_mask_wej][:-1:2], b_time[b_mask_wej][1::2]):
        l_wej, = ax.plot([x0, x1], [0, 0], color='tab:red')
    for x0, x1 in zip(b_time[b_mask_eej][:-1:2], b_time[b_mask_eej][1::2]):
        l_eej, = ax.plot([x0, x1], [0, 0], color='tab:purple')
    
    l_eb, = ax.plot(b_time[b_mask_eb], b_qd_j_df[b_mask_eb]*0, '<', color='tab:green')
    l_pb, = ax.plot(b_time[b_mask_pb], b_qd_j_df[b_mask_pb]*0, '>', color='tab:orange')
    l_min, = ax.plot(b_time[b_mask_min], b_qd_j_df[b_mask_min], 'v', color='tab:red')
    l_max, = ax.plot(b_time[b_mask_max], b_qd_j_df[b_mask_max], '^', color='tab:purple')
  
    ax.grid()
    
    ax.set_xlabel('time')
    ax.set_ylabel('J_SemiQD / A/km')

    for x, y in zip(b_time[b_mask_min], b_qd_j_df[b_mask_min]):
        ax.text(x, y, ' %.3g A/km' % y, color='tab:red', va='top', ha='left')
    for x, y in zip(b_time[b_mask_max], b_qd_j_df[b_mask_max]):
        ax.text(x, y, ' %.3g A/km' % y, color='tab:purple', va='bottom', ha='left')
    
    #ymax = aabs(ax.get_ylim()).max();
    ax.set_ylim((-ymax, ymax)) # center y axis around 0
    
    ax = ax.twinx()
    
    def _plot(ascending, north):
        p_mask = get_quadrant_mask(profile, ascending=ascending, north=north)
        l_r, = ax.plot(p_time[p_mask], p_qd_j_r[p_mask], '-', markersize=4, label='J_R', color='tab:cyan')
        return l_r
    
    l_r = _plot(ascending=True, north=True)
    l_r = _plot(ascending=False, north=True)
    l_r = _plot(ascending=False, north=False)
    l_r = _plot(ascending=True, north=False)

    ax.set_ylabel('J_C / A/km2')
    ax.legend((l_wej, l_eej, l_df, l_cf, l_r, l_min, l_max, l_eb, l_pb), ('WEJ', 'EEJ', 'J_DF', 'J_CF', 'J_R', 'MIN', 'MAX', 'EB', 'P' ))

    
    #if not ascending:
    #    xmin, xmax = ax.get_xlim(); ax.set_xlim((xmax, xmin)) # flip x axis
    ymax = aabs(ax.get_ylim()).max(); ax.set_ylim((-ymax, ymax)) # center y axis around 0

fig = figure(figsize=(36, 6), dpi=100)
plot(subplot(1, 1, 1))

from numpy import isnan, logical_and, abs as aabs, cumsum, stack, pi, concatenate, timedelta64
from matplotlib.pyplot import figure, subplot, show
from scipy.interpolate import interp1d
from cartopy.feature import LAND, OCEAN, COASTLINE
from cartopy.crs import Mollweide, Orthographic, PlateCarree
%matplotlib inline

time_threshold = timedelta64(1, 's')

# latitude profile
p_time = profile['Timestamp'].values
p_lat = profile['Latitude'].values
p_lon = profile['Longitude'].values
p_lat_qd = profile['Latitude_QD'].values
p_lon_qd = profile['Longitude_QD'].values
p_mlt = profile['MLT_QD'].values

# boubdaries
b_time = boundaries['Timestamp'].values
b_lat = boundaries['Latitude'].values
b_lon = boundaries['Longitude'].values
b_lat_qd = boundaries['Latitude_QD'].values
b_lon_qd = boundaries['Longitude_QD'].values
b_mlt = boundaries['MLT_QD'].values
point_type = boundaries['PointType'].values

b_mask_eb = point_type&(MASK_POINT_TYPE|MASK_BOUNDARY_TYPE) == PT_BOUNDARY|PT_EQUATORIAL
b_mask_pb = point_type&(MASK_POINT_TYPE|MASK_BOUNDARY_TYPE) == PT_BOUNDARY|PT_POLAR
b_mask_peak_eej = point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_EEJ|PT_PEAK
b_mask_peak_wej = point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_WEJ|PT_PEAK
b_mask_eej = point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_EEJ|PT_BOUNDARY
b_mask_wej = point_type&(MASK_EJ_TYPE|MASK_POINT_TYPE) == PT_WEJ|PT_BOUNDARY
b_eej_pair_start = zeros(point_type[b_mask_eej].shape, 'bool')
b_eej_pair_start[::2] = True
b_wej_pair_start = zeros(point_type[b_mask_wej].shape, 'bool')
b_wej_pair_start[::2] = True

# orbit from MAGx_LR products
o_time = orbit['Timestamp'].values
o_lat = orbit['Latitude'].values
o_lon = orbit['Longitude'].values
o_lat_qd = orbit['QDLat'].values
o_lon_qd = orbit['QDLon'].values
o_mlt = orbit['MLT'].values

# ground magnetic disturbance
g_lat = ground_mag_dist['Latitude'].values
g_lon = ground_mag_dist['Longitude'].values

# EJ mask -> True if point is inside WEJ or EEJ
o_mask_eej = interp1d(
    b_time[b_mask_eej].astype('int64'), b_eej_pair_start, kind='zero', bounds_error=False, fill_value=0.0
)(o_time.astype('int64')).astype('bool')
o_mask_wej = interp1d(
    b_time[b_mask_wej].astype('int64'), b_wej_pair_start, kind='zero', bounds_error=False, fill_value=0.0
)(o_time.astype('int64')).astype('bool')

p_mask_eej = interp1d(
    b_time[b_mask_eej].astype('int64'), b_eej_pair_start, kind='zero', bounds_error=False, fill_value=0.0
)(p_time.astype('int64')).astype('bool')
p_mask_wej = interp1d(
    b_time[b_mask_wej].astype('int64'), b_wej_pair_start, kind='zero', bounds_error=False, fill_value=0.0
)(p_time.astype('int64')).astype('bool')


def split_array(data, threshold, start=None, end=None):
    if start is None:
        start = 0
    if end is None:
        end = len(data)
    split = concatenate((
        [start],
        (abs(data[start+1:end] - data[start:end-1]) > threshold).nonzero()[0] + 1 + start,
        [end]
    ))
    return zip(split[:-1], split[1:])

# -----------------------------------------------------------------------------

def plot_qdlat_vs_time(ax, ylim, label): 
    is_north = min(ylim) > 0
    l_orb, = ax.plot(o_time, o_lat_qd, '-', color='silver')      
    for start, end in split_array(p_time, time_threshold):
        l_lp, = ax.plot(p_time[start:end], p_lat_qd[start:end], '-', color='tab:blue')
    for start, end in split_array(p_time[p_mask_eej], time_threshold):
        l_lp_eej, = ax.plot(p_time[p_mask_eej][start:end], p_lat_qd[p_mask_eej][start:end], '-', color='purple')
    for start, end in split_array(p_time[p_mask_wej], time_threshold):
        l_lp_wej, = ax.plot(p_time[p_mask_wej][start:end], p_lat_qd[p_mask_wej][start:end], '-', color='red')
    l_peak_eej, = ax.plot(b_time[b_mask_peak_eej], b_lat_qd[b_mask_peak_eej], '^', color='purple')
    l_peak_wej, = ax.plot(b_time[b_mask_peak_wej], b_lat_qd[b_mask_peak_wej], 'v', color='red')
    l_eb, = ax.plot(b_time[b_mask_eb], b_lat_qd[b_mask_eb], '<', color='tab:green')
    l_pb, = ax.plot(b_time[b_mask_pb], b_lat_qd[b_mask_pb], '>', color='tab:orange')
    
    ax.set_ylim(ylim)
    ax.grid()
    ax.set_title('AEJxLPS/AEJxPBS - QD Latitudes - %s' % label)
    ax.set_ylabel('QD Latitude / deg')
    ax.set_xlabel('time')
    ax.legend(
        (l_peak_wej, l_peak_eej, l_eb, l_pb, l_lp, l_lp_wej, l_lp_eej, l_orb),
        ('WEJ/MIN', 'EEJ/MAX', 'EB', 'PB', 'LP', 'WEJ', 'EEJ', 'MAG')
    )

fig = figure(figsize=(18, 12), dpi=100)
plot_qdlat_vs_time(subplot(211), (45, 90), 'North')
plot_qdlat_vs_time(subplot(212), (-90, -45), 'South')

# -----------------------------------------------------------------------------

def plot_qdlon_vs_time(ax):
    ylim = (-180, 180)
    
    for start, end in split_array(o_lon_qd, 180.0):
        l_orb, = ax.plot(o_time[start:end], o_lon_qd[start:end], '-', color='silver')
    for start0, end0 in split_array(p_time, time_threshold):
        for start, end in split_array(p_lon_qd, 180.0, start0, end0):
            l_lp, = ax.plot(p_time[start:end], p_lon_qd[start:end], '-', color='tab:blue')    
    for start0, end0 in split_array(p_time[p_mask_eej], time_threshold):
        for start, end in split_array(p_lon_qd[p_mask_eej], 180.0, start0, end0):
            l_lp_eej, = ax.plot(p_time[p_mask_eej][start:end], p_lon_qd[p_mask_eej][start:end], '-', color='purple')
    for start0, end0 in split_array(p_time[p_mask_wej], time_threshold):
        for start, end in split_array(p_lon_qd[p_mask_wej], 180.0, start0, end0):
            l_lp_wej, = ax.plot(p_time[p_mask_wej][start:end], p_lon_qd[p_mask_wej][start:end], '-', color='red')
    
    l_peak_eej, = ax.plot(b_time[b_mask_peak_eej], b_lon_qd[b_mask_peak_eej], '^', color='purple')
    l_peak_wej, = ax.plot(b_time[b_mask_peak_wej], b_lon_qd[b_mask_peak_wej], 'v', color='red')
    l_eb, = ax.plot(b_time[b_mask_eb], b_lon_qd[b_mask_eb], '<', color='tab:green')
    l_pb, = ax.plot(b_time[b_mask_pb], b_lon_qd[b_mask_pb], '>', color='tab:orange')
    
    ax.set_ylim(ylim)
    ax.set_yticks([-180, -135, -90, -45, 0, 45, 90, 135, 180])
    ax.grid()
    ax.set_title('AEJxLPL/AEJxPBL - QD Longitude')
    ax.set_ylabel('QD Longitude / deg')
    ax.set_xlabel('time')
    ax.legend(
        (l_peak_wej, l_peak_eej, l_eb, l_pb, l_lp, l_lp_wej, l_lp_eej, l_orb),
        ('WEJ/MIN', 'EEJ/MAX', 'EB', 'PB', 'LP', 'WEJ', 'EEJ', 'MAG')
    )

fig = figure(figsize=(18, 12), dpi=100)
plot_qdlon_vs_time(subplot(111))

# -----------------------------------------------------------------------------

def plot_mlt_vs_time(ax):
    ylim = (0, 24)
    
    for start, end in split_array(o_mlt, 12.0):
        l_orb, = ax.plot(o_time[start:end], o_mlt[start:end], '-', color='silver')
    for start0, end0 in split_array(p_time, time_threshold):
        for start, end in split_array(p_mlt, 12.0, start0, end0):
            l_lp, = ax.plot(p_time[start:end], p_mlt[start:end], '-', color='tab:blue')    
    for start0, end0 in split_array(p_time[p_mask_eej], time_threshold):
        for start, end in split_array(p_mlt[p_mask_eej], 12.0, start0, end0):
            l_lp_eej, = ax.plot(p_time[p_mask_eej][start:end], p_mlt[p_mask_eej][start:end], '-', color='purple')
    for start0, end0 in split_array(p_time[p_mask_wej], time_threshold):
        for start, end in split_array(p_mlt[p_mask_wej], 12.0, start0, end0):
            l_lp_wej, = ax.plot(p_time[p_mask_wej][start:end], p_mlt[p_mask_wej][start:end], '-', color='red')
    
    l_peak_eej, = ax.plot(b_time[b_mask_peak_eej], b_mlt[b_mask_peak_eej], '^', color='purple')
    l_peak_wej, = ax.plot(b_time[b_mask_peak_wej], b_mlt[b_mask_peak_wej], 'v', color='red')
    l_eb, = ax.plot(b_time[b_mask_eb], b_mlt[b_mask_eb], '<', color='tab:green')
    l_pb, = ax.plot(b_time[b_mask_pb], b_mlt[b_mask_pb], '>', color='tab:orange')
    
    ax.set_ylim(ylim)
    ax.set_yticks([0, 3, 6, 9, 12, 15, 18, 21, 24])
    ax.grid()
    ax.set_title('AEJxLPS/AEJxPBS - Magnetic Local Time')
    ax.set_ylabel('MLT / hour')
    ax.set_xlabel('time')
    ax.legend(
        (l_peak_wej, l_peak_eej, l_eb, l_pb, l_lp, l_lp_wej, l_lp_eej, l_orb),
        ('WEJ/MIN', 'EEJ/MAX', 'EB', 'PB', 'LP', 'WEJ', 'EEJ', 'MAG')
    )

fig = figure(figsize=(18, 12), dpi=100)
plot_mlt_vs_time(subplot(111))

# -----------------------------------------------------------------------------

def plot_maps(ax, north):
    h = 1 if north else -1
    gl = ax.gridlines(crs=PlateCarree(), draw_labels=False, linewidth=0, color='silver', alpha=0.5, linestyle='--')
    ax.set_xlim([-6378073.21863, 6378073.21863])
    ax.set_ylim([-6378073.21863, 6378073.21863])
    ax.add_feature(LAND, facecolor=(1.0, 1.0, 0.9))
    ax.add_feature(OCEAN, facecolor=(0.9, 1.0, 1.0))
    ax.add_feature(COASTLINE, edgecolor='silver')

    #ax.plot(g_lon, g_lat, 'o', color='tab:brown', markersize=20, transform=PlateCarree(), mfc='none')
    ax.plot(g_lon, g_lat, 'x', color='tab:brown', markersize=14, transform=PlateCarree(), mfc='none')
    ax.plot(o_lon[o_lat*h > 0], o_lat[o_lat*h > 0], '-', color='silver', transform=PlateCarree())
    for start, end in split_array(p_time, time_threshold):
        ax.plot(p_lon[start:end], p_lat[start:end], '-', color='tab:blue', transform=PlateCarree())
    for start, end in split_array(p_time[p_mask_eej], time_threshold):
        ax.plot(p_lon[p_mask_eej][start:end], p_lat[p_mask_eej][start:end], '-', color='purple', transform=PlateCarree())
    for start, end in split_array(p_time[p_mask_wej], time_threshold):
        ax.plot(p_lon[p_mask_wej][start:end], p_lat[p_mask_wej][start:end], '-', color='red', transform=PlateCarree())
    ax.plot(b_lon[b_mask_peak_eej], b_lat[b_mask_peak_eej], '^', color='purple', transform=PlateCarree())
    ax.plot(b_lon[b_mask_peak_wej], b_lat[b_mask_peak_wej], 'v', color='red', transform=PlateCarree())
    ax.plot(b_lon[b_mask_eb], b_lat[b_mask_eb], '<', color='tab:green', transform=PlateCarree())
    ax.plot(b_lon[b_mask_pb], b_lat[b_mask_pb], '>', color='tab:orange', transform=PlateCarree())

fig = figure(figsize=(18, 9), dpi=100)
plot_maps(subplot(1, 2, 1, projection=Orthographic(0, 90)), north=True)
plot_maps(subplot(1, 2, 2, projection=Orthographic(0, -90)), north=False)

# -----------------------------------------------------------------------------

def plot_mlt(ax, north=True):
    ax.set_theta_zero_location('S')
    h = 1 if north else -1
    
    def _plot(x, y, *args, **kwargs):
        return ax.plot(x*(pi/12), 90 - y*h, *args, **kwargs)
    
    o_mask = o_lat_qd*h > 0
    for start, end in split_array(o_time[o_mask], time_threshold):
        _plot(o_mlt[o_mask][start:end], o_lat_qd[o_mask][start:end], '-', color='silver', markersize=0.15)
    for start, end in split_array(p_time, time_threshold):
        _plot(p_mlt[start:end], p_lat_qd[start:end], '-', color='tab:blue')
    for start, end in split_array(p_time[p_mask_eej], time_threshold):
        _plot(p_mlt[p_mask_eej][start:end], p_lat_qd[p_mask_eej][start:end], '-', color='purple')
    for start, end in split_array(p_time[p_mask_wej], time_threshold):
        _plot(p_mlt[p_mask_wej][start:end], p_lat_qd[p_mask_wej][start:end], '-', color='red')
    _plot(b_mlt[b_mask_peak_eej], b_lat_qd[b_mask_peak_eej], '^', color='purple')
    _plot(b_mlt[b_mask_peak_wej], b_lat_qd[b_mask_peak_wej], 'v', color='red')
    _plot(b_mlt[b_mask_eb], b_lat_qd[b_mask_eb], '<', color='tab:green')
    _plot(b_mlt[b_mask_pb], b_lat_qd[b_mask_pb], '>', color='tab:orange')
    
    ax.set_ylim(0, 50)
    ax.set_xticklabels(['%2.2i' % (x*12/pi) for x in ax.get_xticks()])
    ax.set_yticklabels(['%.0f' % ((90 - y)*h) for y in ax.get_yticks()])

fig = figure(figsize=(18, 9), dpi=100)
plot_mlt(subplot(1, 2, 1, projection='polar'), north=True)
plot_mlt(subplot(1, 2, 2, projection='polar'), north=False)