Fix clock normalisation and convert sp3 position nodata values to NaNs

README.md

@@ -156,7 +156,7 @@ E03     -0.1216 0.042   0.1286
 ```
 These can also be stacked:
 ```bash
-clkq -i IGS2R03FIN_20191990000_01D_30S_CLK.CLK COD0R03FIN_20191990000_01D_30S_CLK.CLK --norm "daily"
+clkq -i IGS2R03FIN_20191990000_01D_30S_CLK.CLK COD0R03FIN_20191990000_01D_30S_CLK.CLK --norm "daily" --norm "epoch"
 ```
 Result:
 ```

gnssanalysis/filenames.py

@@ -588,7 +588,7 @@ def determine_sp3_name_props(file_path: pathlib.Path) -> Dict[str, Any]:
     """
     name_props = {}
     try:
-        sp3_df = gn_io.sp3.read_sp3(file_path)
+        sp3_df = gn_io.sp3.read_sp3(file_path, nodata_to_nan=False)
         props_from_existing_name = determine_name_props_from_filename(file_path.name)
         logging.debug(f"props_from_existing_name =\n{props_from_existing_name}")
         # name_props["analysis_center"] = sp3_df.attrs["HEADER"].HEAD.AC[0:3].upper().ljust(3,"X")

gnssanalysis/gn_diffaux.py

@@ -364,6 +364,9 @@ def compare_clk(
                 clk_a_unst=clk_a_unst, clk_b_unst=clk_b_unst
             )
 
+        clk_a_unst[clk_b_unst.isna()] = _np.nan  # replace corresponding values in clk_a_unst with NaN where clk_b_unst is NaN
+        clk_b_unst[clk_a_unst.isna()] = _np.nan  # replace corresponding values in clk_b_unst with NaN where clk_a_unst is NaN
+
             # get the sv to use for norm and overwrite norm_type value with sv prn code
         _logging.info("---removing common mode from clk 1---")
         _gn_io.clk.rm_clk_bias(clk_a_unst, norm_types=norm_types)
@@ -521,10 +524,10 @@ def sisre(
 
 
 def diffsp3(
-    sp3_a_path, sp3_b_path, tol, log_lvl, clk_a_path, clk_b_path, hlm_mode=None, plot=False, write_rac_file=False
+    sp3_a_path, sp3_b_path, tol, log_lvl, clk_a_path, clk_b_path, nodata_to_nan=True, hlm_mode=None, plot=False, write_rac_file=False
 ):
     """Compares two sp3 files and outputs a dataframe of differences above tolerance if such were found"""
-    sp3_a, sp3_b = _gn_io.sp3.read_sp3(sp3_a_path), _gn_io.sp3.read_sp3(sp3_b_path)
+    sp3_a, sp3_b = _gn_io.sp3.read_sp3(sp3_a_path, nodata_to_nan=nodata_to_nan), _gn_io.sp3.read_sp3(sp3_b_path, nodata_to_nan=nodata_to_nan)
 
     as_sisre = False  # the switch only needed for logging msg
     clk_a = clk_b = None

gnssanalysis/gn_io/sp3.py

@@ -39,19 +39,41 @@
 SP3_CLOCK_NODATA_STRING = " 999999.999999"  # Not used for reading, as fractional components are optional
 SP3_CLOCK_NODATA_NUMERIC = 999999
 SP3_POS_NODATA_STRING = "      0.000000"
+SP3_POS_NODATA_NUMERIC = 0
 SP3_CLOCK_STD_NODATA = -1000
 SP3_POS_STD_NODATA = -100
 
 
-def sp3_clock_nodata_to_nan(sp3_df):
+def sp3_pos_nodata_to_nan(
+    sp3_df: _pd.DataFrame
+) -> None:
     """
-    Converts the SP3 Clock column's nodata value (999999 or 999999.999999 - the fractional component optional) to NaNs.
+    Converts the SP3 Positional column's nodata values (0.000000) to NaNs.
     See https://files.igs.org/pub/data/format/sp3_docu.txt
