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

Jump to bottom

NPI-3429 - Fix mergesp3 utility #37

Merged
merged 7 commits into from
Aug 2, 2024
Merged

NPI-3429 - Fix mergesp3 utility #37

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
9bcf882
fix: resolve failure with sp3merge - diff attrs on DataFrames means p…
ronaldmaj Jul 30, 2024
e3988c6
update: merge_attrs function to output proper start datetime and n_ep…
ronaldmaj Jul 31, 2024
798e27a
add: unit test for sp3merge function
ronaldmaj Aug 1, 2024
ba606bc
chore: update numpy v2 NINF removed. [NPI-3429]
seballgeyer Aug 2, 2024
c7225df
dummy: intro time py package
ronaldmaj Aug 2, 2024
fbb2b4c
dummy: revert - remove time py package
ronaldmaj Aug 2, 2024
6a2d0c1
Update tests/test_sp3.py
ronaldmaj Aug 2, 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
28 changes: 24 additions & 4 deletions gnssanalysis/gn_datetime.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
"""Base time conversion functions"""

from datetime import datetime as _datetime
from datetime import timedelta as _timedelta
from io import StringIO as _StringIO
Expand Down Expand Up @@ -293,6 +294,27 @@ def j20002rnxdt(j2000secs: _np.ndarray) -> _np.ndarray:
return date_y + date_m + date_d + time_h + time_m + time_s


def rnxdt_to_datetime(rnxdt: str) -> _datetime:
"""
Transform str in RNX / SP3 format to datetime object

:param str rnxdt: String of the datetime in RNX / SP3 format: "YYYY MM DD HH mm ss.ssssssss"
:return _datetime: Tranformed python datetime object: equivalent of input rnxdt string
"""
return _datetime.strptime(rnxdt, "%Y %m %d %H %M %S.00000000")


def datetime_to_rnxdt(dt: _datetime) -> str:
"""
Transform datetime object to str of RNX / SP3 format

:param _datetime dt: Python datetime object to be transformed
:return str: Transformed str of RNX / SP3 format: "YYYY MM DD HH mm ss.ssssssss"
"""
zero_padded_str = dt.strftime("%Y %m %d %H %M %S.00000000") # Initially str is zero padded - removed in next line
return " ".join([(" " + block[1:]) if block[0] == "0" else block for block in zero_padded_str.split(" ")])


def strdatetime2datetime(dt_arr, as_j2000=True):
"""conversion of IONEX map headers ndarray Y M D h m s to datetime64"""
datetime = (
Expand Down Expand Up @@ -323,15 +345,13 @@ def snx_time_to_pydatetime(snx_time: str) -> _datetime:
@_overload
def round_timedelta(
delta: _timedelta, roundto: _timedelta, *, tol: float = ..., abs_tol: _Optional[_timedelta]
) -> _timedelta:
...
) -> _timedelta: ...


@_overload
def round_timedelta(
delta: _np.timedelta64, roundto: _np.timedelta64, *, tol: float = ..., abs_tol: _Optional[_np.timedelta64]
) -> _np.timedelta64:
...
) -> _np.timedelta64: ...


def round_timedelta(delta, roundto, *, tol=0.5, abs_tol=None):
Expand Down
35 changes: 27 additions & 8 deletions gnssanalysis/gn_io/sp3.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -508,9 +508,9 @@ def clk_std_formatter(x):

# POS nodata formatting
# Fill +/- infinity values with SP3 nodata value for POS columns
epoch_vals["X"].replace(to_replace=[_np.inf, _np.NINF], value=SP3_POS_NODATA_STRING, inplace=True)
epoch_vals["Y"].replace(to_replace=[_np.inf, _np.NINF], value=SP3_POS_NODATA_STRING, inplace=True)
epoch_vals["Z"].replace(to_replace=[_np.inf, _np.NINF], value=SP3_POS_NODATA_STRING, inplace=True)
epoch_vals["X"].replace(to_replace=[_np.inf, -_np.inf], value=SP3_POS_NODATA_STRING, inplace=True)
Copy link
Collaborator

Choose a reason for hiding this comment

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

This will need updating again for Pandas 3. I'll try to replace it with format() calls at the point of writing the data out.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

This will need updating again for Pandas 3. I'll try to replace it with format() calls at the point of writing the data out.

Yeah, there's a few changes that we'll need to make in order to make us Pandas 3 compliant. Running the sp3merge utility already has a whole bunch of warnings around this

