Swarm LP Ion Drift and Effective Mass Products in VirES¶

SPACECRAFT = "A" 
START_TIME = "2015-12-15T00:00:00Z"
END_TIME = "2016-01-30T00:00:00Z"
V_MAX = 25

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

SERVER_URL = "https://vires.services/ows"
COLLECTION = f"SW_PREL_EFI{SPACECRAFT}IDM_2_"

request = SwarmRequest(SERVER_URL)

start_orbit_number = request.get_orbit_number(SPACECRAFT, START_TIME, mission = 'Swarm')
end_orbit_number = request.get_orbit_number(SPACECRAFT, END_TIME)
start_time, end_time =  request.get_times_for_orbits(start_orbit_number, end_orbit_number, mission = 'Swarm', spacecraft=SPACECRAFT)
start_time, end_time = request.get_times_for_orbits(start_orbit_number, end_orbit_number, spacecraft=SPACECRAFT,  mission = 'Swarm')

import numpy as np

def get_orbit_latitude(latitude, orbit_direction):
    """ Unwrap latitude along the orbit.
    
    Returns:
        orbit_latitude ... unwrapped latitude values
    Args:
        latitude ... latitude values in degrees between -90 and 90.
        orbit_direction ... orbit direction flags
            +1 ascending direction
            -1 descnding direction
            0 means orbit direction is not available
    """
    orbit_latitude = np.array(latitude, copy=True)
    mask = (np.asarray(orbit_direction) == -1)
    orbit_latitude[mask] = 180 - orbit_latitude[mask] 
    orbit_latitude[orbit_latitude >= 180] -= 360 
    return orbit_latitude

import numpy as np

EPOCH = np.datetime64("2000-01-01")


def times_to_seconds(times, epoch=None, precision="ns", scale=1e-9):
    """ Convert times (datetime64) to seconds (float).
    
    The epoch (zero value) is set by default to 2000-01-01T00:00:00.
    """
    if epoch is None:
        epoch = EPOCH
    return dtimes_to_seconds(np.asarray(times, "datetime64[ns]") - epoch, precision, scale)


def seconds_to_times(seconds, epoch=None, precision="ns", scale=1e9):
    """ Convert seconds (float) to times (datetime64).
    
    The epoch (zero value) is set by default to 2000-01-01T00:00:00.
    """
    if epoch is None:
        epoch = EPOCH
    return seconds_to_dtimes(seconds, precision, scale) + epoch


def dtimes_to_seconds(dtimes, precision="ns", scale=1e-9):
    """ Convert time differences (timedelta64) to seconds (float)
    with nanosecond precision.
    """
    return scale * np.asarray(dtimes, f"timedelta64[{precision}]").astype("int64")


def seconds_to_dtimes(seconds, precision="ns", scale=1e9):
    """ Convert seconds to time differences (timedelta64)
    with nanosecond precision.
    """
    return (scale * np.asarray(seconds)).astype("int64").astype(f"timedelta64[{precision}]")

from scipy.stats import binned_statistic_2d

def bin_data_2d(x, y, v, x_step, y_step, x_extent=None, y_extent=None, statistic='mean'):
    """ Conventience wrapper around the scipy.stats.binned_statistic_2d.
    
    See also https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binned_statistic_2d.html    
    
    Returns:
        result ... the calculated binned statictics (Ny x Nx array)
        x_edge ... edges of the bins in the x coordinate (Nx + 1 array)
        y_edge ... edges of the bins in the y coordinate (Ny + 1 array)
        
    Args:
        x ... 1st coordinate values
        y ... 2nd coordinate values
        v ... binned values
        x_step ... x proposed binning step (adjusted to fit the x_extent)
        y_step ... y proposed binning step (adjusted to fit the y_extent)
        x_extent ... optional pair of min. and max. boulds for the x coordinate
        y_extent ... optional pait of min. and max. boulds for the y coordinate
        staticstics ... type of the statistic to be calculated. For the possible
            options see the scipy.stats.binned_statistic_2d documentation.
    """
    def _get_bins(start, end, step):
        nstep = abs(int(np.ceil((end - start)/step))) + 1
        return np.linspace(start, end, nstep)
    
    if x_extent is None:
        x_extent = [x.min(), x.max()]

    if y_extent is None:
        y_extent = [y.min(), y.max()]

    result, x_edges, y_edges, _ = binned_statistic_2d(
        x, y, v, statistic=statistic, bins=(
            _get_bins(x_extent[0], x_extent[1], x_step),
            _get_bins(y_extent[0], y_extent[1], y_step),
        ),
    )
    return result.T, x_edges, y_edges


