Skip to content

Commit

Permalink
Merge pull request #34 from GeoscienceAustralia/jami/trop
Browse files Browse the repository at this point in the history 
modify to add capability to read both Ginan and Bernese troposphere .…
  • Loading branch information
@ronaldmaj
ronaldmaj authored Jul 29, 2024
2 parents bdead64 + 2ce2eb6 commit 6d54062
Show file tree
Hide file tree
Showing 4 changed files with 140 additions and 20 deletions.
8 changes: 6 additions & 2 deletions gnssanalysis/gn_io/sp3.py
View file
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
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"""
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
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
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))

0 comments on commit 6d54062

Please sign in to comment.