epoch_vals["Y"].replace(to_replace=[_np.inf, -_np.inf], value=SP3_POS_NODATA_STRING, inplace=True)
epoch_vals["Z"].replace(to_replace=[_np.inf, -_np.inf], value=SP3_POS_NODATA_STRING, inplace=True)
# Now do the same for NaNs
epoch_vals["X"].fillna(value=SP3_POS_NODATA_STRING, inplace=True)
epoch_vals["Y"].fillna(value=SP3_POS_NODATA_STRING, inplace=True)
Expand All @@ -522,7 +522,7 @@ def clk_std_formatter(x):
# CLK nodata formatting
# Throw both +/- infinity, and NaN values to the SP3 clock nodata value.
# See https://stackoverflow.com/a/17478495
epoch_vals["CLK"].replace(to_replace=[_np.inf, _np.NINF], value=SP3_CLOCK_NODATA_STRING, inplace=True)
epoch_vals["CLK"].replace(to_replace=[_np.inf, -_np.inf], value=SP3_CLOCK_NODATA_STRING, inplace=True)
epoch_vals["CLK"].fillna(value=SP3_CLOCK_NODATA_STRING, inplace=True)

# Now invoke DataFrame to_string() to write out the values, leveraging our formatting functions for the
Expand Down Expand Up @@ -560,9 +560,25 @@ def merge_attrs(df_list: List[_pd.DataFrame]) -> _pd.Series:
"""
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"])
heads = df.loc["HEAD"].values
# Determine earliest start epoch from input files (and format as output string) - DATETIME in heads[2]:
dt_objs = [_gn_datetime.rnxdt_to_datetime(dt_str) for dt_str in heads[2]]
out_dt_str = _gn_datetime.datetime_to_rnxdt(min(dt_objs))
# Version Spec, PV Flag, AC label when mixed
version_str = "X" # To specify two different version, we use 'X'
pv_flag_str = "P" # If both P and V files merged, specify minimum spec - POS only
ac_str = "".join([pv[:2] for pv in sorted(set(heads[7]))])[
:4
] # Use all 4 chars assigned in spec - combine 2 char from each
# Assign values when mixed:
values_if_mixed = _np.asarray([version_str, pv_flag_str, out_dt_str, None, "M", None, "MIX", ac_str, "MX", "MIX"])
head = df[0].loc["HEAD"].values
head[mask_mixed] = values_if_mixed[mask_mixed]
# total_num_epochs needs to be assigned manually - length can be the same but have different epochs in each file
# Determine number of epochs combined dataframe will contain) - N_EPOCHS in heads[3]:
first_set_of_epochs = set(df_list[0].index.get_level_values("J2000"))
total_num_epochs = len(first_set_of_epochs.union(*[set(df.index.get_level_values("J2000")) for df in df_list[1:]]))
head[3] = str(total_num_epochs)
sv_info = df.loc["SV_INFO"].max(axis=1).values.astype(int)
return _pd.Series(_np.concatenate([head, sv_info]), index=df.index)

Expand All @@ -579,13 +595,16 @@ def sp3merge(
:return pd.DataFrame: The merged sp3 DataFrame.
"""
sp3_dfs = [read_sp3(sp3_file, nodata_to_nan=nodata_to_nan) for sp3_file in sp3paths]
# Create a new attrs dictionary to be used for the output DataFrame
merged_attrs = merge_attrs(sp3_dfs)
ronaldmaj marked this conversation as resolved.
Show resolved Hide resolved
# If attrs of two DataFrames are different, pd.concat will fail - set them to empty dict instead
for df in sp3_dfs:
df.attrs = {}
merged_sp3 = _pd.concat(sp3_dfs)
merged_sp3.attrs["HEADER"] = merge_attrs(sp3_dfs)

merged_sp3.attrs["HEADER"] = merged_attrs
if clkpaths is not None:
clk_dfs = [_gn_io.clk.read_clk(clk_file) for clk_file in clkpaths]
merged_sp3.EST.CLK = _pd.concat(clk_dfs).EST.AS * 1000000

return merged_sp3


Expand Down
20 changes: 18 additions & 2 deletions gnssanalysis/gn_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -309,7 +309,7 @@ def snxmap(sinexpaths, outdir):
@_click.option(
"-o",
"--output",
type=_click.Path(exists=True),
type=_click.Path(),
help="output path",
default=_os.curdir + "/merge.sp3",
)
Expand All @@ -326,6 +326,8 @@ def sp3merge(sp3paths, clkpaths, output, nodata_to_nan):
from .gn_io import sp3

