SciQLop / speasy

Space Physics made EASY! A simple Python package to deal with main Space Physics WebServices (CDA,SSC,AMDA,..)
Other
24 stars 7 forks source link

3DView support #80

Open jeandet opened 1 year ago

jeandet commented 1 year ago

@vgenot I've added some ssc comparisons here, does the absolute and relative errors thresholds look OK to you? Not sure why MMS1 fails this test in GSM.

jeandet commented 1 year ago

Here are two plots for MMS1 in GSE and GSM: MMS1_GSE MMS1_GSM

GSE seems pretty OK to me with < 10 km accuracy, while trajectories in GSM coordinates system seems to have a bigger error (up to 200 km). I guess this comes from differences in magnetic pole position determination since Z is the most affected axis and Y a bit less (X being really accurate).

jeandet commented 1 year ago

just for the record the script used to produce plots (using speasy private API):


from datetime import timedelta

import numpy as np
import matplotlib.pyplot as plt

import speasy as spz
from speasy.webservices.cdpp_3dview.ws import _WS_impl

_ws = _WS_impl()

def find_body(name: str):
    bodies = _ws.get_spacecraft_list()
    for body in bodies:
        if body.name == name or body.name.lower() == name.lower():
            return body
    return None

def find_frame(name: str):
    frames = _ws.get_frame_list()
    for frame in frames:
        if frame.name == name or frame.name.lower() == name.lower():
            return frame
    return None

def compare_ssc_trajectory(ssc_index, coordinate_system, cdpp_coordinate_system=None):
    cdpp_coordinate_system = cdpp_coordinate_system or coordinate_system
    start_time = spz.core.make_utc_datetime("2021-01-01")
    stop_time = spz.core.make_utc_datetime("2021-01-05")
    ssc_data = spz.get_data(ssc_index, start_time, stop_time,
                            coordinate_system=coordinate_system)
    ssc_data = ssc_data.to_dataframe()

    cdpp_3DView_data = _ws.get_orbit_data(body=find_body(ssc_index.Id),
                                          start_time=start_time,
                                          stop_time=stop_time,
                                          frame=find_frame(coordinate_system),
                                          time_vector=list(
                                              map(spz.core.make_utc_datetime, ssc_data.index.to_pydatetime())))

    cdpp_3DView_data = cdpp_3DView_data.to_dataframe()

    ssc_data.columns = cdpp_3DView_data.columns
    abs_error = abs(ssc_data - cdpp_3DView_data)
    percent_error = 100 * abs((ssc_data - cdpp_3DView_data) / ssc_data)

    fig, axs = plt.subplots(3, 1, sharex=True)
    axs[0].set_title(f"{ssc_index.Id} Trajectories from SSCWeb and 3DView in {coordinate_system} frame")
    axs[0].plot(ssc_data, label=ssc_data.columns)
    axs[0].plot(cdpp_3DView_data, label=cdpp_3DView_data.columns)
    axs[0].set_ylabel('Spacecraft position (km)')
    axs[0].legend()

    axs[1].set_title(f"Relative error of 3DView trajectory vs SSC")
    axs[1].plot(percent_error, label=percent_error.columns)
    axs[1].set_ylabel('Error (%)')
    axs[1].set_ylim(0, 10)
    axs[1].legend()

    axs[2].set_title(f"Absolute error of 3DView trajectory vs SSC")
    axs[2].plot(abs_error, label=abs_error.columns)
    axs[2].set_ylabel('Error (km)')
    axs[2].legend()

    plt.show()

compare_ssc_trajectory(spz.inventories.tree.ssc.Trajectories.mms1, "GSM")
compare_ssc_trajectory(spz.inventories.tree.ssc.Trajectories.mms1, "GSE")
sonarcloud[bot] commented 2 months ago

Quality Gate Failed Quality Gate failed

Failed conditions
2 Security Hotspots

See analysis details on SonarCloud