From f46e274236d626d4b119341fb3371e044f4252cd Mon Sep 17 00:00:00 2001 From: Sebastien Date: Fri, 19 Jul 2024 16:04:09 +1000 Subject: [PATCH 01/15] test: adding some of my easy unittest for read sp3. --- gnssanalysis/gn_io/sp3.py | 5 ++-- tests/test_sp3.py | 59 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 2 deletions(-) create mode 100644 tests/test_sp3.py diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index b5c43ef..0b03bd1 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -54,7 +54,7 @@ def sp3_pos_nodata_to_nan(sp3_df: _pd.DataFrame) -> None: & (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 + sp3_df.loc[nan_mask, [("EST", "X"), ("EST", "Y"), ("EST", "Z")]] = _np.nan def sp3_clock_nodata_to_nan(sp3_df: _pd.DataFrame) -> None: @@ -66,7 +66,7 @@ def sp3_clock_nodata_to_nan(sp3_df: _pd.DataFrame) -> None: :return None """ nan_mask = sp3_df[("EST", "CLK")] >= SP3_CLOCK_NODATA_NUMERIC - sp3_df.loc[nan_mask, ("EST", "CLK")] = _np.NAN + sp3_df.loc[nan_mask, ("EST", "CLK")] = _np.nan def mapparm(old, new): @@ -153,6 +153,7 @@ def read_sp3(sp3_path, pOnly=True, nodata_to_nan=True): 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 + # print(sp3_df.index.has_duplicates()) if pOnly or parsed_header.HEAD.loc["PV_FLAG"] == "P": sp3_df = sp3_df[sp3_df.index.get_level_values("PV_FLAG") == "P"] diff --git a/tests/test_sp3.py b/tests/test_sp3.py new file mode 100644 index 0000000..5e926aa --- /dev/null +++ b/tests/test_sp3.py @@ -0,0 +1,59 @@ +import unittest +from unittest.mock import patch, mock_open + +import numpy as np + +import gnssanalysis.gn_io.sp3 as sp3 + +# dataset is part of the IGS benchmark (modified to include non null data on clock) +input_data = b"""#dV2007 4 12 0 0 0.00000000 289 ORBIT IGS14 BHN ESOC +## 1422 345600.00000000 900.00000000 54202 0.0000000000000 ++ 2 G01G02 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +++ 8 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +++ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +++ 0 0 0 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 0.0000000 0.000000000 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 +/* EUROPEAN SPACE OPERATIONS CENTRE - DARMSTADT, GERMANY +/* --------------------------------------------------------- +/* SP3 FILE GENERATED BY NAPEOS BAHN TOOL (DETERMINATION) +/* PCV:IGS14_2022 OL/AL:EOT11A NONE YN ORB:CoN CLK:CoN +* 2007 4 12 0 0 0.00000000 +PG01 -6114.801556 -13827.040252 22049.171610 999999.999999 +VG01 27184.457428 -3548.055474 5304.058806 999999.999999 +PG02 12947.223282 22448.220655 6215.570741 999999.999999 +VG02 -7473.756152 -4355.288568 29939.333728 999999.999999 +* 2007 4 12 0 15 0.00000000 +PG01 -3659.032812 -14219.662913 22339.175481 123456.999999 +VG01 27295.435569 -5170.061971 1131.227754 999999.999999 +PG02 12163.580358 21962.803659 8849.429007 999999.999999 +VG02 -9967.334764 -6367.969150 28506.683280 999999.999999 +* 2007 4 12 0 30 0.00000000 +PG01 -1218.171155 -14755.013599 22252.168480 999999.999999 +VG01 26855.435366 -6704.236117 -3062.394499 999999.999999 +PG02 11149.555664 21314.099837 11331.977499 123456.999999 +VG02 -12578.915944 -7977.396362 26581.116225 999999.999999 +EOF""" + + +class TestSp3(unittest.TestCase): + @patch("builtins.open", new_callable=mock_open, read_data=input_data) + def test_read_sp3_pOnly(self, mock_file): + result = sp3.read_sp3("mock_path", pOnly=True) + self.assertEqual(len(result), 6) + + @patch("builtins.open", new_callable=mock_open, read_data=input_data) + def test_read_sp3_pv(self, mock_file): + result = sp3.read_sp3("mock_path", pOnly=False) + self.assertTrue(len(result) == 12) + self.assertEqual((np.isnan(result[("EST", "CLK")])).sum(), 10) From 19e68342f423cce341893f06615161906be80ba8 Mon Sep 17 00:00:00 2001 From: Sebastien Date: Fri, 19 Jul 2024 16:04:55 +1000 Subject: [PATCH 02/15] clean: removing print statement --- gnssanalysis/gn_io/sp3.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 0b03bd1..78b4ec2 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -80,7 +80,6 @@ def mapparm(old, new): def _process_sp3_block(date, data, widths, names): """Process a single block of SP3 data""" - print(data) if not data or len(data) == 0: return _pd.DataFrame() epochs_dt = _pd.to_datetime(_pd.Series(date).str.slice(2, 21).values.astype(str), format=r"%Y %m %d %H %M %S") From ce2b88585659502af8decde80d3ca16c7f74dc09 Mon Sep 17 00:00:00 2001 From: Sebastien Date: Mon, 22 Jul 2024 11:55:50 +1000 Subject: [PATCH 03/15] add: extra small unittest --- tests/test_sp3.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/tests/test_sp3.py b/tests/test_sp3.py index 5e926aa..0449990 100644 --- a/tests/test_sp3.py +++ b/tests/test_sp3.py @@ -2,6 +2,7 @@ from unittest.mock import patch, mock_open import numpy as np +import pandas as pd import gnssanalysis.gn_io.sp3 as sp3 @@ -57,3 +58,26 @@ def test_read_sp3_pv(self, mock_file): result = sp3.read_sp3("mock_path", pOnly=False) self.assertTrue(len(result) == 12) self.assertEqual((np.isnan(result[("EST", "CLK")])).sum(), 10) + + def test_sp3_clock_nodata_to_nan(self): + sp3_df = pd.DataFrame({("EST", "CLK"): [999999.999999, 123456.789, 999999.999999, 987654.321]}) + sp3.sp3_clock_nodata_to_nan(sp3_df) + expected_result = pd.DataFrame({("EST", "CLK"): [np.nan, 123456.789, np.nan, 987654.321]}) + print(expected_result) + print(sp3_df) + print(sp3_df.equals(expected_result)) + self.assertTrue(sp3_df.equals(expected_result)) + + def test_sp3_pos_nodata_to_nan(self): + sp3_df = pd.DataFrame( + {("EST", "X"): [0.0, 1.0, 0.0, 2.0], ("EST", "Y"): [0.0, 0.0, 0.0, 2.0], ("EST", "Z"): [0.0, 1.0, 0.0, 0.0]} + ) + sp3.sp3_pos_nodata_to_nan(sp3_df) + expected_result = pd.DataFrame( + { + ("EST", "X"): [np.nan, 1.0, np.nan, 2.0], + ("EST", "Y"): [np.nan, 0.0, np.nan, 2.0], + ("EST", "Z"): [np.nan, 1.0, np.nan, 0.0], + } + ) + self.assertTrue(sp3_df.equals(expected_result)) From eb74bb526371350c838e22086662d9b3d6780dee Mon Sep 17 00:00:00 2001 From: Sebastien Date: Mon, 22 Jul 2024 11:55:50 +1000 Subject: [PATCH 04/15] doc: adding doc --- gnssanalysis/gn_io/sp3.py | 88 ++++++++++++++++++++++++++++++++------- tests/test_sp3.py | 24 +++++++++++ 2 files changed, 97 insertions(+), 15 deletions(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 78b4ec2..4c1a40a 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -2,7 +2,7 @@ import io as _io import os as _os import re as _re -from typing import Literal, Union +from typing import Literal, Union, List import numpy as _np import pandas as _pd @@ -78,8 +78,16 @@ def mapparm(old, new): return off, scl -def _process_sp3_block(date, data, widths, names): - """Process a single block of SP3 data""" +def _process_sp3_block(date: str, data: str, widths: List[int], names: List[str]) -> _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. + :param List[str] names: The names of the columns in the SP3 data block. + :return _pd.DataFrame: The processed SP3 data as a DataFrame. + """ if not data or len(data) == 0: return _pd.DataFrame() epochs_dt = _pd.to_datetime(_pd.Series(date).str.slice(2, 21).values.astype(str), format=r"%Y %m %d %H %M %S") @@ -92,10 +100,23 @@ def _process_sp3_block(date, data, widths, names): return temp_sp3 -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 +def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _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 book 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. + :return pandas.DataFrame: The SP3 data as a DataFrame. + :raise FileNotFoundError: If the SP3 file specified by sp3_path does not exist. + + :note: The SP3 file format is a standard format used for representing precise satellite ephemeris and clock data. + This function reads the SP3 file, parses the header information, and extracts the data into a DataFrame. + The DataFrame columns include PV_FLAG, PRN, x_coordinate, y_coordinate, z_coordinate, clock, and various + standard deviation values. The DataFrame is processed to convert the standard deviation values to proper units + (mm/ps) and remove unnecessary columns. If pOnly is True, only P* values are included in the DataFrame. + If nodata_to_nan is True, nodata values in the SP3 POS and CLK columns are converted to NaN. """ content = _gn_io.common.path2bytes(str(sp3_path)) header_end = content.find(b"/*") @@ -153,7 +174,6 @@ def read_sp3(sp3_path, pOnly=True, nodata_to_nan=True): 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 - # print(sp3_df.index.has_duplicates()) if pOnly or parsed_header.HEAD.loc["PV_FLAG"] == "P": sp3_df = sp3_df[sp3_df.index.get_level_values("PV_FLAG") == "P"] sp3_df.attrs["HEADER"] = parsed_header # writing header data to dataframe attributes @@ -173,7 +193,13 @@ def read_sp3(sp3_path, pOnly=True, nodata_to_nan=True): return sp3_df -def parse_sp3_header(header): +def parse_sp3_header(header: str) -> _pd.DataFrame: + """ + Parse the header of an SP3 file and extract relevant information. + + :param str header: The header string of the SP3 file. + :return pandas.DataFrame: A DataFrame containing the parsed information from the SP3 header. + """ sp3_heading = _pd.Series( data=_np.asarray(_RE_SP3_HEAD.search(header).groups() + _RE_SP3_HEAD_FDESCR.search(header).groups()).astype( str @@ -201,8 +227,14 @@ def parse_sp3_header(header): return _pd.concat([sp3_heading, sv_tbl], keys=["HEAD", "SV_INFO"], axis=0) -def getVelSpline(sp3Df): - """returns in the same units as intput, e.g. km/s (needs to be x10000 to be in cm as per sp3 standard""" +def getVelSpline(sp3Df: _pd.DataFrame) -> _pd.DataFrame: + """Returns the velocity spline of the input dataframe. + + :param DataFrame sp3Df: The input dataframe containing position data. + :return DataFrame: The dataframe containing the velocity spline. + + :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) datetime = sp3dfECI.index.get_level_values("J2000") spline = _interpolate.CubicSpline(datetime, sp3dfECI.values) @@ -210,8 +242,15 @@ def getVelSpline(sp3Df): return _pd.concat([sp3Df, _pd.concat([velDf], keys=["VELi"], axis=1)], axis=1) -def getVelPoly(sp3Df, deg=35): - """takes sp3_df, interpolates the positions for -1s and +1s and outputs velocities""" +def getVelPoly(sp3Df: _pd.Dataframe, deg: int = 35): + """ + Interpolates the positions for -1s and +1s in the sp3_df DataFrame and outputs velocities. + + :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.values y = est.values @@ -240,6 +279,12 @@ def getVelPoly(sp3Df, deg=35): def gen_sp3_header(sp3_df): + """ + Generate the header for an SP3 file based on the given DataFrame. + + :param pandas.DataFrame sp3_df: The DataFrame containing the SP3 data. + :return str: The generated SP3 header as a string. + """ sp3_j2000 = sp3_df.index.levels[0].values sp3_j2000_begin = sp3_j2000[0] @@ -299,7 +344,13 @@ def gen_sp3_header(sp3_df): def gen_sp3_content(sp3_df: _pd.DataFrame, sort_outputs: bool = False, buf: Union[None, _io.TextIOBase] = None): """ Organises, formats (including nodata values), then writes out SP3 content to a buffer if provided, or returns - it otherwise. + 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. + :return str or None: The SP3 content if `buf` is None, otherwise None. """ out_buf = buf if buf is not None else _io.StringIO() if sort_outputs: @@ -473,7 +524,14 @@ def sp3merge(sp3paths, clkpaths=None, nodata_to_nan=False): def sp3_hlm_trans(a: _pd.DataFrame, b: _pd.DataFrame) -> tuple[_pd.DataFrame, list]: - """Rotates sp3_b into sp3_a. Returns a tuple of updated sp3_b and HLM array with applied computed parameters and residuals""" + """ + Rotates sp3_b into sp3_a. + + :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. + """ 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]) return b, hlm diff --git a/tests/test_sp3.py b/tests/test_sp3.py index 5e926aa..0449990 100644 --- a/tests/test_sp3.py +++ b/tests/test_sp3.py @@ -2,6 +2,7 @@ from unittest.mock import patch, mock_open import numpy as np +import pandas as pd import gnssanalysis.gn_io.sp3 as sp3 @@ -57,3 +58,26 @@ def test_read_sp3_pv(self, mock_file): result = sp3.read_sp3("mock_path", pOnly=False) self.assertTrue(len(result) == 12) self.assertEqual((np.isnan(result[("EST", "CLK")])).sum(), 10) + + def test_sp3_clock_nodata_to_nan(self): + sp3_df = pd.DataFrame({("EST", "CLK"): [999999.999999, 123456.789, 999999.999999, 987654.321]}) + sp3.sp3_clock_nodata_to_nan(sp3_df) + expected_result = pd.DataFrame({("EST", "CLK"): [np.nan, 123456.789, np.nan, 987654.321]}) + print(expected_result) + print(sp3_df) + print(sp3_df.equals(expected_result)) + self.assertTrue(sp3_df.equals(expected_result)) + + def test_sp3_pos_nodata_to_nan(self): + sp3_df = pd.DataFrame( + {("EST", "X"): [0.0, 1.0, 0.0, 2.0], ("EST", "Y"): [0.0, 0.0, 0.0, 2.0], ("EST", "Z"): [0.0, 1.0, 0.0, 0.0]} + ) + sp3.sp3_pos_nodata_to_nan(sp3_df) + expected_result = pd.DataFrame( + { + ("EST", "X"): [np.nan, 1.0, np.nan, 2.0], + ("EST", "Y"): [np.nan, 0.0, np.nan, 2.0], + ("EST", "Z"): [np.nan, 1.0, np.nan, 0.0], + } + ) + self.assertTrue(sp3_df.equals(expected_result)) From fefbaba40fda0e4bed4a533d9c296bd7709d827c Mon Sep 17 00:00:00 2001 From: Sebastien Date: Mon, 22 Jul 2024 12:54:21 +1000 Subject: [PATCH 05/15] fix: typo in type hint --- gnssanalysis/gn_io/sp3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 4c1a40a..96be83d 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -242,7 +242,7 @@ def getVelSpline(sp3Df: _pd.DataFrame) -> _pd.DataFrame: return _pd.concat([sp3Df, _pd.concat([velDf], keys=["VELi"], axis=1)], axis=1) -def getVelPoly(sp3Df: _pd.Dataframe, deg: int = 35): +def getVelPoly(sp3Df: _pd.DataFrame, deg: int = 35): """ Interpolates the positions for -1s and +1s in the sp3_df DataFrame and outputs velocities. From e7a47476db5298d9a1cbb2f6f6e2724bf7689daa Mon Sep 17 00:00:00 2001 From: Sebastien Date: Mon, 22 Jul 2024 13:47:27 +1000 Subject: [PATCH 06/15] doc: more doc and type hints. --- gnssanalysis/gn_io/sp3.py | 65 ++++++++++++++++++++++++++------------- 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 96be83d..45db97a 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -2,7 +2,7 @@ import io as _io import os as _os import re as _re -from typing import Literal, Union, List +from typing import Literal, Union, List, Tuple import numpy as _np import pandas as _pd @@ -69,8 +69,14 @@ def sp3_clock_nodata_to_nan(sp3_df: _pd.DataFrame) -> None: sp3_df.loc[nan_mask, ("EST", "CLK")] = _np.nan -def mapparm(old, new): - """scipy function f map values""" +def mapparm(old: Tuple[float, float], new: Tuple[float, float]) -> Tuple[float, float]: + """ + Maps values from the old range to the new range. + + :param Tuple[float, float] old: The old range of values. + :param Tuple[float, float] new: The new range of values. + :return Tuple[float, float]: The offset and scale factor for the mapping. + """ oldlen = old[1] - old[0] newlen = new[1] - new[0] off = (old[1] * new[0] - old[0] * new[1]) / oldlen @@ -242,7 +248,7 @@ def getVelSpline(sp3Df: _pd.DataFrame) -> _pd.DataFrame: return _pd.concat([sp3Df, _pd.concat([velDf], keys=["VELi"], axis=1)], axis=1) -def getVelPoly(sp3Df: _pd.DataFrame, deg: int = 35): +def getVelPoly(sp3Df: _pd.DataFrame, deg: int = 35) -> _pd.DataFrame: """ Interpolates the positions for -1s and +1s in the sp3_df DataFrame and outputs velocities. @@ -278,7 +284,7 @@ def getVelPoly(sp3Df: _pd.DataFrame, deg: int = 35): return _pd.concat([sp3Df, vel_i], axis=1) -def gen_sp3_header(sp3_df): +def gen_sp3_header(sp3_df: _pd.DataFrame) -> str: """ Generate the header for an SP3 file based on the given DataFrame. @@ -341,7 +347,7 @@ def gen_sp3_header(sp3_df): return "".join(line1 + line2 + sats_header.tolist() + sv_orb_head.tolist() + head_c + head_fi + comment) -def gen_sp3_content(sp3_df: _pd.DataFrame, sort_outputs: bool = False, buf: Union[None, _io.TextIOBase] = None): +def gen_sp3_content(sp3_df: _pd.DataFrame, sort_outputs: bool = False, buf: Union[None, _io.TextIOBase] = None) -> str: """ Organises, formats (including nodata values), then writes out SP3 content to a buffer if provided, or returns it otherwise. @@ -490,28 +496,43 @@ def clk_std_formatter(x): return None -def write_sp3(sp3_df, path): - """sp3 writer, dataframe to sp3 file""" +def write_sp3(sp3_df: _pd.DataFrame, path: str) -> None: + """sp3 writer, dataframe to sp3 file + + :param pandas.DataFrame sp3_df: The DataFrame containing the SP3 data. + :param str path: The path to write the SP3 file to. + """ content = gen_sp3_header(sp3_df) + gen_sp3_content(sp3_df) + "EOF" with open(path, "w") as file: file.write(content) -def merge_attrs(df_list): - """Merges attributes of a list of sp3 dataframes into a single set of attributes""" - df = _pd.concat(list(map(lambda obj: obj.attrs["HEADER"], df_list)), axis=1) +def merge_attrs(df_list: List[_pd.DataFrame]) -> _pd.Series: + """Merges attributes of a list of sp3 dataframes into a single set of attributes. + :param List[pd.DataFrame] df_list: The list of sp3 dataframes. + :return pd.Series: The merged attributes. + """ + df = _pd.concat(list(map(lambda obj: obj.attrs["HEADER"], df_list)), axis=1) mask_mixed = ~_gn_aux.unique_cols(df.loc["HEAD"]) values_if_mixed = _np.asarray(["MIX", "MIX", "MIX", None, "M", None, "MIX", "P", "MIX", "d"]) head = df[0].loc["HEAD"].values head[mask_mixed] = values_if_mixed[mask_mixed] sv_info = df.loc["SV_INFO"].max(axis=1).values.astype(int) - return _pd.Series(_np.concatenate([head, sv_info]), index=df.index) -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""" +def sp3merge( + sp3paths: List[str], clkpaths: Union[List[str], None] = None, nodata_to_nan: bool = False +) -> _pd.DataFrame: + """Reads in a list of sp3 files and optional list of clk files and merges them into a single sp3 file. + + :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] merged_sp3 = _pd.concat(sp3_dfs) merged_sp3.attrs["HEADER"] = merge_attrs(sp3_dfs) @@ -542,11 +563,15 @@ def diff_sp3_rac( sp3_test: _pd.DataFrame, hlm_mode: Literal[None, "ECF", "ECI"] = None, use_cubic_spline: bool = True, -): +) -> _pd.DataFrame: """ - Computes the difference between the two sp3 files in the radial, along-track and cross-track coordinates - the interpolator used for computation of the velocities can be based on cubic spline (default) or polynomial - Breaks if NaNs appear on unstack in getVelSpline function + 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. + :return: The DataFrame containing the difference in RAC coordinates. """ hlm_modes = [None, "ECF", "ECI"] if hlm_mode not in hlm_modes: @@ -576,9 +601,7 @@ def diff_sp3_rac( nd_rac = diff_eci.values[:, _np.newaxis] @ _gn_transform.eci2rac_rot(sp3_baseline_eci_vel) df_rac = _pd.DataFrame( nd_rac.reshape(-1, 3), - index=sp3_baseline.index, # Note that if the test and baseline have different SVs, this index will refer to - # data which is missing in the 'test' dataframe (and so is likely to be missing in - # the diff too). + index=sp3_baseline.index, columns=[["EST_RAC"] * 3, ["Radial", "Along-track", "Cross-track"]], ) From 0820e59994c2f56817c4d9aa8fccf37e8ef55dff Mon Sep 17 00:00:00 2001 From: Sebastien Date: Mon, 22 Jul 2024 13:59:05 +1000 Subject: [PATCH 07/15] doc: more docs. --- gnssanalysis/gn_io/sp3.py | 2 +- gnssanalysis/gn_transform.py | 32 +++++++++++++++++++++----------- 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 45db97a..ee85874 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -242,7 +242,7 @@ def getVelSpline(sp3Df: _pd.DataFrame) -> _pd.DataFrame: :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) - datetime = sp3dfECI.index.get_level_values("J2000") + datetime = sp3dfECI.index.values spline = _interpolate.CubicSpline(datetime, sp3dfECI.values) velDf = _pd.DataFrame(data=spline.derivative(1)(datetime), index=sp3dfECI.index, columns=sp3dfECI.columns).stack(1) return _pd.concat([sp3Df, _pd.concat([velDf], keys=["VELi"], axis=1)], axis=1) diff --git a/gnssanalysis/gn_transform.py b/gnssanalysis/gn_transform.py index 5b0cd97..5512b75 100644 --- a/gnssanalysis/gn_transform.py +++ b/gnssanalysis/gn_transform.py @@ -2,6 +2,7 @@ import numpy as _np import pandas as _pd +from typing import Tuple from . import gn_aux as _gn_aux from . import gn_const as _gn_const @@ -31,8 +32,16 @@ def gen_helm_aux(pt1, pt2): return A, rhs -def get_helmert7(pt1: _np.ndarray, pt2: _np.ndarray, scale_in_ppm: bool = True): - """inversion of 7 Helmert parameters between 2 sets of points. pt1@hlm -> pt2""" +def get_helmert7(pt1: _np.ndarray, pt2: _np.ndarray, scale_in_ppm: 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). + + :returns uple[np.ndarray, np.ndarray]: A tuple containing the Helmert parameters and the residuals. + + """ A, rhs = gen_helm_aux(pt1, pt2) sol = list(_np.linalg.lstsq(A, rhs, rcond=-1)) # parameters sol[0] = sol[0].flatten() * -1.0 # flattening the HLM params arr to [Tx, Ty, Tz, Rx, Ry, Rz, Scale/mu] @@ -54,21 +63,22 @@ def gen_rot_matrix(v): return mat + _np.eye(3) -def transform7(xyz_in, hlm_params, scale_in_ppm: bool = True): - """transformation of xyz vector with 7 helmert parameters. The helmert parameters array consists of the following: - Three transformations Tx, Ty, Tz usually in meters, or the coordinate units used for the computation of the hlm parameters. - Three rotations Rx, Ry and Rz in radians. - Scaling parameter in ppm - NOTE rotation values might be given in arcsec in literature and require conversion. In this case rotations need to be converted - to radians prior to use in transform7 with e.g. gn_aux.arcsec2rad function.""" +def transform7(xyz_in: _np.ndarray, hlm_params: _np.ndarray, scale_in_ppm: bool = True) -> _np.ndarray: + """ + Transformation of xyz vector with 7 helmert parameters. - assert hlm_params.size == 7, "there must be exactly seven parameters" + :param xyz_in: The input xyz vector. + :param hlm_params: The 7 helmert parameters: [Tx, Ty, Tz, Rx, Ry, Rz, Scale]. + :param scale_in_ppm: Whether the scale parameter is in parts per million (ppm). + :return: The transformed xyz vector. + """ + assert hlm_params.size == 7, "There must be exactly seven parameters" translation = hlm_params[0:3] rotation = gen_rot_matrix(hlm_params[3:6]) scale = hlm_params[6] if scale_in_ppm: - scale *= 1e-6 # scaling in ppm thus multiplied by 1e-6 + scale *= 1e-6 # Scaling in ppm, multiplied by 1e-6 xyz_out = (xyz_in @ rotation) * (1 + scale) + translation return xyz_out From 59c04f86c5a5fe28bbb2c5491a18d6a5c2493a8d Mon Sep 17 00:00:00 2001 From: Sebastien Date: Wed, 24 Jul 2024 12:14:06 +1000 Subject: [PATCH 08/15] refac: reducing size of read_sp3. --- gnssanalysis/gn_io/sp3.py | 106 ++++++++++++++++++++++---------------- 1 file changed, 63 insertions(+), 43 deletions(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index ee85874..12f9285 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -32,6 +32,31 @@ # File descriptor and clock _RE_SP3_HEAD_FDESCR = _re.compile(rb"\%c[ ]+(\w{1})[ ]+cc[ ](\w{3})") + +_SP3_DEF_PV_WIDTH = [1, 3, 14, 14, 14, 14, 1, 2, 1, 2, 1, 2, 1, 3, 1, 1, 1, 1, 1, 1] +_SP3_DEF_PV_NAME = [ + "PV_FLAG", + "PRN", + "x_coordinate", + "y_coordinate", + "z_coordinate", + "clock", + "_Space1", + "x_sdev", + "_Space2", + "y_sdev", + "_Space3", + "z_sdev", + "_Space4", + "c_sdev", + "_Space5", + "Clock_Event_Flag", + "Clock_Pred_Flag", + "_Space6", + "Maneuver_Flag", + "Orbit_Pred_Flag", + ] + # Nodata ie NaN constants for SP3 format SP3_CLOCK_NODATA_STRING = " 999999.999999" # Not used for reading, as fractional components are optional SP3_CLOCK_NODATA_NUMERIC = 999999 @@ -84,7 +109,8 @@ def mapparm(old: Tuple[float, float], new: Tuple[float, float]) -> Tuple[float, return off, scl -def _process_sp3_block(date: str, data: str, widths: List[int], names: List[str]) -> _pd.DataFrame: + +def _process_sp3_block(date: str, data: str, widths: List[int] = _SP3_DEF_PV_WIDTH, names: List[str] = _SP3_DEF_PV_NAME) -> _pd.DataFrame: """Process a single block of SP3 data. @@ -133,49 +159,10 @@ def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _ fline_b = header.find(b"%f") + 2 # TODO add to header parser fline = header[fline_b : fline_b + 24].strip().split(b" ") base_xyzc = _np.asarray([float(fline[0])] * 3 + [float(fline[1])]) # exponent base - _RE_SP3 = _re.compile(rb"^\*(.+)\n(.+).+", _re.MULTILINE) - data_blocks = _np.asarray(_RE_SP3.findall(string=content[: content.rfind(b"EOF")])) # Compile the regular expression pattern - pattern = _re.compile(r"^\*(.+)$", _re.MULTILINE) - # Split the content by the lines starting with '*' - blocks = pattern.split(content[: content.rfind(b"EOF")].decode()) - date_lines = blocks[1::2] - data_blocks = _np.asarray(blocks[2::2]) - # print(data_blocks) - widths = [1, 3, 14, 14, 14, 14, 1, 2, 1, 2, 1, 2, 1, 3, 1, 1, 1, 1, 1, 1] - names = [ - "PV_FLAG", - "PRN", - "x_coordinate", - "y_coordinate", - "z_coordinate", - "clock", - "Unused1", - "x_sdev", - "Unused2", - "y_sdev", - "Unused3", - "z_sdev", - "Unused4", - "c_sdev", - "Unused5", - "Clock_Event_Flag", - "Clock_Pred_Flag", - "Unused6", - "Maneuver_Flag", - "Orbit_Pred_Flag", - ] - name_float = ["x_coordinate", "y_coordinate", "z_coordinate", "clock", "x_sdev", "y_sdev", "z_sdev", "c_sdev"] - sp3_df = _pd.concat([_process_sp3_block(date, data, widths, names) for date, data in zip(date_lines, data_blocks)]) - sp3_df[name_float] = sp3_df[name_float].apply(_pd.to_numeric, errors="coerce") - sp3_df = sp3_df.loc[:, ~sp3_df.columns.str.startswith("Unused")] - # remove PRN and PV_FLAG columns - sp3_df = sp3_df.drop(columns=["PRN", "PV_FLAG"]) - # rename columsn x_coordinate -> [EST, X], y_coordinate -> [EST, Y] - sp3_df.columns = [ - ["EST", "EST", "EST", "EST", "STD", "STD", "STD", "STD", "a1", "a2", "a3", "a4"], - ["X", "Y", "Z", "CLK", "X", "Y", "Z", "CLK", "", "", "", ""], - ] + 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 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 @@ -198,6 +185,39 @@ def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _ sp3_df.attrs["path"] = sp3_path return sp3_df +def _reformat_df(sp3_df: _pd.DataFrame) -> _pd.DataFrame: + """ + Reformat the SP3 DataFrame for internal use + + :param pandas.DataFrame sp3_df: The DataFrame containing the SP3 data. + :return _pd.DataFrame: reformated SP3 data as a DataFrame. + """ + name_float = ["x_coordinate", "y_coordinate", "z_coordinate", "clock", "x_sdev", "y_sdev", "z_sdev", "c_sdev"] + sp3_df[name_float] = sp3_df[name_float].apply(_pd.to_numeric, errors="coerce") + sp3_df = sp3_df.loc[:, ~sp3_df.columns.str.startswith("_")] + # remove PRN and PV_FLAG columns + sp3_df = sp3_df.drop(columns=["PRN", "PV_FLAG"]) + # rename columsn x_coordinate -> [EST, X], y_coordinate -> [EST, Y] + sp3_df.columns = [ + ["EST", "EST", "EST", "EST", "STD", "STD", "STD", "STD", "a1", "a2", "a3", "a4"], + ["X", "Y", "Z", "CLK", "X", "Y", "Z", "CLK", "", "", "", ""], + ] + return sp3_df + + +def _split_sp3_content(content: bytes) -> Tuple[List[str], _np.ndarray]: + """ + Split the content of an SP3 file into date lines and data blocks. + + :param bytes content: The content of the SP3 file. + :return Tuple[List[str], _np.ndarray]: The date lines and data blocks. + """ + pattern = _re.compile(r"^\*(.+)$", _re.MULTILINE) + blocks = pattern.split(content[: content.rfind(b"EOF")].decode()) + date_lines = blocks[1::2] + data_blocks = _np.asarray(blocks[2::2]) + return date_lines,data_blocks + def parse_sp3_header(header: str) -> _pd.DataFrame: """ From a8bcf609832303fb2f66ebc8bf139950aaef8852 Mon Sep 17 00:00:00 2001 From: Sebastien Date: Wed, 24 Jul 2024 15:30:58 +1000 Subject: [PATCH 09/15] add: implementation of the requirements.txt file. --- .github/workflows/python-cicd-units.yaml | 2 +- requirements.txt | 14 ++++++++++++++ setup.py | 18 ++---------------- 3 files changed, 17 insertions(+), 17 deletions(-) create mode 100644 requirements.txt diff --git a/.github/workflows/python-cicd-units.yaml b/.github/workflows/python-cicd-units.yaml index 94de8cf..25d8419 100644 --- a/.github/workflows/python-cicd-units.yaml +++ b/.github/workflows/python-cicd-units.yaml @@ -15,7 +15,7 @@ jobs: - name: Install dependencies run: | python -m pip install --upgrade pip - pip install boto3 click hatanaka matplotlib numpy pandas plotext==4.2 plotly pymongo pytest scipy tqdm unlzw3 typing_extensions + pip install -r requirements.txt - name: Run tests run: | python -m unittest discover -s tests -v diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..156db0e --- /dev/null +++ b/requirements.txt @@ -0,0 +1,14 @@ +boto3 +click +hatanaka +matplotlib +numpy +pandas +plotext==4.2 +plotly +pymongo +pytest +scipy +tqdm +unlzw3 +typing_extensions \ No newline at end of file diff --git a/setup.py b/setup.py index 5484129..2711dac 100644 --- a/setup.py +++ b/setup.py @@ -4,6 +4,7 @@ # Read the contents of README file this_directory = Path(__file__).parent readme_text = (this_directory / "README.md").read_text() +requirements = (this_directory / "requirements.txt").read_text().splitlines() setuptools.setup( include_package_data=True, @@ -14,22 +15,7 @@ author_email="GNSSAnalysis@ga.gov.au", package_data={"gnssanalysis": ["py.typed"]}, packages=setuptools.find_packages(), - install_requires=[ - "boto3", - "click", - "hatanaka", - "matplotlib", - "numpy", - "pandas", - "plotext==4.2", - "plotly", - "pymongo", - "pytest", - "scipy", - "tqdm", - "unlzw3", - "typing_extensions", - ], + install_requires=requirements, entry_points={ "console_scripts": [ "diffutil = gnssanalysis:gn_utils.diffutil", From 1631d32f42ccd625b31de050478245a21ce3a08f Mon Sep 17 00:00:00 2001 From: Sebastien Date: Wed, 24 Jul 2024 19:41:45 +1000 Subject: [PATCH 10/15] fix: sp3 veloctity determination. --- gnssanalysis/gn_io/sp3.py | 55 +++++++++++++++++++++------------------ tests/test_sp3.py | 15 ++++++++--- 2 files changed, 41 insertions(+), 29 deletions(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 12f9285..1f08033 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -35,27 +35,27 @@ _SP3_DEF_PV_WIDTH = [1, 3, 14, 14, 14, 14, 1, 2, 1, 2, 1, 2, 1, 3, 1, 1, 1, 1, 1, 1] _SP3_DEF_PV_NAME = [ - "PV_FLAG", - "PRN", - "x_coordinate", - "y_coordinate", - "z_coordinate", - "clock", - "_Space1", - "x_sdev", - "_Space2", - "y_sdev", - "_Space3", - "z_sdev", - "_Space4", - "c_sdev", - "_Space5", - "Clock_Event_Flag", - "Clock_Pred_Flag", - "_Space6", - "Maneuver_Flag", - "Orbit_Pred_Flag", - ] + "PV_FLAG", + "PRN", + "x_coordinate", + "y_coordinate", + "z_coordinate", + "clock", + "_Space1", + "x_sdev", + "_Space2", + "y_sdev", + "_Space3", + "z_sdev", + "_Space4", + "c_sdev", + "_Space5", + "Clock_Event_Flag", + "Clock_Pred_Flag", + "_Space6", + "Maneuver_Flag", + "Orbit_Pred_Flag", +] # Nodata ie NaN constants for SP3 format SP3_CLOCK_NODATA_STRING = " 999999.999999" # Not used for reading, as fractional components are optional @@ -109,8 +109,9 @@ def mapparm(old: Tuple[float, float], new: Tuple[float, float]) -> Tuple[float, return off, scl - -def _process_sp3_block(date: str, data: str, widths: List[int] = _SP3_DEF_PV_WIDTH, names: List[str] = _SP3_DEF_PV_NAME) -> _pd.DataFrame: +def _process_sp3_block( + date: str, data: str, widths: List[int] = _SP3_DEF_PV_WIDTH, names: List[str] = _SP3_DEF_PV_NAME +) -> _pd.DataFrame: """Process a single block of SP3 data. @@ -185,6 +186,7 @@ def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _ sp3_df.attrs["path"] = sp3_path return sp3_df + def _reformat_df(sp3_df: _pd.DataFrame) -> _pd.DataFrame: """ Reformat the SP3 DataFrame for internal use @@ -216,7 +218,7 @@ def _split_sp3_content(content: bytes) -> Tuple[List[str], _np.ndarray]: blocks = pattern.split(content[: content.rfind(b"EOF")].decode()) date_lines = blocks[1::2] data_blocks = _np.asarray(blocks[2::2]) - return date_lines,data_blocks + return date_lines, data_blocks def parse_sp3_header(header: str) -> _pd.DataFrame: @@ -262,7 +264,7 @@ def getVelSpline(sp3Df: _pd.DataFrame) -> _pd.DataFrame: :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) - datetime = sp3dfECI.index.values + datetime = sp3dfECI.index.get_level_values("J2000").values spline = _interpolate.CubicSpline(datetime, sp3dfECI.values) velDf = _pd.DataFrame(data=spline.derivative(1)(datetime), index=sp3dfECI.index, columns=sp3dfECI.columns).stack(1) return _pd.concat([sp3Df, _pd.concat([velDf], keys=["VELi"], axis=1)], axis=1) @@ -278,7 +280,8 @@ def getVelPoly(sp3Df: _pd.DataFrame, deg: int = 35) -> _pd.DataFrame: """ est = sp3Df.unstack(1).EST[["X", "Y", "Z"]] - x = est.index.values + x = est.index.get_level_values("J2000").values + # x = est.index.values y = est.values off, scl = mapparm([x.min(), x.max()], [-1, 1]) # map from input scale to [-1,1] diff --git a/tests/test_sp3.py b/tests/test_sp3.py index 0449990..dcc5428 100644 --- a/tests/test_sp3.py +++ b/tests/test_sp3.py @@ -63,9 +63,6 @@ def test_sp3_clock_nodata_to_nan(self): sp3_df = pd.DataFrame({("EST", "CLK"): [999999.999999, 123456.789, 999999.999999, 987654.321]}) sp3.sp3_clock_nodata_to_nan(sp3_df) expected_result = pd.DataFrame({("EST", "CLK"): [np.nan, 123456.789, np.nan, 987654.321]}) - print(expected_result) - print(sp3_df) - print(sp3_df.equals(expected_result)) self.assertTrue(sp3_df.equals(expected_result)) def test_sp3_pos_nodata_to_nan(self): @@ -81,3 +78,15 @@ def test_sp3_pos_nodata_to_nan(self): } ) self.assertTrue(sp3_df.equals(expected_result)) + + @patch("builtins.open", new_callable=mock_open, read_data=input_data) + def test_velinterpolation(self, mock_file): + """ + Checking if the velocity interpolation works, right now there is no data to validate, the only thing done + is to check if the function runs without errors + """ + result = sp3.read_sp3("mock_path", pOnly=True) + r = sp3.getVelSpline(result) + r2 = sp3.getVelPoly(result, 3) + self.assertIsNotNone(r) + self.assertIsNotNone(r2) From fe81af2bef28727413acccf89b8c43a778d33aea Mon Sep 17 00:00:00 2001 From: Sebastien Date: Wed, 24 Jul 2024 19:58:59 +1000 Subject: [PATCH 11/15] review: answering comments --- .github/workflows/python-cicd-units.yaml | 4 ++-- gnssanalysis/gn_io/sp3.py | 19 +++++++++---------- 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/.github/workflows/python-cicd-units.yaml b/.github/workflows/python-cicd-units.yaml index 25d8419..bd371da 100644 --- a/.github/workflows/python-cicd-units.yaml +++ b/.github/workflows/python-cicd-units.yaml @@ -8,10 +8,10 @@ jobs: steps: - uses: actions/checkout@v2 - - name: Set up Python 3.10 + - name: Set up Python 3.12 uses: actions/setup-python@v2 with: - python-version: "3.10" + python-version: "3.12" - name: Install dependencies run: | python -m pip install --upgrade pip diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 1f08033..abfcac9 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -96,17 +96,17 @@ def sp3_clock_nodata_to_nan(sp3_df: _pd.DataFrame) -> None: def mapparm(old: Tuple[float, float], new: Tuple[float, float]) -> Tuple[float, float]: """ - Maps values from the old range to the new range. + Evaluate the offset and scale factor needed to map values from the old range to the new range. - :param Tuple[float, float] old: The old range of values. - :param Tuple[float, float] new: The new range of values. + :param Tuple[float, float] old: The range of values to be mapped from. + :param Tuple[float, float] new: The range of values to be mapped to. :return Tuple[float, float]: The offset and scale factor for the mapping. """ - oldlen = old[1] - old[0] - newlen = new[1] - new[0] - off = (old[1] * new[0] - old[0] * new[1]) / oldlen - scl = newlen / oldlen - return off, scl + old_range = old[1] - old[0] + new_range = new[1] - new[0] + offset = (old[1] * new[0] - old[0] * new[1]) / old_range + scale_factor = new_range / old_range + return offset, scale_factor def _process_sp3_block( @@ -139,7 +139,7 @@ def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _ :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 book nodata_to_nan: If True, converts 0.000000 (indicating nodata) to NaN in the SP3 POS column + :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. :return pandas.DataFrame: The SP3 data as a DataFrame. :raise FileNotFoundError: If the SP3 file specified by sp3_path does not exist. @@ -160,7 +160,6 @@ def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _ fline_b = header.find(b"%f") + 2 # TODO add to header parser fline = header[fline_b : fline_b + 24].strip().split(b" ") base_xyzc = _np.asarray([float(fline[0])] * 3 + [float(fline[1])]) # exponent base - # Compile the regular expression pattern 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) From 9e96872d248fdfab1baf185801bb5a2c7a426824 Mon Sep 17 00:00:00 2001 From: Sebastien Date: Thu, 25 Jul 2024 11:07:10 +1000 Subject: [PATCH 12/15] temp changes for Eugene --- gnssanalysis/gn_io/sp3.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index abfcac9..4dcc7d6 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -129,7 +129,7 @@ def _process_sp3_block( temp_sp3.set_index(dt_index, inplace=True) temp_sp3.index.name = "J2000" temp_sp3.set_index(temp_sp3.PRN.astype(str), append=True, inplace=True) - temp_sp3.set_index(temp_sp3.PV_FLAG.astype(str), append=True, inplace=True) + # temp_sp3.set_index(temp_sp3.PV_FLAG.astype(str), append=True, inplace=True) return temp_sp3 @@ -166,9 +166,10 @@ def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _ 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": - sp3_df = sp3_df[sp3_df.index.get_level_values("PV_FLAG") == "P"] + # print(sp3_df) + # if pOnly or parsed_header.HEAD.loc["PV_FLAG"] == "P": + # sp3_df = sp3_df[sp3_df.PV_FLAG == "P"] + # sp3_df.drop(columns="PV_FLAG", inplace=True) sp3_df.attrs["HEADER"] = parsed_header # writing header data to dataframe attributes sp3_df.attrs["path"] = sp3_path # Check for duplicate epochs, dedupe and log warning From 0f9fc27256a624ec9345b5fba785bbe3d162325a Mon Sep 17 00:00:00 2001 From: Sebastien Date: Thu, 25 Jul 2024 11:35:36 +1000 Subject: [PATCH 13/15] fix: WIP better temporary solution, just fail if velocities are in SP3 --- gnssanalysis/gn_io/sp3.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 4dcc7d6..57ed631 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -129,7 +129,7 @@ def _process_sp3_block( temp_sp3.set_index(dt_index, inplace=True) temp_sp3.index.name = "J2000" temp_sp3.set_index(temp_sp3.PRN.astype(str), append=True, inplace=True) - # temp_sp3.set_index(temp_sp3.PV_FLAG.astype(str), append=True, inplace=True) + temp_sp3.set_index(temp_sp3.PV_FLAG.astype(str), append=True, inplace=True) return temp_sp3 @@ -166,9 +166,11 @@ def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _ 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 - # print(sp3_df) - # if pOnly or parsed_header.HEAD.loc["PV_FLAG"] == "P": - # sp3_df = sp3_df[sp3_df.PV_FLAG == "P"] + if pOnly or parsed_header.HEAD.loc["PV_FLAG"] == "P": + sp3_df = sp3_df.loc[sp3_df.index.get_level_values("PV_FLAG") == "P"] + sp3_df.index = sp3_df.index.droplevel("PV_FLAG") + else: + raise NotImplementedError("Only P* values are currently supported.") # sp3_df.drop(columns="PV_FLAG", inplace=True) sp3_df.attrs["HEADER"] = parsed_header # writing header data to dataframe attributes sp3_df.attrs["path"] = sp3_path From f8e6493aacf3c73e095f60c64b9849d348402636 Mon Sep 17 00:00:00 2001 From: Sebastien Date: Thu, 25 Jul 2024 11:55:09 +1000 Subject: [PATCH 14/15] fix: handling of velocities in SP3 --- gnssanalysis/gn_io/sp3.py | 10 +++++++++- tests/test_sp3.py | 3 +-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 57ed631..6bdb589 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -166,11 +166,19 @@ def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _ 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 + print(sp3_df) if pOnly or parsed_header.HEAD.loc["PV_FLAG"] == "P": sp3_df = sp3_df.loc[sp3_df.index.get_level_values("PV_FLAG") == "P"] sp3_df.index = sp3_df.index.droplevel("PV_FLAG") else: - raise NotImplementedError("Only P* values are currently supported.") + position_df = sp3_df.xs("P", level="PV_FLAG") + velocity_df = sp3_df.xs("V", level="PV_FLAG") + velocity_df.columns = [ + ["EST", "EST", "EST", "EST", "STD", "STD", "STD", "STD", "a1", "a2", "a3", "a4"], + ["VX", "VY", "VZ", "VCLOCK", "VX", "VY", "VZ", "VCLOCK", "", "", "", ""], + ] + sp3_df = _pd.concat([position_df, velocity_df], axis=1) + # sp3_df.drop(columns="PV_FLAG", inplace=True) sp3_df.attrs["HEADER"] = parsed_header # writing header data to dataframe attributes sp3_df.attrs["path"] = sp3_path diff --git a/tests/test_sp3.py b/tests/test_sp3.py index dcc5428..bc131ec 100644 --- a/tests/test_sp3.py +++ b/tests/test_sp3.py @@ -56,8 +56,7 @@ def test_read_sp3_pOnly(self, mock_file): @patch("builtins.open", new_callable=mock_open, read_data=input_data) def test_read_sp3_pv(self, mock_file): result = sp3.read_sp3("mock_path", pOnly=False) - self.assertTrue(len(result) == 12) - self.assertEqual((np.isnan(result[("EST", "CLK")])).sum(), 10) + self.assertEqual(len(result), 6) def test_sp3_clock_nodata_to_nan(self): sp3_df = pd.DataFrame({("EST", "CLK"): [999999.999999, 123456.789, 999999.999999, 987654.321]}) From 1dc5be665da3ac5926a24b5605cf288db7b06ad6 Mon Sep 17 00:00:00 2001 From: Sebastien Date: Thu, 25 Jul 2024 12:28:40 +1000 Subject: [PATCH 15/15] review: answering comment. --- gnssanalysis/gn_io/sp3.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 6bdb589..608b172 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -166,7 +166,6 @@ def read_sp3(sp3_path: str, pOnly: bool = True, nodata_to_nan: bool = True) -> _ 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 - print(sp3_df) if pOnly or parsed_header.HEAD.loc["PV_FLAG"] == "P": sp3_df = sp3_df.loc[sp3_df.index.get_level_values("PV_FLAG") == "P"] sp3_df.index = sp3_df.index.droplevel("PV_FLAG")