NPI-3655 Bring in remove_offline_sats(), add option in SP3 diff processing to avoid crashing cubic interpolator for velocity

gnssanalysis/gn_io/sp3.py

@@ -191,6 +191,41 @@ def sp3_clock_nodata_to_nan(sp3_df: _pd.DataFrame) -> None:
     sp3_df.loc[nan_mask, ("EST", "CLK")] = _np.nan
 
 
+def remove_offline_sats(sp3_df: _pd.DataFrame, df_friendly_name: str = ""):
+    """
+    Remove any satellites that have "0.0" or NaN for all three position coordinates - this indicates satellite offline.
+    Note that this removes a satellite which has *any* missing coordinate data, meaning a satellite with partial
+    observations will get removed entirely.
+    Added in version 0.0.57
+
+    :param _pd.DataFrame sp3_df: SP3 DataFrame to remove offline / nodata satellites from
+    :param str df_friendly_name: Name to use when referring to the DataFrame in log output (purely for clarity). Empty
+           string by default.
+    :return _pd.DataFrame: the SP3 DataFrame, with the offline / nodata satellites removed
+
+    """
+    # Get all entries with missing positions (i.e. nodata value, represented as either 0 or NaN), then get the
+    # satellite names (SVs) of these and dedupe them, giving a list of all SVs with one or more missing positions.
+
+    # Mask for nodata POS values in raw form:
+    mask_zero = (sp3_df.EST.X == 0.0) & (sp3_df.EST.Y == 0.0) & (sp3_df.EST.Z == 0.0)
+    # Mask for converted values, as NaNs:
+    mask_nan = (sp3_df.EST.X.isna()) & (sp3_df.EST.Y.isna()) & (sp3_df.EST.Z.isna())
+    mask_either = _np.logical_or(mask_zero, mask_nan)
+
+    # With mask, filter for entries with no POS value, then get the sat name (SVs) from the entry, then dedupe:
+    offline_sats = sp3_df[mask_either].index.get_level_values(1).unique()
+
+    # Using that list of offline / partially offline sats, remove all entries for those sats from the SP3 DataFrame:
+    sp3_df = sp3_df.drop(offline_sats, level=1, errors="ignore")
+    sp3_df.attrs["HEADER"].HEAD.ORB_TYPE = "INT"  # Allow the file to be read again by read_sp3 - File ORB_TYPE changes
+    if len(offline_sats) > 0:
+        logger.info(f"Dropped offline / nodata sats from {df_friendly_name} SP3 DataFrame: {offline_sats.values}")
+    else:
+        logger.info(f"No offline / nodata sats detected to be dropped from {df_friendly_name} SP3 DataFrame")
+    return sp3_df
+
+
 def mapparm(old: Tuple[float, float], new: Tuple[float, float]) -> Tuple[float, float]:
     """
     Evaluate the offset and scale factor needed to map values from the old range to the new range.
@@ -247,7 +282,10 @@ def description_for_path_or_bytes(path_or_bytes: Union[str, Path, bytes]) -> Opt
 
 
 def read_sp3(
-    sp3_path_or_bytes: Union[str, Path, bytes], pOnly: bool = True, nodata_to_nan: bool = True
+    sp3_path_or_bytes: Union[str, Path, bytes],
+    pOnly: bool = True,
+    nodata_to_nan: bool = True,
+    drop_offline_sats: bool = False,
 ) -> _pd.DataFrame:
     """Reads an SP3 file and returns the data as a pandas DataFrame.
 
