In [1]:

import datetime as dt

import holoviews as hv
import hvplot.pandas
import hvplot.xarray
import pandas as pd
import panel as pn
import pyspedas
from pytplot import get_data

pn.extension(sizing_mode="stretch_width")
from loguru import logger
from utils import *

THEMIS Level 2 FGM file info:

Fluxgate Magnetometer (FGM) data: - The Level 2 file includes FGE (engineering) magnetic field, FGH (high-resolution) magnetic field, FGL (low-resolution) magnetic field, and FGS (spin-resolution) magnetic field. - The FGH, FGL, and FGE data are given in Spinning Spacecraft (SSL), Despun Spacecraft (DSL), Geocentric Solar Ecliptic (GSE) and Geocentric Solar Magnetospheric (GSM) coordinates. The FGS data is given in DSL, GSE, and GSM coordinates. - Units are nanotesla.

In [38]:

def calculate_PVI_xr(vec, tau, resample_frequency=None, interval_of_averaging=None):
    """
    This function calculates the Partial Variance of Increments (PVI) series for a given time series.

    Parameters:
    vec (np.array): The input time series, represented as a numpy array.
    tau (int): The time lag, in unit `s`,typically selected to lie in the inertial range of the fluctuations.
    resample_frequency (int): The resample frequency, in unit `s`. If None, defaults to tau.
    interval_of_averaging (int): The number of samples over which to compute the trailing average. It's often chosen to be comparable to, or greater than, a correlation length (or time) of the signal.

    Returns:
    PVI_series (np.array): The resulting PVI series.
    """

    # Sample the vector at the given time lag (tau)
    # vec_sampled = vec.resample(time=tau).mean(dim='time')
    if resample_frequency is None:
        resample_frequency = tau

    # Interpolate to a regular time grid
    # vec_sampled = vec.resample(time=pd.Timedelta(resample_frequency, unit='s')).interpolate('linear')
    vec_sampled = vec.resample(
        time=pd.Timedelta(resample_frequency, unit="s")
    ).interpolate("linear")

    # Calculate the increments in the vector
    n = int(tau / resample_frequency)
    increments = vec_sampled.diff(dim="time", n=n, label="lower")
    # Calculate the magnitudes of these increments
    mag_increments = linalg.norm(increments, dims="v_dim")

    # Square the magnitudes of the increments, and compute a moving average over the specified interval
    if interval_of_averaging is None:
        normalized_factor = np.sqrt(np.mean(np.square(mag_increments)))
        PVI_series = mag_increments / normalized_factor
    else:
        PVI_series = mag_increments
        ...
        # TODO: Implement moving average
    if "units" in PVI_series.attrs:
        del PVI_series.attrs["units"]
    PVI_series["time"] = PVI_series["time"] + pd.Timedelta(tau / 2, unit="s")
    return PVI_series.rename("PVI")


def PVI_map(vec, tau_range, resample_frequency=None):
    """_summary_

    Args:
        vec (_type_): _description_
    """
    if resample_frequency == None:
        PVI_series = xr.concat(
            [calculate_PVI_xr(vec, tau) for tau in tau_range], dim="tau"
        )
    else:
        PVI_series = xr.concat(
            [calculate_PVI_xr(vec, tau, resample_frequency) for tau in tau_range],
            dim="tau",
        )
    return PVI_series

In [2]:

probe = "c"
start = dt.datetime(2019, 1, 7)
end = dt.date.today()
datatype = "fgs"

thx_fgs = f"th{probe}_fgs_gsm"
thx_fgl = f"th{probe}_fgl_gsm"

tstart = start.strftime("%Y-%m-%d")
tstop = (start + dt.timedelta(1)).strftime("%Y-%m-%d")
trange = [tstart, tstop]

pyspedas.themis.fgm(probe=probe, trange=trange, varnames=[thx_fgs, thx_fgl])
thx_fgs_data = get_data(thx_fgs, xarray=True)
# thx_fgl_data = get_data(thx_fgl, xarray=True)

