Comparing electron temperatures from EFI IDM and EFI LP products¶

SPACECRAFT = "B" 
START_TIME = '2022-06-30T00:00:00Z'
END_TIME = '2022-07-30T00:00:00Z'

from viresclient import SwarmRequest

SERVER_URL = "https://vires.services/ows"
request = SwarmRequest(SERVER_URL)
collection_efi_idm = f'SW_PREL_EFI{SPACECRAFT}IDM_2_'
collection_efi_lp =f'SW_OPER_EFI{SPACECRAFT}_LP_1B'

start_orbit_number = request.get_orbit_number(SPACECRAFT, START_TIME)
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)

def request_data(url, collection, measurements,start_time, end_time):
    """
    Make a request for data from the Swarm API.
    """
    request = SwarmRequest(url)
    request.set_collection(collection)
    request.set_products(
        measurements=measurements,
        auxiliaries=['OrbitDirection', 'OrbitNumber', 'QDLat', 'QDOrbitDirection'],
    )

    data = request.get_between(
        start_time=start_time,
        end_time=end_time,
        asynchronous=True,
    ).as_xarray()
    return data

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 viresclient._wps.time_util import parse_duration
from scipy.interpolate import interp1d
import numpy as np


def match_times(times1, times2, atol):
    assert times1.ndim == 1 # assert 1D array
    assert times2.ndim == 1 # assert 1D array
    
    epoch = times1[0] if times1.size > 0 else None
    seconds1 = times_to_seconds(times1, epoch=epoch)
    seconds2 = times_to_seconds(times2, epoch=epoch)
    seconds_atol = dtimes_to_seconds(atol)
    
    # find the nearest neigbour
    assert ((seconds2[1:] - seconds2[:-1]) >= 0).all() # assert sorted values
    idx2 = interp1d(
        seconds2, np.arange(seconds2.size),
        kind='nearest',
        copy=False,
        fill_value="extrapolate",
        assume_sorted=True,
    )(seconds1).astype("int64")

    # pick values within the given tolerance
    idx1 = (abs(seconds2[idx2] - seconds1) <= seconds_atol).nonzero()[0]
    idx2 = idx2[idx1]
    
    return idx1, idx2

idx_lp, idx_idm = match_times(
    data_lp["Timestamp"].values,
    data_idm["Timestamp"].values,
    atol=parse_duration("PT0.5S"),
)

te_lp = data_lp["Te"].values[idx_lp]
te_idm = data_idm["T_e"].values[idx_idm]

def _get_lim(values1, values2):
    return (
        min(values1.min(), values2.min()),
        max(values1.max(), values2.max()),
    )

import numpy as np
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