+
+    :param _pd.DataFrame sp3_df: SP3 data frame to filter nodata values for
+    :return None
     """
-    nan_mask = sp3_df.iloc[:, 1:5].values >= SP3_CLOCK_NODATA_NUMERIC
-    nans = nan_mask.astype(float)
-    nans[nan_mask] = _np.NAN
-    sp3_df.iloc[:, 1:5] = sp3_df.iloc[:, 1:5].values + nans
+    nan_mask = (
+        (sp3_df[("EST", "X")] == SP3_POS_NODATA_NUMERIC)
+        & (sp3_df[("EST", "Y")] == SP3_POS_NODATA_NUMERIC)
+        & (sp3_df[("EST", "Z")] == SP3_POS_NODATA_NUMERIC)
+    )
+    sp3_df.loc[nan_mask, [("EST", "X"), ("EST", "Y"), ("EST", "Z")]] = _np.NAN
+
+
+def sp3_clock_nodata_to_nan(
+    sp3_df: _pd.DataFrame
+) -> None:
+    """
+    Converts the SP3 Clock column's nodata values (999999 or 999999.999999 - the fractional component optional) to NaNs.
+    See https://files.igs.org/pub/data/format/sp3_docu.txt
+
+    :param _pd.DataFrame sp3_df: SP3 data frame to filter nodata values for
+    :return None
+    """
+    nan_mask = sp3_df[("EST", "CLK")] >= SP3_CLOCK_NODATA_NUMERIC
+    sp3_df.loc[nan_mask, ("EST", "CLK")] = _np.NAN
 
 
 def mapparm(old, new):
@@ -63,7 +85,7 @@ def mapparm(old, new):
     return off, scl
 
 
-def read_sp3(sp3_path, pOnly=True):
+def read_sp3(sp3_path, pOnly=True, nodata_to_nan=True):
     """Read SP3 file
     Returns STD values converted to proper units (mm/ps) also if present.
     by default leaves only P* values (positions), removing the P key itself
@@ -144,7 +166,10 @@ def read_sp3(sp3_path, pOnly=True):
         ]
         sp3_df.STD = base_xyzc**sp3_df.STD.values
 
-    sp3_clock_nodata_to_nan(sp3_df)  # Convert 999999* (which indicates nodata in the SP3 CLK column) to NaN
+    if nodata_to_nan:
+        sp3_pos_nodata_to_nan(sp3_df)  # Convert 0.000000 (which indicates nodata in the SP3 POS column) to NaN
+        sp3_clock_nodata_to_nan(sp3_df)  # Convert 999999* (which indicates nodata in the SP3 CLK column) to NaN
+
     if pOnly or parsed_header.HEAD.loc["PV_FLAG"] == "P":
         pMask = series.astype("S1") == b"P"
         sp3_df = sp3_df[pMask].set_index([dt_index[pMask], series.str[1:4].values[pMask].astype(str).astype(object)])
@@ -458,9 +483,9 @@ def merge_attrs(df_list):
     return _pd.Series(_np.concatenate([head, sv_info]), index=df.index)
 
 
-def sp3merge(sp3paths, clkpaths=None):
+def sp3merge(sp3paths, clkpaths=None, nodata_to_nan=False):
     """Reads in a list of sp3 files and optianl list of clk file and merges them into a single sp3 file"""
-    sp3_dfs = [read_sp3(sp3_file) for sp3_file in sp3paths]
+    sp3_dfs = [read_sp3(sp3_file, nodata_to_nan=nodata_to_nan) for sp3_file in sp3paths]
     merged_sp3 = _pd.concat(sp3_dfs)
     merged_sp3.attrs["HEADER"] = merge_attrs(sp3_dfs)
 

gnssanalysis/gn_utils.py

@@ -213,6 +213,13 @@ def clk(ctx, norm):
     help="path to aux2 file",
     show_default=True,
 )