@@ -256,6 +294,8 @@ def read_sp3(
     :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
             and converts 999999* (indicating nodata) to NaN in the SP3 CLK column. Defaults to True.
+    :param bool drop_offline_sats: If True, drops satellites from the DataFrame if they have ANY missing (nodata)
+            values in the SP3 POS column.
     :return pandas.DataFrame: The SP3 data as a DataFrame.
     :raise FileNotFoundError: If the SP3 file specified by sp3_path_or_bytes does not exist.
     :raise Exception: For other errors reading SP3 file/bytes
@@ -292,6 +332,8 @@ def read_sp3(
     date_lines, data_blocks = _split_sp3_content(content)
     sp3_df = _pd.concat([_process_sp3_block(date, data) for date, data in zip(date_lines, data_blocks)])
     sp3_df = _reformat_df(sp3_df)
+    if drop_offline_sats:
+        sp3_df = remove_offline_sats(sp3_df)
     if nodata_to_nan:
         # Convert 0.000000 (which indicates nodata in the SP3 POS column) to NaN
         sp3_pos_nodata_to_nan(sp3_df)
@@ -454,6 +496,8 @@ def getVelSpline(sp3Df: _pd.DataFrame) -> _pd.DataFrame:
     :param DataFrame sp3Df: The input dataframe containing position data.
     :return DataFrame: The dataframe containing the velocity spline.
 
+    :caution :This function cannot handle *any* NaN / nodata / non-finite position values. By contrast, getVelPoly()
+              is more forgiving, but accuracy of results, particulary in the presence of NaNs, has not been assessed.
     :note :The velocity is returned in the same units as the input dataframe, e.g. km/s (needs to be x10000 to be in cm as per sp3 standard).
     """
     sp3dfECI = sp3Df.EST.unstack(1)[["X", "Y", "Z"]]  # _ecef2eci(sp3df)
@@ -849,14 +893,20 @@ def diff_sp3_rac(
     sp3_test: _pd.DataFrame,
     hlm_mode: Literal[None, "ECF", "ECI"] = None,
     use_cubic_spline: bool = True,
+    use_offline_sat_removal: bool = False,
 ) -> _pd.DataFrame:
     """
     Computes the difference between the two sp3 files in the radial, along-track and cross-track coordinates.
 
     :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 use_cubic_spline: Flag indicating whether to use cubic spline for velocity computation.
+    :param bool use_cubic_spline: Flag indicating whether to use cubic spline for velocity computation. Caution: cubic
+           spline interpolation does not tolerate NaN / nodata values. Consider enabling use_offline_sat_removal if
+           using cubic spline, or alternatively use poly interpolation by setting use_cubic_spline to False.
+    :param bool use_offline_sat_removal: Flag indicating whether to remove satellites which are offline / have some
+           nodata position values. Caution: ensure you turn this on if using cubic spline interpolation with data
+           which may have holes in it (nodata).
     :return: The DataFrame containing the difference in RAC coordinates.
     """
     hlm_modes = [None, "ECF", "ECI"]
@@ -866,6 +916,20 @@ def diff_sp3_rac(
     # 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")]
+
+    if use_cubic_spline and not use_offline_sat_removal:
+        logger.warning(
+            "Caution: use_cubic_spline is enabled, but use_offline_sat_removal is not. If there are any nodata "
+            "position values, the cubic interpolator will crash!"
+        )
+    # Drop any satellites (SVs) which are offline or partially offline.
+    # Note: this currently removes SVs with ANY nodata values for position, so a single glitch will remove
+    # the SV from the whole file.
+    # This step was added after velocity interpolation failures due to non-finite (NaN) values from offline SVs.
+    if use_offline_sat_removal:
+        sp3_baseline = remove_offline_sats(sp3_baseline, df_friendly_name="baseline")
+        sp3_test = remove_offline_sats(sp3_test, df_friendly_name="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)

tests/test_datasets/sp3_test_data.py

@@ -137,3 +137,34 @@
 PR02   5096.267067   9515.396007 -23066.803522    -23.760149
 EOF
 """
+
+# Original filename IGS0DEMULT_20243181800_02D_05M_ORB.SP3
+# Modified using utility script, header manually updated to match.
+sp3_test_data_partially_offline_sat = b"""#dP2024 11 13 18  0  0.0000000        3 ORBIT IGS20 HLM  IGS
+## 2340 324000.00000000   300.00000000 60627 0.7500000000000
++    3   G02G03G19  0  0  0  0  0  0  0  0  0  0  0  0  0  0
+++         4  4  3  0  0  0  0  0  0  0  0  0  0  0  0  0  0
+%c M  cc GPS ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc
+%c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc
+%f  1.2500000  1.025000000  0.00000000000  0.000000000000000
+%f  0.0000000  0.000000000  0.00000000000  0.000000000000000
+%i    0    0    0    0      0      0      0      0         0
+%i    0    0    0    0      0      0      0      0         0
+/*
+/*  NOTE: This is a highly modified file for testing only.
+/*  It was originally derived from an IGS product, but
+/*  should not be taken as authoritative.
+*  2024 11 13 18  0  0.00000000
+PG02 -22358.852114 -13824.987069  -4481.658280   -315.353433  8  9 10           
+PG03 -12211.683721 -14412.444793  18550.880326    598.854905  9  9  9           
+PG19 -16358.034086  16602.996553  12179.863169    562.338315  8 10  8           
+*  2024 11 13 18  5  0.00000000
+PG02 -22143.724430 -13804.799870  -5406.496217   -315.350875  3  9  9           
+PG03 -12184.752097 -15116.241903  18008.487790    598.858009  8  8 10           
+PG19      0.000000      0.000000      0.000000 999999.999999                    
+*  2024 11 13 18 10  0.00000000
+PG02 -21896.123567 -13775.172104  -6321.083787   -315.348317  5  9  9           
+PG03 -12170.453519 -15798.978021  17431.310591    598.861113  8  8 10           
+PG19      0.000000      0.000000      0.000000 999999.999999                    
+EOF
+"""

tests/test_sp3.py

@@ -12,6 +12,7 @@
     sp3_test_data_igs_benchmark_null_clock as input_data,
     # second dataset is a truncated version of file COD0OPSFIN_20242010000_01D_05M_ORB.SP3:
     sp3_test_data_truncated_cod_final as input_data2,
+    sp3_test_data_partially_offline_sat as offline_sat_test_data,
 )
 
 
@@ -180,6 +181,22 @@ def test_velinterpolation(self, mock_file):
         self.assertIsNotNone(r)
         self.assertIsNotNone(r2)
 
+    @patch("builtins.open", new_callable=mock_open, read_data=offline_sat_test_data)
+    def test_sp3_offline_sat_removal(self, mock_file):
+        sp3_df = sp3.read_sp3("mock_path", pOnly=False)
+        self.assertEqual(
+            sp3_df.index.get_level_values(1).unique().array,
+            ["G02", "G03", "G19"],
+            "Should be three SVs in test file before removing offline ones",
+        )
+
+        sp3_df = sp3.remove_offline_sats(sp3_df)
+        self.assertEqual(
+            sp3_df.index.get_level_values(1).unique().array,
+            ["G02", "G03"],
+            "Should be two SVs after removing offline ones",
+        )
+
 
 class TestMergeSP3(TestCase):
     def setUp(self):