fgl_tstart, fgl_tstop = pytplot.get_timespan(thx_fgl)
fgl_tstart = pd.Timestamp(fgl_tstart, unit="s")
fgl_tstop = pd.Timestamp(fgl_tstop, unit="s")

22-Jul-23 16:17:39: Downloading remote index: http://themis.ssl.berkeley.edu/data/themis/thc/l2/fgm/2019/
22-Jul-23 16:17:39: File is current: /Users/zijin/data/themis/thc/l2/fgm/2019/thc_l2_fgm_20190107_v01.cdf

In [3]:

# Time index may contain duplicate values
thx_fgs_data = thx_fgs_data.groupby("time").first()

# Interpolate to a regular time grid
thx_fgs_data = thx_fgs_data.resample(time=pd.Timedelta(8, unit="s")).interpolate(
    "linear"
)

In [42]:

tau_range = [8, 16, 32, 64, 128, 256]
pvi = PVI_map(thx_fgs_data, tau_range, tau_range[0])

In [43]:

pvi = pvi.assign_coords({"tau": tau_range})

In [44]:

pvi.hvplot.quadmesh(x="time", y="tau", logy=True)

22-Jul-23 16:41:56: /Users/zijin/mambaforge/envs/cool_solar_wind/lib/python3.10/site-packages/bokeh/core/property/bases.py:259: DeprecationWarning: elementwise comparison failed; this will raise an error in the future.
  return new == old

In [2]:

date_slider = pn.widgets.DateSlider(
    name="Date Slider",
    start=start,
    end=end,
    step=1000
    * 60
    * 60
    * 24,  # (number): The selected step i the slider in milliseconds
)

overview_pane = pn.pane.PNG(
    themis_image_url(start, probe, "overview_moms"), name="overview_moms"
)
orbit_pane = pn.pane.GIF(
    themis_image_url(start, probe, "orbit_moon"), name="orbit_moon"
)


def update_pane_image_callback(target, event):
    date = event.new
    image_url = themis_image_url(date, probe, target.name)
    target.object = image_url


date_slider.link(overview_pane, callbacks={"value": update_pane_image_callback})
date_slider.link(orbit_pane, callbacks={"value": update_pane_image_callback})


def index_plot(probe, date):
    thx_fgm = f"th{probe}_fgs_gsm"
    thx_fgl = f"th{probe}_fgl_gsm"

    tstart = date.strftime("%Y-%m-%d")
    tstop = (date + dt.timedelta(1)).strftime("%Y-%m-%d")
    trange = [tstart, tstop]

    pyspedas.themis.fgm(probe=probe, trange=trange, varnames=[thx_fgm, thx_fgl])
    thx_fgm_data = get_data(thx_fgm, xarray=True)

    fgl_tstart, fgl_tstop = pytplot.get_timespan(thx_fgl)
    fgl_tstart = pd.Timestamp(fgl_tstart, unit="s")
    fgl_tstop = pd.Timestamp(fgl_tstop, unit="s")

    pvi = calculate_PVI_xr(thx_fgm_data, "10s")
    mag_change = 10 * calculate_mag_change_xr(thx_fgm_data, "10s")

    pvi_plot = pvi.hvplot(x="time")
    threshold_line = hv.Curve([(fgl_tstart, 5), (fgl_tstop, 5)])
    thx_fgm_plot = thx_fgm_data.hvplot(x="time", by="v_dim")

    # thx_fgl_data = get_data(thx_fgl, xarray=True)
    # thx_fgl_plot = thx_fgl_data.hvplot(x="time", by="v_dim", kind="scatter")

    pvi_plot.opts(ylim=(0, 10), alpha=0.3)
    threshold_line.opts(alpha=0.3)

    index_plot = pvi_plot * mag_change.hvplot(x="time", alpha=0.3) * threshold_line

    markdown = pn.pane.Markdown()

    def update_widget(event):
        temp_tstart = pd.Timestamp(event.new[0])
        temp_tstop = pd.Timestamp(event.new[1])
        temp_datatype = "fgs"

        if temp_tstart > fgl_tstart and temp_tstop < fgl_tstop:
            temp_datatype = "fgl"

        markdown.object = f"""
        ```
        {{"probe": "{probe}", "tstart": "{temp_tstart.strftime('%Y-%m-%dT%H:%M:%S')}", "tstop": "{temp_tstop.strftime('%Y-%m-%dT%H:%M:%S')}", "datatype": "{temp_datatype}"}},
        ```
        """

    rng = hv.streams.RangeX(source=index_plot)
    rng.param.watch(update_widget, "x_range")

    return pn.Column(markdown, pn.Row(index_plot, thx_fgm_plot))