+@_click.option(
+    "--nodata-to-nan",
+    type=bool,
+    help="convert nodata values (0.000000 for POS, 999999 or 999999.999999 for CLK) to NaNs. Default: True",
+    default=True,
+    show_default=True,
+)
 @_click.option(
     "--hlm_mode",
     type=_click.Choice(["ECF", "ECI"], case_sensitive=False),
@@ -225,7 +232,7 @@ def clk(ctx, norm):
     is_flag=True,
     help="outputs Radial/Along-track/Cross-track into a file",
 )
-def sp3(ctx, aux1, aux2, hlm_mode, rac):  # no coef
+def sp3(ctx, aux1, aux2, nodata_to_nan, hlm_mode, rac):  # no coef
     from .gn_diffaux import diffsp3
 
     diffutil_verify_input(ctx.parent.params["input"])
@@ -236,6 +243,7 @@ def sp3(ctx, aux1, aux2, hlm_mode, rac):  # no coef
         clk_b_path=aux2,
         tol=ctx.parent.params["atol"],
         log_lvl=ctx.parent.params["log_lvl"],
+        nodata_to_nan=nodata_to_nan,
         hlm_mode=hlm_mode,
         plot=ctx.parent.params["plot"],
         write_rac_file=rac,
@@ -305,14 +313,20 @@ def snxmap(sinexpaths, outdir):
     help="output path",
     default=_os.curdir + "/merge.sp3",
 )
-def sp3merge(sp3paths, clkpaths, output):
+@_click.option(
+    "--nodata-to-nan",
+    type=bool,
+    help="convert nodata values (0.000000 for POS, 999999 or 999999.999999 for CLK) to NaNs. Default: False",
+    default=False,
+)
+def sp3merge(sp3paths, clkpaths, output, nodata_to_nan):
     """
     sp3 files paths to merge, Optional clock files which is useful to insert clk offset values into sp3 file.
     """
     from .gn_io import sp3
 
     _logging.info(msg=output)
-    merged_df = sp3.sp3merge(sp3paths=sp3paths, clkpaths=clkpaths)
+    merged_df = sp3.sp3merge(sp3paths=sp3paths, clkpaths=clkpaths, nodata_to_nan=nodata_to_nan)
     sp3.write_sp3(sp3_df=merged_df, path=output)
 
 
@@ -488,6 +502,13 @@ def trace2mongo(trace_paths, db_name):
     default="table",
     help="If JSON format chosen, choose how the output JSON schema is formated. Default is 'table'. Options: 'table', 'split', 'records', 'index', 'columns', 'values'",
 )
+@_click.option(
+    "--nodata-to-nan",
+    type=bool,
+    help="convert nodata values (0.000000 for POS, 999999 or 999999.999999 for CLK) to NaNs. Default: True",
+    default=True,
+    show_default=True,
+)
 @_click.option(
     "-h",
     "--hlm_mode",
@@ -535,7 +556,7 @@ def trace2mongo(trace_paths, db_name):
     default=None,
     show_default=True,
 )
-def orbq(input, output_path, format, csv_separation, json_format, hlm_mode, satellite, constellation, header, index, reject_re):
+def orbq(input, output_path, format, csv_separation, json_format, nodata_to_nan, hlm_mode, satellite, constellation, header, index, reject_re):
     """
     A simple utility to assess pairs of sp3 files
     """
@@ -547,8 +568,8 @@ def orbq(input, output_path, format, csv_separation, json_format, hlm_mode, sate
     # else:
     #     _logging.disable()
 
-    sp3_a = gn_io.sp3.read_sp3(input[0])
-    sp3_b = gn_io.sp3.read_sp3(input[1])
+    sp3_a = gn_io.sp3.read_sp3(input[0], nodata_to_nan=nodata_to_nan)
+    sp3_b = gn_io.sp3.read_sp3(input[1], nodata_to_nan=nodata_to_nan)
     if reject_re is not None:
         logger.log(msg=f"Excluding satellites based on regex expression: '{reject_re}'", level=_logging.INFO)
         reject_mask = sp3_a.index.get_level_values(1).str.match(reject_re)