def bin_time_data_2d(t, y, v, t_step, y_step, t_extent=None, y_extent=None, statistic='mean'):
    """ Conventience wrapper around the scipy.stats.binned_statistic_2d with time
    as the first coordinate.
    
    See also https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binned_statistic_2d.html    
    
    Returns:
        result ... the calculated binned statictics (Ny x Nt array)
        t_edge ... edges of the bins in the t coordinate (Nt + 1 array)
        y_edge ... edges of the bins in the y coordinate (Ny + 1 array)
        
    Args:
        t ... 1st coordinate values
        y ... 2nd coordinate values
        v ... binned values
        t_step ... x proposed binning step (adjusted to fit the t_extent)
        y_step ... y proposed binning step (adjusted to fit the y_extent)
        t_extent ... optional pair of min. and max. boulds for the t coordinate
        y_extent ... optional pait of min. and max. boulds for the y coordinate
        staticstics ... type of the statistic to be calculated. For the possible
            options see the scipy.stats.binned_statistic_2d documentation.
    """    
    if t_extent is None:
        t_extent = [t.min(), t.max()]
    
    result, t_edges, y_edges = bin_data_2d(
        times_to_seconds(t), y, v,
        x_step=dtimes_to_seconds(t_step),
        y_step=y_step,
        x_extent=times_to_seconds(t_extent),
        y_extent=y_extent,
        statistic=statistic,
    )
    return result, seconds_to_times(t_edges), y_edges

import numpy as np

mask = data['M_i_eff_Flags'].values & 0x01 == 0x00 # extract values for flag bit-0 set to 0
m_i_eff = data["M_i_eff"].values
m_tbt_model = data["M_i_eff_tbt_model"].values
timestamp = data["Timestamp"].values
orbit_number = data["OrbitNumber"].values
latitude = data["Latitude"].values
orbit_direction = data["OrbitDirection"].values
qd_latitude = data["QDLat"].values
qd_orbit_direction = data["QDOrbitDirection"].values


def _get_data_range(values):
    """ Get minimum and maximum values of an array. """
    values = values[~np.isnan(values)] # strip possible NaN values
    return values.min(), values.max()

from matplotlib.colors import ListedColormap

orbit_latitude = get_orbit_latitude(qd_latitude, qd_orbit_direction) % 360

def plot_by_orbit(ax,v_max,orbit_number,orbit_latitude,color_data, mask, label):
    v_min=0
    # v_max = v_max

    values_bin_median, orbit_number_bin, orbit_latitude_bin = bin_data_2d(
        x=orbit_number[mask],
        y=orbit_latitude[mask],
        v=color_data[mask],
        x_step=2, # orbits per bin
        y_step=1, # 1 degree y-bins
        x_extent=[ 
            orbit_number.min(),
            orbit_number.max() + 1,
        ],
        y_extent=[0, 360],
        statistic='median',
    )

    values_bin_count, _, _ = bin_data_2d(
        x=orbit_number,
        y=orbit_latitude,
        v=np.ones(timestamp.shape),
        x_step=2, # orbits per bin
        y_step=1, # 1 degree y-bins
        x_extent=[ 
            orbit_number.min(),
            orbit_number.max() + 1,
        ],
        y_extent=[0, 360],
        statistic='count',
    )


    # remap orbits to times 
    time_bin = seconds_to_times(np.linspace(
        times_to_seconds(timestamp.min()),
        times_to_seconds(timestamp.max()),
        orbit_number_bin.size,
    ))

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

    intensity = 0.85
    data_occurance_cmap = ListedColormap(np.asarray([
        (0.0, 0.0, 0.0, 0.0),
        (intensity, intensity, intensity, 1.0),
    ]))



    # plot rejected data with a light gray color
    ax.pcolormesh(
        time_bin,
        orbit_latitude_bin,
        values_bin_count > 0,
        cmap=data_occurance_cmap,
    )

    # plot the actual values
    im = ax.pcolormesh(
        time_bin,
        orbit_latitude_bin,
        values_bin_median,
        cmap='viridis_r',
        vmin=v_min,
        vmax=v_max
    )

    ax.set_xlabel('Timestamp')
    ax.set_ylabel('Latitude')
    cbar = plt.colorbar(im,ax=ax)
    cbar.set_label(f'{label} / a.m.u.')
    ax.set_ylim(0, 360)
    ax.set_yticks(
        [0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360],
        ["0°↑", "+30°↑", "+60°↑", "+90°", "+60°↓", "+30°↓", "0°↓", "−30°↓", "−60°↓", "−90°", "−60°↑", "−30°↑", "0°↑"]
    )
    ax.grid()