Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

modify to add capability to read both Ginan and Bernese troposphere .… #34

Merged
merged 4 commits into from
Jul 29, 2024
Merged

modify to add capability to read both Ginan and Bernese troposphere .… #34

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
41887c9
modify to add capability to read both Ginan and Bernese troposphere .…
umma-zannat Jul 26, 2024
ebba943
add unit tests and docs to .tro file reader
umma-zannat Jul 26, 2024
90a7a6f
refact: using unittest package
seballgeyer Jul 26, 2024
2ce2eb6
fix: removing future warning in sp3 tests.
seballgeyer Jul 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 6 additions & 2 deletions gnssanalysis/gn_io/sp3.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)


Expand Down Expand Up @@ -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()]

Expand Down
84 changes: 67 additions & 17 deletions gnssanalysis/gn_io/trop.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,41 +9,91 @@
from .. import gn_io as _gn_io


def _read_tro_solution(path: str, recenter: bool = True) -> _pd.DataFrame:
"""Parses tro snx file into a dataframe.
Enabling recenter overrides the default SOD values to 43200 s"""
def _read_tro_solution(path: str, recenter: bool = True, trop_mode="Ginan") -> _pd.DataFrame:
"""For backwards compatibility"""
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Docstring should be one-line sphinx (missing the description of input parameters, errors raise, what's returned, etc.: :param :return)
https://github.com/NilsJPWerner/autoDocstring/blob/master/docs/one-line-sphinx.md

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function is up for deletion if no one uses it (as there is the read_tro_solution function below this

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.

: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)
tro_estimate = _gn_io._snx_extract_blk(snx_bytes=snx_bytes, blk_name="TROP/SOLUTION", remove_header=True)
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=" | ")
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=" | ")
return None

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
2 changes: 1 addition & 1 deletion tests/test_sp3.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
66 changes: 66 additions & 0 deletions tests/test_trop.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
import unittest
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
"""


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(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))
Loading