diff --git a/gnssanalysis/gn_datetime.py b/gnssanalysis/gn_datetime.py index b9fbc7b..dcee083 100644 --- a/gnssanalysis/gn_datetime.py +++ b/gnssanalysis/gn_datetime.py @@ -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 @@ -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 = ( @@ -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): diff --git a/gnssanalysis/gn_io/sp3.py b/gnssanalysis/gn_io/sp3.py index 2717fa3..17c19ba 100644 --- a/gnssanalysis/gn_io/sp3.py +++ b/gnssanalysis/gn_io/sp3.py @@ -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) + 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) @@ -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 @@ -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) @@ -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) + # 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 diff --git a/gnssanalysis/gn_utils.py b/gnssanalysis/gn_utils.py index 93a8989..dd7062b 100644 --- a/gnssanalysis/gn_utils.py +++ b/gnssanalysis/gn_utils.py @@ -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", ) @@ -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) @@ -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: diff --git a/requirements.txt b/requirements.txt index 156db0e..7f412e7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,6 +6,7 @@ numpy pandas plotext==4.2 plotly +pyfakefs pymongo pytest scipy diff --git a/tests/test_sp3.py b/tests/test_sp3.py index 08518a2..12c0af8 100644 --- a/tests/test_sp3.py +++ b/tests/test_sp3.py @@ -1,5 +1,6 @@ import unittest from unittest.mock import patch, mock_open +from pyfakefs.fake_filesystem_unittest import TestCase import numpy as np import pandas as pd @@ -7,7 +8,7 @@ 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 @@ -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) @@ -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) + self.assertEqual(result.attrs["HEADER"].HEAD.PV_FLAG, "P")