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,47 @@ def format_index(
 def sp3_difference(
     base_sp3_file: _Path,
     test_sp3_file: _Path,
+    svs: list[str],
+    orb_ref_frame: Literal["ECF", "ECI"] = "ECF",
+    orb_hlm_mode: Literal[None, "ECF", "ECI"] = None,
+    epochwise_hlm: bool = False,
+    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 bool epochwise_hlm: Epochwise Helmert transformation
+    :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, ref_frame=orb_ref_frame, hlm_mode=orb_hlm_mode, epochwise_hlm=epochwise_hlm)  # 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 +759,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,22 @@ 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
+
+    :param _pd.DataFrame clk_df: CLK or SP3 dataframe where to retreive satellite clocks from
+    :raises IndexError: Raise error if the dataframe is not indexed correctly
+    :return _pd.Series: Retrieved satellite clocks
+    """
+    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

@@ -214,7 +214,6 @@ def _process_sp3_block(
 ) -> _pd.DataFrame:
     """Process a single block of SP3 data.
 
-
     :param    str date: The date of the SP3 data block.
     :param    str data: The SP3 data block.
     :param    List[int] widths: The widths of the columns in the SP3 data block.
@@ -251,7 +250,6 @@ def read_sp3(
 ) -> _pd.DataFrame:
     """Reads an SP3 file and returns the data as a pandas DataFrame.
 
-
     :param str sp3_path: The path to the SP3 file.
     :param bool pOnly: If True, only P* values (positions) are included in the DataFrame. Defaults to True.
     :param bool nodata_to_nan: If True, converts 0.000000 (indicating nodata) to NaN in the SP3 POS column
@@ -448,6 +446,20 @@ def parse_sp3_header(header: bytes, warn_on_negative_sv_acc_values: bool = True)
     return _pd.concat([sp3_heading, sv_tbl], keys=["HEAD", "SV_INFO"], axis=0)
 
 
+def clean_sp3_orb(sp3_df: _pd.DataFrame) -> _pd.DataFrame:
+    sp3_df = sp3_df.filter(items=[("EST", "X"), ("EST", "Y"), ("EST", "Z")])
+
+    # Drop any duplicates in the index
+    sp3_df = sp3_df[~sp3_df.index.duplicated(keep="first")]
+
+    # Drop satellites that have any NaN values
+    nan_mask = sp3_df.isna()
+    nan_prns = nan_mask.groupby("PRN").any().any(axis=1)
+    sp3_df_cleaned = sp3_df[~sp3_df.index.get_level_values("PRN").isin(nan_prns[nan_prns].index)]
+
+    return sp3_df_cleaned
+
+
 def getVelSpline(sp3Df: _pd.DataFrame) -> _pd.DataFrame:
     """Returns the velocity spline of the input dataframe.
 
@@ -474,7 +486,6 @@ def getVelPoly(sp3Df: _pd.DataFrame, deg: int = 35) -> _pd.DataFrame:
     :param DataFrame sp3Df: A pandas DataFrame containing the sp3 data.
     :param int deg: Degree of the polynomial fit. Default is 35.
     :return DataFrame: A pandas DataFrame with the interpolated velocities added as a new column.
-
     """
     est = sp3Df.unstack(1).EST[["X", "Y", "Z"]]
     x = est.index.get_level_values("J2000").values
@@ -580,7 +591,6 @@ def gen_sp3_content(
     Organises, formats (including nodata values), then writes out SP3 content to a buffer if provided, or returns
     it otherwise.
 
-    Args:
     :param pandas.DataFrame sp3_df: The DataFrame containing the SP3 data.
     :param bool sort_outputs: Whether to sort the outputs. Defaults to False.
     :param io.TextIOBase  buf: The buffer to write the SP3 content to. Defaults to None.
@@ -792,7 +802,6 @@ def sp3merge(
     :param List[str] sp3paths: The list of paths to the sp3 files.
     :param Union[List[str], None] clkpaths: The list of paths to the clk files, or None if no clk files are provided.
     :param bool nodata_to_nan: Flag indicating whether to convert nodata values to NaN.
-
     :return pd.DataFrame: The merged sp3 DataFrame.
     """
     sp3_dfs = [read_sp3(sp3_file, nodata_to_nan=nodata_to_nan) for sp3_file in sp3paths]
@@ -809,25 +818,46 @@ def sp3merge(
     return merged_sp3
 
 
-def sp3_hlm_trans(a: _pd.DataFrame, b: _pd.DataFrame) -> tuple[_pd.DataFrame, list]:
+def hlm_trans(pt1: _np.ndarray, pt2: _np.ndarray) -> tuple[_np.ndarray, list]:
     """
-     Rotates sp3_b into sp3_a.
+    Rotates a set of points pt1 into pt2.
+
+    :param _np.ndarray pt1: The first set of points.
+    :param _np.ndarray pt2: The second set of points.
+    :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, pt2)
+    xyz_out = _gn_transform.transform7(xyz_in=pt2, hlm_params=hlm[0])
+    return xyz_out, hlm
 
-     :param DataFrame a: The sp3_a DataFrame.
-     :param DataFrame b : The sp3_b DataFrame.
 
-    :returntuple[pandas.DataFrame, list]: A tuple containing the updated sp3_b DataFrame and the HLM array with applied computed parameters and residuals.
+def sp3_hlm_trans(a: _pd.DataFrame, b: _pd.DataFrame, epochwise: bool = False) -> tuple[_pd.DataFrame, list]:
     """
-    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])
+    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.
+    """
+    hlm = []
+    if epochwise:
+        for epoch in b.index.get_level_values("J2000").unique():
+            b.loc[epoch, [("EST", "X"), ("EST", "Y"), ("EST", "Z")]], hlm_epoch = hlm_trans(a.loc[epoch].EST[["X", "Y", "Z"]].values, b.loc[epoch].EST[["X", "Y", "Z"]].values)
+            hlm.append(hlm_epoch)
+    else:
+        b.iloc[:, :3], hlm = hlm_trans(a.EST[["X", "Y", "Z"]].values, b.EST[["X", "Y", "Z"]].values)
+
     return b, hlm
 
 
 # TODO: move to gn_diffaux.py (and other associated functions as well)?
 def diff_sp3_rac(
     sp3_baseline: _pd.DataFrame,
     sp3_test: _pd.DataFrame,
+    ref_frame: Literal["ECF", "ECI"] = "ECF",
     hlm_mode: Literal[None, "ECF", "ECI"] = None,
+    epochwise_hlm: bool = False,
     use_cubic_spline: bool = True,
 ) -> _pd.DataFrame:
     """
