Moving antenna test

pyproject.toml

@@ -27,6 +27,8 @@ setuptools = "^70.0.0"
 viztracer = "^0.16.3"
 plotly = "^5.22.0"
 line-profiler = "^4.1.3"
+requests = "^2.32.3"
+numba = "^0.60.0"
 
 
 [build-system]

src/demos/main_gsdc2021.py

@@ -13,7 +13,7 @@
 @helpers.timeit
 def process_gsdc2021_dataset(path_dataset: Path):
     """
-    This script creates a local copy of gunzipped RINEX files from the GSDC2021 data set and process them with prx
+    This script creates a local copy of gunzipped RINEX files from the GSDC2021 data set and processes them with prx
 
     The filepath to the original GSDC2021 folder is specified in `path_dataset`.
     This has to be downloaded from the Kaggle website (requires login):

src/demos/rover_base_differencing/rover_base_differencing.py

@@ -0,0 +1,45 @@
+from pathlib import Path
+from prx import helpers
+from prx.converters import anything_to_rinex_3
+from prx.user import parse_prx_csv_file
+from prx.main import process as prx_process
+
+log = helpers.get_logger(__name__)
+
+
+def parse(obs_file: Path):
+    rnx_file = anything_to_rinex_3(obs_file)
+    assert rnx_file, f"Failed to convert {obs_file} to RINEX 3"
+    if not rnx_file.exists():
+        prx_process(rnx_file)
+    df, metadata = parse_prx_csv_file(rnx_file.with_suffix(".csv"))
+    # Receiver clocks are off by at most a few tens of milliseconds, so let's round for matching rover-base obs
+    df["time_of_reception_rounded"] = df.time_of_reception_in_receiver_time.dt.round(
+        "50ms"
+    )
+    df["sv"] = df["constellation"].astype(str) + df["prn"].astype(str).str.pad(
+        2, fillchar="0"
+    )
+    return df
+
+
+# Provides and example of how to compute between-receivers single
+# differences and double differences after processing obs files with prx.
+def main():
+    # We'll use two IGS stations that are just roughly a kilometer apart
+    df_rover, df_base = (
+        parse(Path(__file__).parent / "obs/TLSE00FRA_R_20241900000_15M_01S_MO.crx.gz"),
+        parse(Path(__file__).parent / "obs/TLSG00FRA_R_20241900000_15M_01S_MO.crx.gz"),
+    )
+    df = df_rover.merge(
+        df_base,
+        on=["time_of_reception_rounded", "sv", "rnx_obs_identifier"],
+        suffixes=("_rover", "_base"),
+    ).reset_index()
+    df["C_obs_m_sd"] = df.C_obs_m_rover - df.C_obs_m_base
+
+    print(df.info())
+
+
+if __name__ == "__main__":
+    main()

src/prx/helpers.py

@@ -4,8 +4,8 @@
 from pathlib import Path
 import logging
 
-from prx.util import is_rinex_3_obs_file, is_rinex_3_nav_file
 from prx.rinex_obs.parser import parse as prx_obs_parse
+from prx import rinex_nav as rinex_nav
 import xarray
 from imohash import imohash
 from prx import constants
@@ -14,7 +14,6 @@
 import subprocess
 import math
 import joblib
-import georinex
 import os
 from astropy.utils import iers
 from astropy import time as astrotime
@@ -73,8 +72,8 @@ def hash_of_file_content(file: Path, use_sampling: bool = False):
     hash_time = pd.Timestamp.now() - t0
     if hash_time > pd.Timedelta(seconds=1):
         log.info(
-            f"Hashing file content took {hash_time}, we might want want to think about partially hashing the"
-            f" file, e.g. using https://github.com/kalafut/py-imohash"
+            f"Hashing file content took {hash_time}, you might want want to think about partially hashing the"
+            f" file by passing sampling=True to this function."
         )
 
     return hash_string
@@ -105,6 +104,14 @@ def timestamp_2_timedelta(timestamp: pd.Timestamp, time_scale):
     assert False, f"Time scale {time_scale} not supported."
 
 
+def timestamp_to_mid_day(ts):
+    return (
+        pd.Timestamp(year=ts.year, month=1, day=1)
+        + pd.Timedelta(days=ts.day_of_year - 1)
+        + pd.Timedelta(hours=12)
+    )
+
+
 def timedelta_2_weeks_and_seconds(time_delta: pd.Timedelta):
     if pd.isnull(time_delta):
         return np.nan, np.nan
@@ -394,16 +401,6 @@ def obs_dataset_to_obs_dataframe(ds: xarray.Dataset):
     return flat_obs
 
 
-def parse_rinex_file(rinex_file_path: Path):
-    if is_rinex_3_obs_file(rinex_file_path):
-        return parse_rinex_obs_file(rinex_file_path)
-    elif is_rinex_3_nav_file(rinex_file_path):
-        return parse_rinex_nav_file(rinex_file_path)
-    assert (
-        False
-    ), f"File {rinex_file_path} appears to be neither RINEX 3 OBS nor NAV file."
-
-
 @timeit
 def parse_rinex_obs_file(rinex_file_path: Path):
     @disk_cache.cache(ignore=["rinex_file"])
@@ -417,52 +414,6 @@ def parse_rinex_obs_file_cached(rinex_file: Path, file_hash: str):
     )
 
 
-@timeit
-def parse_rinex_nav_file(rinex_file_path: Path):
-    @disk_cache.cache(ignore=["rinex_file"])
-    def parse_rinex_nav_file_cached(rinex_file: Path, file_hash: str):
-        log.info(f"Parsing {rinex_file} (hash {file_hash}) ...")
-        repair_with_gfzrnx(rinex_file)
-        return georinex.load(rinex_file)
-
-    return parse_rinex_nav_file_cached(
-        rinex_file_path, hash_of_file_content(rinex_file_path)
-    )
-
-
-def get_gpst_utc_leap_seconds(rinex_file: Path):
-    leap_seconds_astropy = compute_gps_utc_leap_seconds(
-        yyyy=int(rinex_file.name[12:16]), doy=int(rinex_file.name[16:19])
-    )
-
-    # sanity check if leap second information is in the header of the RNX NAV file
-    # compare astropy leap and RNX NAV leap seconds
-    header = georinex.rinexheader(rinex_file)
-    if "LEAP SECONDS" in header:
-        ls_before = header["LEAP SECONDS"][0:6].strip()
-        assert (
-            0 < len(ls_before) < 3
-        ), f"Unexpected leap seconds {ls_before} in {rinex_file}"
-
-        ls_after = header["LEAP SECONDS"][6:12].strip()
-        if ls_after == "":
-            return int(ls_before)
-        assert (
-            0 < len(ls_after) < 3
-        ), f"Unexpected leap seconds {ls_after} in {rinex_file}"
-
-        assert (
-            ls_after == ls_before
-        ), f"Leap second change announcement in {rinex_file}, this case is not tested, aborting."
-
-        leap_seconds_rnx = int(ls_before)
-        assert (
-            leap_seconds_rnx == leap_seconds_astropy
-        ), "leap second computed from astropy is different from RINEX NAV header"
-
-    return leap_seconds_astropy
-
-
 def is_sorted(iterable):
     return all(iterable[i] <= iterable[i + 1] for i in range(len(iterable) - 1))
 

src/prx/main.py

@@ -9,10 +9,11 @@
 import git
 
 from prx import atmospheric_corrections as atmo
-from prx.helpers import parse_rinex_file
+from prx.helpers import parse_rinex_obs_file
 from prx.rinex_nav import nav_file_discovery
 from prx import constants, helpers, converters, user
 from prx.rinex_nav import evaluate as rinex_evaluate
+from prx.rinex_nav.evaluate import parse_rinex_nav_file
 
 log = helpers.get_logger(__name__)
 
@@ -27,7 +28,9 @@ def write_prx_file(
     output_writers = {"jsonseq": write_json_text_sequence_file, "csv": write_csv_file}
     if output_format not in output_writers.keys():
         assert False, f"Output format {output_format} not supported,  we can do {list(output_writers.keys())}"
-    output_writers[output_format](prx_header, prx_records, file_name_without_extension)
+    return output_writers[output_format](
+        prx_header, prx_records, file_name_without_extension
+    )
 
 
 def write_json_text_sequence_file(
@@ -78,9 +81,6 @@ def write_csv_file(
     output_file = Path(
         f"{str(file_name_without_extension)}.{constants.cPrxCsvFileExtension}"
     )
-    # write header
-    with open(output_file, "w", encoding="utf-8") as file:
-        file.write(f"# {json.dumps(prx_header)}\n")
     flat_records["sat_elevation_deg"] = np.rad2deg(
         flat_records.elevation_rad.to_numpy()
     )
@@ -134,6 +134,15 @@ def write_csv_file(
     )
     # Keep only records with valid sat states
     records = records[records.sat_clock_offset_m.notna()]
+    # write header
+    prx_header["processing_time"] = str(
+        pd.Timestamp.now() - prx_header["processing_start_time"]
+    )
+    prx_header["processing_start_time"] = prx_header["processing_start_time"].strftime(
+        "%Y-%m-%d %H:%M:%S.%f3"
+    )
+    with open(output_file, "w", encoding="utf-8") as file:
+        file.write(f"# {json.dumps(prx_header)}\n")
     records.to_csv(
         path_or_buf=output_file,
         index=False,
@@ -142,6 +151,7 @@ def write_csv_file(
         date_format="%Y-%m-%d %H:%M:%S.%f",
     )
     log.info(f"Generated CSV prx file: {file}")
+    return output_file
 
 
 def build_metadata(input_files):
@@ -192,8 +202,9 @@ def check_assumptions(
     ), "Handling of observation files using time scales other than GPST not implemented yet."
 
 
-def warm_up_parser_cache(rinex_files):
-    _ = [parse_rinex_file(file) for file in rinex_files]
+def warm_up_parser_cache(obs_files, nav_files):
+    _ = [parse_rinex_obs_file(file) for file in obs_files]
+    _ = [parse_rinex_nav_file(file) for file in nav_files]
 
 
 @helpers.timeit
@@ -202,7 +213,8 @@ def build_records(
     rinex_3_ephemerides_files,
     approximate_receiver_ecef_position_m,
 ):
-    warm_up_parser_cache([rinex_3_obs_file] + rinex_3_ephemerides_files)
+    warm_up_parser_cache([rinex_3_obs_file], rinex_3_ephemerides_files)
+
     approximate_receiver_ecef_position_m = np.array(
         approximate_receiver_ecef_position_m
     )
@@ -283,31 +295,10 @@ def build_records(
         },
     )
 
-    sat_states_per_day = []
-    for file in rinex_3_ephemerides_files:
-        # get year and doy from NAV filename
-        year = int(file.name[12:16])
-        doy = int(file.name[16:19])
-        day_query = query.loc[
-            (
-                query.query_time_isagpst
-                >= pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy - 1)
-            )
-            & (
-                query.query_time_isagpst
-                < pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy)
-            )
-        ]
-        if day_query.empty:
-            continue
-        log.info(f"Computing satellite states for {year}-{doy:03d}")
-        sat_states_per_day.append(
-            rinex_evaluate.compute_parallel(
-                file,
-                day_query,
-            )
-        )
-    sat_states = pd.concat(sat_states_per_day)
+    sat_states = rinex_evaluate.compute_parallel(
+        rinex_3_ephemerides_files,
+        query,
+    )
     sat_states = sat_states.rename(
         columns={
             "sv": "satellite",
@@ -405,7 +396,7 @@ def signal_2_carrier_frequency(row):
         frequency_slot = row["frequency_slot"]
         if np.isnan(row["frequency_slot"]):
             return np.nan
-        # GLONASS satellites with both FDMA and CDMA signals have a frequency slot for FDMA signals,
+        # GLONASS satellites with both FDMA and CDMA signals have a frequency slot for FDMA signals.
         # for CDMA signals we use the common carrier frequency of those signals
         if len(carrier_freqs[constellation][frequency]) == 1:
             return carrier_freqs[constellation][frequency][1]
@@ -415,18 +406,28 @@ def signal_2_carrier_frequency(row):
     flat_obs.loc[:, "carrier_frequency_hz"] = flat_obs.apply(
         signal_2_carrier_frequency, axis=1
     )
+
+    file_2_to_year_and_day = {}
+    for file in rinex_3_ephemerides_files:
+        # get year and doy from NAV filename
+        nav_content = parse_rinex_nav_file(file)
+        year = int(nav_content["time"].dt.year.median())
+        doy = int(nav_content["time"].dt.day_of_year.median())
+        file_2_to_year_and_day[file] = (year, doy)
+    assert len(file_2_to_year_and_day.keys()) == len(rinex_3_ephemerides_files)
     # create a dictionary containing the headers of the different NAV files.
-    # The keys are the "YYYYDDD" (year and day of year) and are located at
-    # [12:19] of the file name using RINEX naming convention
+    # The keys are the "YYYYDDD" (year and day of year)
     nav_header_dict = {
-        file.name[12:19]: georinex.rinexheader(file)
+        f"{file_2_to_year_and_day[file][0]:03d}{file_2_to_year_and_day[file][1]:03d}": georinex.rinexheader(
+            file
+        )
         for file in rinex_3_ephemerides_files
     }
 
     for file in rinex_3_ephemerides_files:
         # get year and doy from NAV filename
-        year = int(file.name[12:16])
-        doy = int(file.name[16:19])
+        year = file_2_to_year_and_day[file][0]
+        doy = file_2_to_year_and_day[file][1]
 
         # Selection criteria: time of emission belonging to the day of the current NAV file
         mask = (
@@ -481,6 +482,7 @@ def signal_2_carrier_frequency(row):
 
 @helpers.timeit
 def process(observation_file_path: Path, output_format="csv"):
+    t0 = pd.Timestamp.now()
     # We expect a Path, but might get a string here:
     observation_file_path = Path(observation_file_path)
     log.info(
@@ -495,13 +497,14 @@ def process(observation_file_path: Path, output_format="csv"):
     metadata = build_metadata(
         {"obs_file": rinex_3_obs_file, "nav_file": aux_files["broadcast_ephemerides"]}
     )
+    metadata["processing_start_time"] = t0
     helpers.repair_with_gfzrnx(rinex_3_obs_file)
     records = build_records(
         rinex_3_obs_file,
         aux_files["broadcast_ephemerides"],
         metadata["approximate_receiver_ecef_position_m"],
     )
-    write_prx_file(
+    return write_prx_file(
         metadata,
         records,
         prx_file,