SP3 clock normalisation & Helmert transformation

gnssanalysis/gn_diffaux.py

@@ -1,6 +1,6 @@
 import logging as _logging
 from pathlib import Path as _Path
-from typing import Union
+from typing import Literal, Union
 
 import numpy as _np
 import pandas as _pd
@@ -324,8 +324,8 @@ def compare_clk(
     :return _pd.DataFrame: clk differences in the same units as input clk dfs (usually seconds)
     """
 
-    clk_a = _gn_io.clk.get_AS_entries(clk_a)
-    clk_b = _gn_io.clk.get_AS_entries(clk_b)
+    clk_a = _gn_io.clk.get_sv_clocks(clk_a)
+    clk_b = _gn_io.clk.get_sv_clocks(clk_b)
 
     if not isinstance(norm_types, list):  # need list for 'sv' to be correctly converted to array of SVs to use for norm
         norm_types = list(norm_types)
@@ -708,32 +708,44 @@ def format_index(
 def sp3_difference(
     base_sp3_file: _Path,
     test_sp3_file: _Path,
+    svs: list[str],
+    orb_hlm_mode: Literal[None, "ECF", "ECI"] = None,
+    clk_norm_types: list = [],
 ) -> _pd.DataFrame:
     """
     Compare two SP3 files to calculate orbit and clock differences. The orbit differences will be represented
     in both X/Y/Z ECEF frame and R/A/C orbit frame, and the clock differences will NOT be normalised.
 
     :param _Path base_sp3_file: Path of the baseline SP3 file
     :param _Path test_sp3_file: Path of the test SP3 file
+    :param list[str] svs: List of satellites to process
+    :param str orb_hlm_mode: Helmert transformation to apply to orbits. Can be None, "ECF", or "ECI".
+    :param list clk_norm_types: Normalizations to apply to clocks. Available options include 'epoch', 'daily', 'sv',
+            any satellite PRN, or any combination of them, defaults to empty list
     :return _pd.DataFrame: The Pandas DataFrame containing orbit and clock differences
     """
     base_sp3_df = _gn_io.sp3.read_sp3(str(base_sp3_file))
+    mask = base_sp3_df.index.get_level_values('PRN').isin(svs)
+    base_sp3_df = base_sp3_df[mask]
+
     test_sp3_df = _gn_io.sp3.read_sp3(str(test_sp3_file))
+    mask = test_sp3_df.index.get_level_values('PRN').isin(svs)
+    test_sp3_df = test_sp3_df[mask]
 
     common_indices = base_sp3_df.index.intersection(test_sp3_df.index)
     diff_est_df = test_sp3_df.loc[common_indices, "EST"] - base_sp3_df.loc[common_indices, "EST"]
-
-    diff_clk_df = diff_est_df["CLK"].to_frame(name="CLK") * 1e3  # TODO: normalise clocks
     diff_xyz_df = diff_est_df.drop(columns=["CLK"]) * 1e3
-    diff_rac_df = _gn_io.sp3.diff_sp3_rac(base_sp3_df, test_sp3_df, hlm_mode=None)  # TODO: hlm_mode
 
+    diff_rac_df = _gn_io.sp3.diff_sp3_rac(base_sp3_df, test_sp3_df, hlm_mode=orb_hlm_mode)  # TODO Eugene: compare orbits by constellation
     diff_rac_df.columns = diff_rac_df.columns.droplevel(0)
-
     diff_rac_df = diff_rac_df * 1e3
 
+    diff_clk_df = compare_clk(test_sp3_df, base_sp3_df, norm_types=clk_norm_types)  # TODO Eugene: compare clocks by constellation
+    diff_clk_df = diff_clk_df.stack(dropna=False).to_frame(name="Clock") * 1e3
+
     diff_sp3_df = diff_xyz_df.join(diff_rac_df)
     diff_sp3_df["3D-Total"] = diff_xyz_df.pow(2).sum(axis=1, min_count=3).pow(0.5)
-    diff_sp3_df["Clock"] = diff_clk_df
+    diff_sp3_df = diff_sp3_df.join(diff_clk_df)
     diff_sp3_df["|Clock|"] = diff_clk_df.abs()
 
     format_index(diff_sp3_df)
@@ -744,23 +756,29 @@ def sp3_difference(
 def clk_difference(
     base_clk_file: _Path,
     test_clk_file: _Path,
-    norm_types: list = [],
+    svs: list[str],
+    clk_norm_types: list = [],
 ) -> _pd.DataFrame:
     """
     Compare two CLK files to calculate clock differences with common mode removed (if specified)
     based on the chosen normalisations.
 
     :param _Path base_clk_file: Path of the baseline CLK file
     :param _Path test_clk_file: Path of the test CLK file
-    :param norm_types list: Normalizations to apply. Available options include 'epoch', 'daily', 'sv',
+    :param list[str] svs: List of satellites to process
+    :param list clk_norm_types: Normalizations to apply to clocks. Available options include 'epoch', 'daily', 'sv',
             any satellite PRN, or any combination of them, defaults to empty list
     :return _pd.DataFrame: The Pandas DataFrame containing clock differences
     """
     base_clk_df = _gn_io.clk.read_clk(base_clk_file)
-    test_clk_df = _gn_io.clk.read_clk(test_clk_file)
+    mask = base_clk_df.index.get_level_values('CODE').isin(svs)
+    base_clk_df = base_clk_df[mask]
 
-    diff_clk_df = compare_clk(test_clk_df, base_clk_df, norm_types=norm_types)
+    test_clk_df = _gn_io.clk.read_clk(test_clk_file)
+    mask = test_clk_df.index.get_level_values('CODE').isin(svs)
+    test_clk_df = test_clk_df[mask]
 
+    diff_clk_df = compare_clk(test_clk_df, base_clk_df, norm_types=clk_norm_types)  # TODO Eugene: compare clocks by constellation
     diff_clk_df = diff_clk_df.stack(dropna=False).to_frame(name="Clock") * 1e9
     diff_clk_df["|Clock|"] = diff_clk_df.abs()
 

gnssanalysis/gn_io/clk.py

@@ -70,11 +70,17 @@ def read_clk(clk_path):
     return clk_df
 
 
-def get_AS_entries(clk_df: _pd.DataFrame) -> _pd.Series:
-    """fastest method to grab a specific category!, same as clk_df.EST.loc['AS'] but >6 times faster"""
-    AS_cat_code = clk_df.index.levels[0].categories.get_loc("AS")
-    mask = clk_df.index.codes[0] == AS_cat_code
-    return _pd.Series(data=clk_df.values[:, 0][mask], index=clk_df.index.droplevel(0)[mask])
+def get_sv_clocks(clk_df: _pd.DataFrame) -> _pd.Series:
+    """Retrieve satellite clocks from a CLK or SP3 dataframe"""
+    if clk_df.index.names == ['A', 'J2000', 'CODE']:
+        # fastest method to grab a specific category!, same as clk_df.EST.loc['AS'] but >6 times faster
+        AS_cat_code = clk_df.index.levels[0].categories.get_loc("AS")
+        mask = clk_df.index.codes[0] == AS_cat_code
+        return _pd.Series(data=clk_df.values[:, 0][mask], index=clk_df.index.droplevel(0)[mask])
+    elif clk_df.index.names == ['J2000', 'PRN']:
+        return _pd.Series(data=clk_df[("EST", "CLK")].values, index=clk_df.index)
+    else:
+        raise IndexError("Incorrect index names of dataframe")
 
 
 def rm_epoch_gnss_bias(clk_df_unst: _pd.DataFrame):

gnssanalysis/gn_io/sp3.py

@@ -809,17 +809,36 @@ def sp3merge(
     return merged_sp3
 
 
-def sp3_hlm_trans(a: _pd.DataFrame, b: _pd.DataFrame) -> tuple[_pd.DataFrame, list]:
+def hlm_trans(a: _np.ndarray, b: _np.ndarray) -> tuple[_np.ndarray, list]:
     """
-     Rotates sp3_b into sp3_a.
+    Rotates b into a.
 
-     :param DataFrame a: The sp3_a DataFrame.
-     :param DataFrame b : The sp3_b DataFrame.
+    :param _np.ndarray a: The a array.
+    :param _np.ndarray b: The b array.
 
-    :returntuple[pandas.DataFrame, list]: A tuple containing the updated sp3_b DataFrame and the HLM array with applied computed parameters and residuals.
+    :return tuple[_np.ndarray, list]: A tuple containing the output array and the HLM array with applied computed parameters and residuals.
     """
-    hlm = _gn_transform.get_helmert7(pt1=a.EST[["X", "Y", "Z"]].values, pt2=b.EST[["X", "Y", "Z"]].values)
-    b.iloc[:, :3] = _gn_transform.transform7(xyz_in=b.EST[["X", "Y", "Z"]].values, hlm_params=hlm[0])
+    hlm = _gn_transform.get_helmert7(pt1=a, pt2=b)
+    xyz_out = _gn_transform.transform7(xyz_in=b, hlm_params=hlm[0])
+    return xyz_out, hlm
+
+
+def sp3_hlm_trans(a: _pd.DataFrame, b: _pd.DataFrame, epochwise: bool = False) -> tuple[_pd.DataFrame, list]:
+    """
+    Rotates sp3_b into sp3_a.
+
+    :param DataFrame a: The sp3_a DataFrame.
+    :param DataFrame b : The sp3_b DataFrame.
+    :param bool epochwise: Epochwise HLM transformation.
+
+    :return tuple[pandas.DataFrame, list]: A tuple containing the updated sp3_b DataFrame and the HLM array with applied computed parameters and residuals.
+    """
+    if epochwise:
+        for epoch in b.index.get_level_values("J2000").unique():
+            b.loc[epoch].iloc[:, :3], hlm = hlm_trans(a.loc[epoch].EST[["X", "Y", "Z"]].values, b.loc[epoch].EST[["X", "Y", "Z"]].values)  # Eugene: hlm will be overwritten in the loop 
+    else:
+        b.iloc[:, :3], hlm = hlm_trans(a.EST[["X", "Y", "Z"]].values, b.EST[["X", "Y", "Z"]].values)
+
     return b, hlm
 
 
@@ -829,6 +848,7 @@ def diff_sp3_rac(
     sp3_test: _pd.DataFrame,
     hlm_mode: Literal[None, "ECF", "ECI"] = None,
     use_cubic_spline: bool = True,
+    epochwise: bool = False,
 ) -> _pd.DataFrame:
     """
     Computes the difference between the two sp3 files in the radial, along-track and cross-track coordinates.
@@ -837,6 +857,7 @@ def diff_sp3_rac(
     :param DataFrame sp3_test: The test sp3 DataFrame.
     :param string hlm_mode: The mode for HLM transformation. Can be None, "ECF", or "ECI".
     :param bool use_cubic_spline: Flag indicating whether to use cubic spline for velocity computation.
+    :param bool epochwise: Epochwise orbit comparison.
     :return: The DataFrame containing the difference in RAC coordinates.
     """
     hlm_modes = [None, "ECF", "ECI"]
@@ -852,11 +873,11 @@ def diff_sp3_rac(
 
     hlm = None  # init hlm var
     if hlm_mode == "ECF":
-        sp3_test, hlm = sp3_hlm_trans(sp3_baseline, sp3_test)
+        sp3_test, hlm = sp3_hlm_trans(sp3_baseline, sp3_test, epochwise)
     sp3_baseline_eci = _gn_transform.ecef2eci(sp3_baseline)
     sp3_test_eci = _gn_transform.ecef2eci(sp3_test)
     if hlm_mode == "ECI":
-        sp3_test_eci, hlm = sp3_hlm_trans(sp3_baseline_eci, sp3_test_eci)
+        sp3_test_eci, hlm = sp3_hlm_trans(sp3_baseline_eci, sp3_test_eci, epochwise)
 
     diff_eci = sp3_test_eci - sp3_baseline_eci
 

gnssanalysis/gn_transform.py

@@ -8,8 +8,13 @@
 from . import gn_const as _gn_const
 
 
-def gen_helm_aux(pt1, pt2):
+def gen_helm_aux(pt1, pt2, dropna = True):
     """aux function for helmert values inversion."""
+    if dropna:
+        mask = ~_np.isnan(pt1).any(axis=1) & ~_np.isnan(pt2).any(axis=1)
+        pt1 = pt1[mask]
+        pt2 = pt2[mask]
+
     pt1 = pt1.astype(float)
     pt2 = pt2.astype(float)
     n_points = pt1.shape[0]
@@ -32,23 +37,23 @@ def gen_helm_aux(pt1, pt2):
     return A, rhs
 
 
-def get_helmert7(pt1: _np.ndarray, pt2: _np.ndarray, scale_in_ppm: bool = True) -> Tuple[_np.ndarray, _np.ndarray]:
+def get_helmert7(pt1: _np.ndarray, pt2: _np.ndarray, scale_in_ppm: bool = True, dropna: bool = True) -> Tuple[_np.ndarray, _np.ndarray]:
     """Inversion of 7 Helmert parameters between 2 sets of points.
 
     :param numpy.ndarray pt1: The first set of points.
     :param numpy.ndarray pt2: The second set of points.
     :param bool scale_in_ppm: Whether the scale parameter is in parts per million (ppm).
+    :param bool dropna: Whether to drop NaN values in input data.
 
     :returns uple[np.ndarray, np.ndarray]: A tuple containing the Helmert parameters and the residuals.
-
     """
-    A, rhs = gen_helm_aux(pt1, pt2)
+    A, rhs = gen_helm_aux(pt1, pt2, dropna)
     sol = list(_np.linalg.lstsq(A, rhs, rcond=-1))  # parameters
+    res = rhs - A @ sol[0]  # computing residuals for pt2
+    sol.append(res.reshape(-1, 3))  # appending residuals
     sol[0] = sol[0].flatten() * -1.0  # flattening the HLM params arr to [Tx, Ty, Tz, Rx, Ry, Rz, Scale/mu]
     if scale_in_ppm:
         sol[0][-1] *= 1e6  # scale in ppm
-    res = rhs - A @ sol[0]  # computing residuals for pt2
-    sol.append(res.reshape(-1, 3))  # appending residuals
     return sol