index_pane = hvplot.bind(lambda date: index_plot(probe, date), date_slider)

app = pn.Column(pn.Row(date_slider), index_pane, pn.Row(overview_pane, orbit_pane))
app

In [34]:

date = dt.datetime(2019, 1, 7)
tstart = date.strftime("%Y-%m-%d")
tstop = (date + dt.timedelta(1)).strftime("%Y-%m-%d")
trange = [tstart, tstop]
pyspedas.themis.fgm(probe=probe, trange=trange)

29-Jun-23 12:28:55: Downloading remote index: http://themis.ssl.berkeley.edu/data/themis/thc/l2/fgm/2019/
29-Jun-23 12:28:55: File is current: /Users/zijin/data/themis/thc/l2/fgm/2019/thc_l2_fgm_20190107_v01.cdf

['thc_fgs_btotal',
 'thc_fgs_gse',
 'thc_fgs_gsm',
 'thc_fgs_dsl',
 'thc_fgl_btotal',
 'thc_fgl_gse',
 'thc_fgl_gsm',
 'thc_fgl_dsl',
 'thc_fgl_ssl',
 'thc_fgh_btotal',
 'thc_fgh_gse',
 'thc_fgh_gsm',
 'thc_fgh_dsl',
 'thc_fgh_ssl']

In [3]:

import holoviews as hv

# Widgets
import pandas as pd
import panel as pn

# (index_plot + thx_fgm_plot).cols(1)


xr = thx_fgm_data
start = xr.time.data.min()
end = xr.time.data.max()

start = pd.Timestamp(start).to_pydatetime()
end = pd.Timestamp(end).to_pydatetime()

# Create a DateRangeSlider widget
date_range_slider = pn.widgets.DateRangeSlider(
    name="Date Range Slider",
    start=start,
    end=end,
    value=(start, end),
    format="%H:%M:%S",
)

# Create a DatetimeRangeInput widget
datetime_range_input = pn.widgets.DatetimeRangeInput(
    name="Datetime Range Input",
    value=(start, end),
    start=start,
    end=end,
)

# Create a Markdown widget to copy
markdown = pn.pane.Markdown()


def markdown_callback(target, event):
    target.object = f"""
    ```
    event = {{
        "probe": "{probe}",
        "tstart": "{event.new[0]}",
        "tstop": "{event.new[1]}",
        "datatype": "{datatype}",
    }}
    ```"""


# Link the values of the DateRangeSlider and DatetimeRangeInput widgets
date_range_slider.link(datetime_range_input, value="value")
date_range_slider.link(markdown, callbacks={"value": markdown_callback})
datetime_range_input.link(date_range_slider, value="value")


def plot(dt_range):
    tstart = str(dt_range[0])
    tstop = str(dt_range[1])
    return plot_event(probe, tstart, tstop, datatype, output="all")


interactive_plot = pn.bind(plot, date_range_slider)

pn.Column(
    pn.Row(date_range_slider, datetime_range_input),
    markdown,
    pn.Row(index_plot, thx_fgm_plot.apply.opts(xlim=date_range_slider)),
    interactive_plot,
)

In [54]:

pd.Timestamp(tstart)

Timestamp('2019-01-01 00:00:00')

In [51]:

pd.Timestamp(_tstart, unit="s").strftime("%Y-%m-%dT%H:%M:%S")

'2019-01-01T00:00:30'

In [18]:

import numpy as np

pd.Timestamp("2002-06-28") > pd.Timestamp("2002-06-29")

False