_logging.info(msg=output)
if clkpaths == ():
clkpaths = None # clkpaths = None is a conditional used in sp3.sp3merge
merged_df = sp3.sp3merge(sp3paths=sp3paths, clkpaths=clkpaths, nodata_to_nan=nodata_to_nan)
sp3.write_sp3(sp3_df=merged_df, path=output)

Expand Down Expand Up @@ -556,11 +558,25 @@ def trace2mongo(trace_paths, db_name):
default=None,
show_default=True,
)
def orbq(input, output_path, format, csv_separation, json_format, nodata_to_nan, 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
"""
from gnssanalysis import gn_io, gn_aux, gn_diffaux

_logging.basicConfig(level="INFO") # seems that logging can only be configured before the first logging call
logger = _logging.getLogger()
# if verbose:
Expand Down
1 change: 1 addition & 0 deletions requirements.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ numpy
pandas
plotext==4.2
plotly
pyfakefs
pymongo
pytest
scipy
Expand Down
127 changes: 126 additions & 1 deletion tests/test_sp3.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,13 +1,14 @@
import unittest
from unittest.mock import patch, mock_open
from pyfakefs.fake_filesystem_unittest import TestCase

import numpy as np
import pandas as pd

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
input_data = b"""#dV2007 4 12 0 0 0.00000000 2 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
Expand Down Expand Up @@ -46,6 +47,103 @@
VG02 -12578.915944 -7977.396362 26581.116225 999999.999999
EOF"""

# second dataset a truncated version of file COD0OPSFIN_20242010000_01D_05M_ORB.SP3
input_data2 = b"""#dP2024 7 19 0 0 0.00000000 2 d+D IGS20 FIT AIUB
## 2323 432000.00000000 300.00000000 60510 0.0000000000000
+ 34 G01G02G03G04G05G06G07G08G09G10G11G12G13G14G15G16G17
+ G18G19G20G21G22G23G24G25G26G27G28G29G30G31G32R01R02
+ 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
++ 10 4 4 5 4 4 4 6 5 4 4 4 4 4 5 4 5
++ 4 5 4 4 4 4 4 4 4 6 4 6 4 4 5 6 5
++ 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 1.2500000 1.025000000 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
/* Center for Orbit Determination in Europe (CODE)
/* Final GNSS orbits and clocks for year-day 2024-2010
/* Middle day of a 3-day long-arc GRE solution
/* Product reference - DOI 10.48350/197025
/* PCV:IGS20 OL/AL:FES2014b NONE YN ORB:CoN CLK:CoN
* 2024 7 19 0 0 0.00000000
PG01 4510.358405 -23377.282442 -11792.723580 239.322216
PG02 -19585.529427 -8704.823858 16358.028672 -396.375750
PG03 -17580.234088 4691.573463 19141.243267 463.949579
PG04 -26617.740573 -946.516862 -357.173231 406.566363
PG05 10519.989071 11642.454057 -21602.301106 -180.212857
PG06 2740.107726 22949.477848 13050.231113 153.243042
PG07 -12152.989068 9271.906499 -21171.448381 -67.284723
PG08 -23564.648181 -7008.524027 -10695.029085 280.133395
PG09 -22670.253623 8239.927695 -11264.465256 278.657733
PG10 9580.520625 -23903.962622 5952.953525 -69.725084
PG11 12172.421489 23643.157258 -326.393541 -707.231762
PG12 17334.922032 -823.708469 19813.865911 -523.374443
PG13 16727.626248 13547.519770 -15930.060856 665.433543
PG14 -12686.532565 23113.995383 2015.864577 452.016655
PG15 24059.592741 4664.732321 -10971.570014 182.639915
PG16 -6877.230588 -15339.991251 -20789.984745 -232.880135
PG17 -10544.199531 13546.167670 20731.026605 674.717955
PG18 11395.390478 -10003.912563 -21773.096162 -637.936101
PG19 -91.221848 15052.899704 21730.850663 512.380157
PG20 4759.352016 20412.066571 -16337.185863 372.619575
PG21 -19969.439191 -14052.082376 12052.083047 108.053447
PG22 -7851.805614 23456.844489 10476.676398 -41.203734
PG23 16034.574849 -20014.311736 -6447.454610 276.983496
PG24 19957.601137 10937.343983 13593.650614 -483.663354
PG25 18503.862277 -12238.621370 14074.861759 498.640109
PG26 1475.737214 -23530.898479 -12009.992564 117.423098
PG27 -14637.611860 -12916.890302 -18487.558552 -31.428894
PG28 -1448.158650 -22094.490708 14658.613933 -336.951588
PG29 24922.894960 -3781.189132 -8612.199695 -588.839802
PG30 -4171.292661 19133.602125 -17670.162520 -355.440356
PG31 -9582.858249 -24150.657329 4174.829416 -226.651916
PG32 7577.135648 -14625.213465 21009.434182 -611.423040
PR01 -11174.717239 13919.992970 -18217.088128 91.816094
PR02 5986.296598 9182.597099 -22986.090345 -23.760098
* 2024 7 19 0 5 0.00000000
PG01 4796.934856 -23696.377197 -10979.751610 239.319708
PG02 -19881.646388 -9206.366139 15702.571850 -396.373498
PG03 -17231.990585 4028.826042 19602.838740 463.954491
PG04 -26611.194412 -1022.842880 595.880633 406.568667
PG05 10136.800208 12350.997157 -21386.843888 -180.212807
PG06 2450.718390 22528.454438 13825.080211 153.236129
PG07 -12588.175387 8544.963490 -21235.315069 -67.282747
PG08 -23153.571120 -7092.421087 -11502.447779 280.138208
PG09 -23081.189929 8177.801481 -10440.335865 278.662254
PG10 9739.714186 -24045.438533 5008.545199 -69.728000
PG11 12100.828204 23674.758365 628.417502 -707.233760
PG12 16927.894794 -156.593450 20175.354086 -523.375506
PG13 16123.821666 13506.143862 -16576.439914 665.434587
PG14 -12806.266956 23108.901725 1071.974280 452.019309
PG15 23667.445329 4810.788715 -11764.325180 182.641029
PG16 -6410.925819 -15949.796500 -20459.301164 -232.877069
PG17 -11294.404827 13477.778140 20374.821327 674.715594
PG18 11808.093049 -9281.574455 -21878.898484 -637.937253
PG19 -926.024078 14963.777562 21760.484912 512.381523
PG20 4496.768024 20968.671114 -15701.186882 372.619665
PG21 -20171.951734 -14415.667672 11259.651245 108.052842
PG22 -8113.879568 23715.467014 9635.144223 -41.205092
PG23 15968.691130 -19748.533800 -7367.767935 276.986193
PG24 20196.206628 11409.190382 12869.197736 -483.663786
PG25 18296.647253 -11700.764982 14777.021547 498.640196
PG26 1756.577155 -23886.875266 -11223.049448 117.421009
PG27 -13946.700542 -12958.004982 -18981.685935 -31.429176
PG28 -1093.607650 -21641.232687 15348.789400 -336.956264
PG29 25224.978997 -3692.145216 -7724.080127 -588.839284
PG30 -4576.427167 18562.149295 -18175.982398 -355.438403
PG31 -9405.141803 -24034.766821 5108.030163 -226.651713
PG32 8358.841176 -14606.846254 20714.894870 -611.422426
PR01 -11838.908273 14235.469200 -17540.792598 91.816562
PR02 5096.267067 9515.396007 -23066.803522 -23.760149
EOF
"""


class TestSp3(unittest.TestCase):
@patch("builtins.open", new_callable=mock_open, read_data=input_data)
Expand Down Expand Up @@ -89,3 +187,30 @@ def test_velinterpolation(self, mock_file):
r2 = sp3.getVelPoly(result, 2)
self.assertIsNotNone(r)
self.assertIsNotNone(r2)


class TestMergeSP3(TestCase):
def setUp(self):
self.setUpPyfakefs()

def test_sp3merge(self):
# Create some fake files
file_paths = ["/fake/dir/file1.sp3", "/fake/dir/file2.sp3"]
self.fs.create_file(file_paths[0], contents=input_data)
self.fs.create_file(file_paths[1], contents=input_data2)

# Call the function to test
result = sp3.sp3merge(sp3paths=file_paths)

# Test that epochs, satellite, attrs data is as expected:
epoch_index = result.index.get_level_values("J2000")
sat_index = result.index.get_level_values("PRN")
# Verify
self.assertEqual(min(epoch_index), 229608000)
self.assertEqual(max(epoch_index), 774619500)
self.assertEqual(sat_index[0], "G01")
self.assertEqual(sat_index[-1], "R02")
self.assertEqual(result.attrs["HEADER"].HEAD.VERSION, "d")
self.assertEqual(result.attrs["HEADER"].HEAD.AC, "AIES")
self.assertEqual(result.attrs["HEADER"].HEAD.COORD_SYS, None)
ronaldmaj marked this conversation as resolved.
Show resolved Hide resolved
self.assertEqual(result.attrs["HEADER"].HEAD.PV_FLAG, "P")
Loading