From 41887c95f3ec297ea82dbf980fe5f02bf375c93a Mon Sep 17 00:00:00 2001 From: Umma Zannat Date: Fri, 26 Jul 2024 00:15:13 +0000 Subject: [PATCH 1/4] modify to add capability to read both Ginan and Bernese troposphere .tro files --- gnssanalysis/gn_io/trop.py | 61 ++++++++++++++++++++++++++++---------- 1 file changed, 45 insertions(+), 16 deletions(-) diff --git a/gnssanalysis/gn_io/trop.py b/gnssanalysis/gn_io/trop.py index 172e4ee..e716e82 100644 --- a/gnssanalysis/gn_io/trop.py +++ b/gnssanalysis/gn_io/trop.py @@ -9,35 +9,64 @@ from .. import gn_io as _gn_io -def _read_tro_solution(path: str, recenter: bool = True) -> _pd.DataFrame: +def _read_tro_solution(path: str, recenter: bool = True, trop_mode="Ginan") -> _pd.DataFrame: + """For backwards compatibility""" + return read_tro_solution(path, recenter=recenter, trop_mode=trop_mode) + + +def read_tro_solution(path: str, recenter: bool = True, trop_mode="Ginan") -> _pd.DataFrame: """Parses tro snx file into a dataframe. - Enabling recenter overrides the default SOD values to 43200 s""" + Enabling recenter overrides the default SOD values to 43200 s + `trop_mode` can be 'Ginan' or 'Bernese'""" snx_bytes = _gn_io.common.path2bytes(path) - tro_estimate = _gn_io._snx_extract_blk(snx_bytes=snx_bytes, blk_name="TROP/SOLUTION", remove_header=True) + tro_estimate = _gn_io.sinex._snx_extract_blk(snx_bytes=snx_bytes, blk_name="TROP/SOLUTION", remove_header=True) if tro_estimate is None: _tqdm.write(f"bounds not found in {path}. Skipping.", end=" | ") return None tro_estimate = tro_estimate[0] # only single block is in tro so bytes only + if trop_mode == "Ginan": + product_headers = ["TGEWET", "TGNWET", "TROTOT", "TROWET"] + column_headers = ["CODE", "REF_EPOCH", "TGEWET", 3, "TGNWET", 5, "TROTOT", 7, "TROWET", 9] + column_dtypes = { + 0: "category", + 1: object, + 2: _np.float32, + 3: _np.float32, + 4: _np.float32, + 5: _np.float32, + 6: _np.float32, + 7: _np.float32, + 8: _np.float32, + 9: _np.float32, + } + elif trop_mode == "Bernese": + product_headers = ["TROTOT", "TGNTOT", "TGETOT"] + column_headers = ["CODE", "REF_EPOCH", "TROTOT", 3, "TGNTOT", 5, "TGETOT", 7] + column_dtypes = { + 0: "category", + 1: object, + 2: _np.float32, + 3: _np.float32, + 4: _np.float32, + 5: _np.float32, + 6: _np.float32, + 7: _np.float32, + } + else: + raise ValueError("trop_mode must be either Ginan or Bernese") + try: solution_df = _pd.read_csv( _BytesIO(tro_estimate), - delim_whitespace=True, + sep='\s+', comment=b"*", index_col=False, header=None, - names=["CODE", "REF_EPOCH", 2, 3, 4, 5, 6, 7], - dtype={ - 0: "category", - 1: object, - 2: _np.float_, - 3: _np.float32, - 4: _np.float32, - 5: _np.float32, - 6: _np.float32, - 7: _np.float32, - }, + names=column_headers, + dtype=column_dtypes, ) + except ValueError as _e: if _e.args[0][:33] == "could not convert string to float": _tqdm.write(f"{path} data corrupted. Skipping", end=" | ") @@ -45,5 +74,5 @@ def _read_tro_solution(path: str, recenter: bool = True) -> _pd.DataFrame: solution_df.REF_EPOCH = _gn_datetime.yydoysec2datetime(solution_df.REF_EPOCH, recenter=recenter, as_j2000=True) solution_df.set_index(["CODE", "REF_EPOCH"], inplace=True) - solution_df.columns = _pd.MultiIndex.from_product([["TROTOT", "TGNTOT", "TGETOT"], ["VAL", "STD"]]) + solution_df.columns = _pd.MultiIndex.from_product([product_headers, ["VAL", "STD"]]) return solution_df From ebba943e6adb20fb41d3a5ab04d9893ce53c6270 Mon Sep 17 00:00:00 2001 From: Umma Zannat Date: Fri, 26 Jul 2024 04:37:27 +0000 Subject: [PATCH 2/4] add unit tests and docs to .tro file reader --- gnssanalysis/gn_io/trop.py | 27 ++++++++++++++-- tests/test_io.py | 64 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 3 deletions(-) create mode 100644 tests/test_io.py diff --git a/gnssanalysis/gn_io/trop.py b/gnssanalysis/gn_io/trop.py index e716e82..35071b1 100644 --- a/gnssanalysis/gn_io/trop.py +++ b/gnssanalysis/gn_io/trop.py @@ -15,10 +15,31 @@ def _read_tro_solution(path: str, recenter: bool = True, trop_mode="Ginan") -> _ def read_tro_solution(path: str, recenter: bool = True, trop_mode="Ginan") -> _pd.DataFrame: - """Parses tro snx file into a dataframe. - Enabling recenter overrides the default SOD values to 43200 s - `trop_mode` can be 'Ginan' or 'Bernese'""" + """ + Parses tro snx file into a dataframe. + + :param path: path to the `.tro` file + :param recenter: recenter overrides day seconds value to midday + :param trop_mode: format of the tropo solution, can be 'Ginan' or 'Bernese' + + :raises ValueError: if `trop_mode` is unsupported + :returns: `pandas.DataFrame` containing the tropospheric solution section data + """ snx_bytes = _gn_io.common.path2bytes(path) + return read_tro_solution_bytes(snx_bytes, recenter=recenter, trop_mode=trop_mode) + + +def read_tro_solution_bytes(snx_bytes: bytes, recenter: bool = True, trop_mode="Ginan") -> _pd.DataFrame: + """ + Parses tro snx file into a dataframe. + + :param snx_bytes: contents of the `.tro` file + :param recenter: recenter overrides day seconds value to midday + :param trop_mode: format of the tropo solution, can be 'Ginan' or 'Bernese' + + :raises ValueError: if `trop_mode` is unsupported + :returns: `pandas.DataFrame` containing the tropospheric solution section data + """ tro_estimate = _gn_io.sinex._snx_extract_blk(snx_bytes=snx_bytes, blk_name="TROP/SOLUTION", remove_header=True) if tro_estimate is None: _tqdm.write(f"bounds not found in {path}. Skipping.", end=" | ") diff --git a/tests/test_io.py b/tests/test_io.py new file mode 100644 index 0000000..28d1906 --- /dev/null +++ b/tests/test_io.py @@ -0,0 +1,64 @@ +from gnssanalysis.gn_io.trop import read_tro_solution_bytes + +GINAN_TROP_FILE = """ +%=TRO 2.00 GAA 2024:185:11916.2 IGN 2024:185:11902 2024:185:11902 P MIX + ++FILE/REFERENCE + DESCRIPTION Geoscience Australia + OUTPUT Solution parameters + SOFTWARE Ginan PEA Version 3.0 + INPUT RINEX +-FILE/REFERENCE + ++TROP/SOLUTION +*STATION__ ____EPOCH_____ TGEWET STDDEV TGNWET STDDEV TROTOT STDDEV TROWET STDDEV + DARW 2024:185:11922 0.15 29.99 0.02 30.00 2443.98 299.88 165.57 299.88 + MAW1 2024:185:11922 0.05 30.00 -0.09 30.00 2252.43 299.96 10.66 299.96 + STR2 2024:185:11922 -0.05 30.00 0.02 29.99 2206.14 299.81 100.30 299.81 + DARW 2024:185:11942 1.13 29.96 -0.06 30.00 2456.94 299.56 176.28 299.56 + MAW1 2024:185:11942 -0.17 29.96 -1.31 29.86 2239.85 297.80 -2.77 297.80 + STR2 2024:185:11942 -0.35 30.00 0.02 29.93 2207.80 298.43 95.24 298.43 + DARW 2024:185:11962 0.57 29.90 0.12 29.99 2448.28 299.05 168.88 299.05 + MAW1 2024:185:11962 -0.15 29.90 -2.08 29.60 2235.75 293.17 -4.30 293.17 + STR2 2024:185:11962 -0.09 29.99 -1.11 29.79 2195.63 295.48 85.57 295.48 + DARW 2024:185:11982 1.10 29.88 0.14 29.99 2451.87 298.94 173.60 298.94 +-TROP/SOLUTION + +%=ENDTRO +""" + +BERNESE_TROP_FILE = """ +%=TRO 0.01 XYZ 24:197:01258 IGS 24:196:00000 24:197:00000 P MIX +*------------------------------------------------------------------------------- ++FILE/REFERENCE +*INFO_TYPE_________ INFO________________________________________________________ + DESCRIPTION My agency/institute + OUTPUT Precise-Point-Positioning solution generated by PPP BPE + CONTACT My e-mail address + SOFTWARE Bernese GNSS Software Version 5.2 +-FILE/REFERENCE ++TROP/SOLUTION +*SITE ____EPOCH___ TROTOT STDDEV TGNTOT STDDEV TGETOT STDDEV + ALIC 24:196:00000 2268.3 2.4 0.296 0.134 -1.446 0.184 + ALIC 24:196:03600 2260.9 1.4 0.355 0.127 -1.399 0.172 + ALIC 24:196:07200 2243.5 1.6 0.414 0.121 -1.353 0.161 + ALIC 24:196:10800 2247.9 1.4 0.473 0.114 -1.306 0.150 + ALIC 24:196:14400 2255.8 1.7 0.532 0.108 -1.259 0.140 + ALIC 24:196:18000 2247.6 1.4 0.591 0.103 -1.213 0.130 + ALIC 24:196:21600 2254.1 1.7 0.650 0.097 -1.166 0.120 + ALIC 24:196:25200 2255.3 1.3 0.709 0.093 -1.120 0.112 + ALIC 24:196:28800 2256.9 1.7 0.768 0.089 -1.073 0.105 + ALIC 24:196:32400 2268.1 1.9 0.827 0.085 -1.027 0.099 +-TROP/SOLUTION +%=ENDTRO +""" + + +def test_ginan_trop(): + # 10 records, 8 fields + assert read_tro_solution_bytes(GINAN_TROP_FILE.encode("utf-8"), trop_mode="Ginan").shape == (10, 8) + + +def test_bernese_trop(): + # 10 records, 6 fields + assert read_tro_solution_bytes(BERNESE_TROP_FILE.encode("utf-8"), trop_mode="Bernese").shape == (10, 6) From 90a7a6f802f582335e5235b06d93e3d0ff65c850 Mon Sep 17 00:00:00 2001 From: Sebastien Date: Fri, 26 Jul 2024 17:25:14 +1000 Subject: [PATCH 3/4] refact: using unittest package --- tests/{test_io.py => test_trop.py} | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) rename tests/{test_io.py => test_trop.py} (77%) diff --git a/tests/test_io.py b/tests/test_trop.py similarity index 77% rename from tests/test_io.py rename to tests/test_trop.py index 28d1906..b383a42 100644 --- a/tests/test_io.py +++ b/tests/test_trop.py @@ -1,3 +1,4 @@ +import unittest from gnssanalysis.gn_io.trop import read_tro_solution_bytes GINAN_TROP_FILE = """ @@ -30,13 +31,13 @@ BERNESE_TROP_FILE = """ %=TRO 0.01 XYZ 24:197:01258 IGS 24:196:00000 24:197:00000 P MIX *------------------------------------------------------------------------------- -+FILE/REFERENCE ++FILE/REFERENCE *INFO_TYPE_________ INFO________________________________________________________ - DESCRIPTION My agency/institute - OUTPUT Precise-Point-Positioning solution generated by PPP BPE - CONTACT My e-mail address - SOFTWARE Bernese GNSS Software Version 5.2 --FILE/REFERENCE + DESCRIPTION My agency/institute + OUTPUT Precise-Point-Positioning solution generated by PPP BPE + CONTACT My e-mail address + SOFTWARE Bernese GNSS Software Version 5.2 +-FILE/REFERENCE +TROP/SOLUTION *SITE ____EPOCH___ TROTOT STDDEV TGNTOT STDDEV TGETOT STDDEV ALIC 24:196:00000 2268.3 2.4 0.296 0.134 -1.446 0.184 @@ -54,11 +55,12 @@ """ -def test_ginan_trop(): - # 10 records, 8 fields - assert read_tro_solution_bytes(GINAN_TROP_FILE.encode("utf-8"), trop_mode="Ginan").shape == (10, 8) - +class TestTroposphere(unittest.TestCase): + def test_ginan_trop(self): + result = read_tro_solution_bytes(GINAN_TROP_FILE.encode("utf-8"), trop_mode="Ginan") + self.assertEqual(result.shape, (10, 8)) -def test_bernese_trop(): - # 10 records, 6 fields - assert read_tro_solution_bytes(BERNESE_TROP_FILE.encode("utf-8"), trop_mode="Bernese").shape == (10, 6) + def test_bernese_trop(self): + # 10 records, 6 fields + result = read_tro_solution_bytes(BERNESE_TROP_FILE.encode("utf-8"), trop_mode="Bernese") + self.assertEqual(result.shape, (10, 6)) From 2ce2eb627d14dbff42bb932c852f0407d82d13d2 Mon Sep 17 00:00:00 2001 From: Sebastien Date: Fri, 26 Jul 2024 17:29:35 +1000 Subject: [PATCH 4/4] fix: removing future warning in sp3 tests. --- gnssanalysis/gn_io/sp3.py | 8 ++++++-- tests/test_sp3.py | 2 +- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 608b172..897d7aa 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -275,7 +275,9 @@ def getVelSpline(sp3Df: _pd.DataFrame) -> _pd.DataFrame: sp3dfECI = sp3Df.EST.unstack(1)[["X", "Y", "Z"]] # _ecef2eci(sp3df) 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) + velDf = _pd.DataFrame(data=spline.derivative(1)(datetime), index=sp3dfECI.index, columns=sp3dfECI.columns).stack( + 1, future_stack=True + ) return _pd.concat([sp3Df, _pd.concat([velDf], keys=["VELi"], axis=1)], axis=1) @@ -309,7 +311,9 @@ def getVelPoly(sp3Df: _pd.DataFrame, deg: int = 35) -> _pd.DataFrame: res_prev = coeff.T.dot(inputs_prev) res_next = coeff.T.dot(inputs_next) - vel_i = _pd.DataFrame((((y - res_prev.T) + (res_next.T - y)) / 2), columns=est.columns, index=est.index).stack() + vel_i = _pd.DataFrame((((y - res_prev.T) + (res_next.T - y)) / 2), columns=est.columns, index=est.index).stack( + future_stack=True + ) vel_i.columns = [["VELi"] * 3] + [vel_i.columns.values.tolist()] diff --git a/tests/test_sp3.py b/tests/test_sp3.py index bc131ec..08518a2 100644 --- a/tests/test_sp3.py +++ b/tests/test_sp3.py @@ -86,6 +86,6 @@ def test_velinterpolation(self, mock_file): """ result = sp3.read_sp3("mock_path", pOnly=True) r = sp3.getVelSpline(result) - r2 = sp3.getVelPoly(result, 3) + r2 = sp3.getVelPoly(result, 2) self.assertIsNotNone(r) self.assertIsNotNone(r2)