@@ -836,27 +866,40 @@ def diff_sp3_rac(
     :param DataFrame sp3_baseline: The baseline sp3 DataFrame.
     :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 epochwise_hlm: Epochwise HLM transformation.
     :param bool use_cubic_spline: Flag indicating whether to use cubic spline for velocity computation.
     :return: The DataFrame containing the difference in RAC coordinates.
     """
+    ref_frames = ["ECF", "ECI"]
+    if ref_frame not in ref_frames:
+        raise ValueError(f"Invalid ref_frame. Expected one of: {ref_frames}")
+
     hlm_modes = [None, "ECF", "ECI"]
     if hlm_mode not in hlm_modes:
         raise ValueError(f"Invalid hlm_mode. Expected one of: {hlm_modes}")
 
-    # Drop any duplicates in the index
-    sp3_baseline = sp3_baseline[~sp3_baseline.index.duplicated(keep="first")]
-    sp3_test = sp3_test[~sp3_test.index.duplicated(keep="first")]
+    sp3_baseline = clean_sp3_orb(sp3_baseline)
+    sp3_test = clean_sp3_orb(sp3_test)
+
     # Ensure the test file is time-ordered so when we align the resulting dataframes will be time-ordered
     sp3_baseline = sp3_baseline.sort_index(axis="index", level="J2000")
     sp3_baseline, sp3_test = sp3_baseline.align(sp3_test, join="inner", axis=0)
 
+    # TODO Eugene: improve following code
     hlm = None  # init hlm var
     if hlm_mode == "ECF":
-        sp3_test, hlm = sp3_hlm_trans(sp3_baseline, sp3_test)
-    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_baseline_ecf = _gn_transform.eci2ecef(sp3_baseline) if ref_frame == "ECI" else sp3_baseline
+        sp3_test_ecf = _gn_transform.eci2ecef(sp3_test) if ref_frame == "ECI" else sp3_test
+        sp3_test_ecf, hlm = sp3_hlm_trans(sp3_baseline_ecf, sp3_test_ecf, epochwise_hlm)
+        sp3_baseline_eci = _gn_transform.ecef2eci(sp3_baseline_ecf)
+        sp3_test_eci = _gn_transform.ecef2eci(sp3_test_ecf)
+    elif hlm_mode == "ECI":
+        sp3_baseline_eci = _gn_transform.ecef2eci(sp3_baseline) if ref_frame == "ECF" else sp3_baseline
+        sp3_test_eci = _gn_transform.ecef2eci(sp3_test) if ref_frame == "ECF" else sp3_test
+        sp3_test_eci, hlm = sp3_hlm_trans(sp3_baseline_eci, sp3_test_eci, epochwise_hlm)
+    else:
+        sp3_baseline_eci = _gn_transform.ecef2eci(sp3_baseline) if ref_frame == "ECF" else sp3_baseline
+        sp3_test_eci = _gn_transform.ecef2eci(sp3_test) if ref_frame == "ECF" else sp3_test
 
     diff_eci = sp3_test_eci - sp3_baseline_eci