From 8ecd5776e7fcd742f4f6db596c71f67cc051766d Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Mon, 1 Jun 2026 18:32:00 +0200 Subject: [PATCH 01/26] allez --- src/prx/test/benchmark.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/prx/test/benchmark.py b/src/prx/test/benchmark.py index f5483e4..a50b7e0 100644 --- a/src/prx/test/benchmark.py +++ b/src/prx/test/benchmark.py @@ -86,7 +86,7 @@ def main(ram: bool, obs_file: Path, warm_parser_cache: bool): configure_logging("DEBUG") cases = generate_inputs( - n_steps=10, + n_steps=6, obs_file=obs_file, root=obs_file.parent / "benchmark_datasets" if obs_file is not None else None, ) From 0a5aeb6edcd882058c99613c980fad187ad3f54d Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 4 Jun 2026 18:24:54 +0200 Subject: [PATCH 02/26] add missing test --- src/prx/test/test_util.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 src/prx/test/test_util.py diff --git a/src/prx/test/test_util.py b/src/prx/test/test_util.py new file mode 100644 index 0000000..6443529 --- /dev/null +++ b/src/prx/test/test_util.py @@ -0,0 +1,19 @@ +import logging +from datetime import timedelta + +import numpy as np + +import pandas as pd +import polars as pl + +from prx.constants import cSecondsPerDay +from prx.util import timedelta_2_seconds + +log = logging.getLogger(__name__) + + +def test_timedelta_2_seconds(): + expected_timedelta_s = cSecondsPerDay + 1.23456789 + assert np.isclose(timedelta_2_seconds(pd.Timedelta(days=1, seconds=1.23456789)), expected_timedelta_s, atol=1e-9) + assert np.isclose(timedelta_2_seconds(pd.Series([pd.Timedelta(days=1, seconds=1.23456789)])).iloc[0], expected_timedelta_s, atol=1e-9) + assert np.isclose(timedelta_2_seconds(pl.Series(values=[timedelta(days=1, seconds=1.23456789)], dtype=pl.Duration(time_unit="ns")))[0], expected_timedelta_s, atol=1e-9) \ No newline at end of file From eb8bae8b6774831e58419a50f36610ef9345e7ad Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 4 Jun 2026 18:25:18 +0200 Subject: [PATCH 03/26] done --- src/prx/rinex_nav/evaluate.py | 81 +++++++++++++++++++++-------------- src/prx/util.py | 26 +++++------ 2 files changed, 59 insertions(+), 48 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index b221cc2..b500640 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -1,9 +1,11 @@ import logging from functools import lru_cache +import polars as pl import pandas as pd import numpy as np from pathlib import Path + import scipy from joblib import Parallel, delayed import georinex @@ -535,21 +537,22 @@ def compute_gal_inav_fnav_indicators(df): """ Based on RINEX 3.05, section A8 """ - df["fnav_or_inav"] = "" is_gal = df.sv.str[0] == "E" - df.loc[is_gal, "fnav_or_inav_indicator"] = np.bitwise_and( + df["fnav_or_inav_int"] = -1 + df.loc[is_gal, "fnav_or_inav_indicator_int"] = np.bitwise_and( df[is_gal].DataSrc.astype(np.uint).to_numpy(), 0b111 ) # We expect only the following navigation message types for Galileo: - indicators = set(df[is_gal].fnav_or_inav_indicator.unique()) + indicators = set(df.loc[is_gal, "fnav_or_inav_indicator_int"].unique()) assert len(indicators.intersection({1, 2, 4, 5})) == len(indicators), ( f"Unexpected Galileo navigation message type: {indicators}" ) - df.loc[is_gal & (df.fnav_or_inav_indicator == 1), "fnav_or_inav"] = "inav" - df.loc[is_gal & (df.fnav_or_inav_indicator == 2), "fnav_or_inav"] = "fnav" - df.loc[is_gal & (df.fnav_or_inav_indicator == 4), "fnav_or_inav"] = "inav" - df.loc[is_gal & (df.fnav_or_inav_indicator == 5), "fnav_or_inav"] = "inav" - return df + df["fnav_or_inav"] = "" + df.loc[is_gal & (df["fnav_or_inav_indicator_int"] == 1), "fnav_or_inav"] = "inav" + df.loc[is_gal & (df["fnav_or_inav_indicator_int"] == 2), "fnav_or_inav"] = "fnav" + df.loc[is_gal & (df["fnav_or_inav_indicator_int"] == 4), "fnav_or_inav"] = "inav" + df.loc[is_gal & (df["fnav_or_inav_indicator_int"] == 5), "fnav_or_inav"] = "inav" + return df.drop(columns=["fnav_or_inav_indicator_int"]) def to_isagpst(time, timescale, gpst_utc_leapseconds): @@ -572,37 +575,49 @@ def to_isagpst(time, timescale, gpst_utc_leapseconds): @timeit -def select_ephemerides(df, query): - df = df[df.ephemeris_reference_time_isagpst.notna()] - query = query.sort_values(by="query_time_isagpst") - df = df.sort_values(by="ephemeris_reference_time_isagpst") +def select_ephemerides(pandas_df: pd.DataFrame, pandas_query: pd.DataFrame): + df = pl.from_pandas(pandas_df) + query = pl.from_pandas(pandas_query) + df = df.drop_nulls(subset=["ephemeris_reference_time_isagpst"]) # Add fnav/inav indicator to query for to select the FNAV ephemeris for E5b signals, and INAV for other signals - query["fnav_or_inav"] = "" - query.loc[ - (query.sv.str[0] == "E") & (query.signal.str[1] == "5"), "fnav_or_inav" - ] = "fnav" - query.loc[ - (query.sv.str[0] == "E") & (query.signal.str[1] != "5"), "fnav_or_inav" - ] = "inav" - query = pd.merge_asof( - query, - df, + query = query.with_columns( + pl.when( + (pl.col("sv").str.starts_with("E")) + & (pl.col("signal").str.slice(1, 1) == "5") + ) + .then(pl.lit("fnav")) + .when( + (pl.col("sv").str.starts_with("E")) + & (pl.col("signal").str.slice(1, 1) != "5") + ) + .then(pl.lit("inav")) + .otherwise(pl.lit("")) + .alias("fnav_or_inav") + ) + query = query.sort("query_time_isagpst").join_asof( + df.sort("ephemeris_reference_time_isagpst"), left_on="query_time_isagpst", right_on="ephemeris_reference_time_isagpst", by=["sv", "fnav_or_inav"], - direction="backward", + strategy="backward", ) # Compute times w.r.t. orbit and clock reference times used by downstream computations - query["query_time_wrt_ephemeris_reference_time_s"] = ( - query["query_time_isagpst"] - query["ephemeris_reference_time_isagpst"] - ).apply(util.timedelta_2_seconds) - query["query_time_wrt_clock_reference_time_s"] = ( - query["query_time_isagpst"] - query["clock_reference_time_isagpst"] - ).apply(util.timedelta_2_seconds) - query["ephemeris_valid"] = (query["query_time_isagpst"] < query["validity_end"]) & ( - query["query_time_isagpst"] > query["validity_start"] - ) - return query + query = query.with_columns( + query_time_wrt_ephemeris_reference_time_s=util.timedelta_2_seconds( + query["query_time_isagpst"] - query["ephemeris_reference_time_isagpst"] + ), + query_time_wrt_clock_reference_time_s=util.timedelta_2_seconds( + query["query_time_isagpst"] - query["clock_reference_time_isagpst"] + ), + ephemeris_valid=( + (query["query_time_isagpst"] < query["validity_end"]) + & (query["query_time_isagpst"] > query["validity_start"]) + ).cast(bool), + ) + + result = query.to_pandas() + result["ephemeris_valid"] = result["ephemeris_valid"].astype(bool) + return result def extract_health_flag_from_query(query): diff --git a/src/prx/util.py b/src/prx/util.py index 57df626..b914ef3 100644 --- a/src/prx/util.py +++ b/src/prx/util.py @@ -13,6 +13,7 @@ import georinex import joblib import numpy as np +import polars as pl import pandas as pd import xarray from imohash import imohash @@ -233,21 +234,16 @@ def week_and_seconds_2_timedelta(weeks, seconds): return pd.Timedelta(weeks * constants.cSecondsPerWeek + seconds, "seconds") -def timedelta_2_seconds(time_delta: pd.Timedelta): - if pd.isnull(time_delta): - return np.nan - assert isinstance(time_delta, pd.Timedelta), ( - "time_delta must be of type pd.Timedelta" - ) - integer_seconds = np.float64(round(time_delta.total_seconds())) - fractional_seconds = ( - np.float64( - timedelta_2_nanoseconds(time_delta) - - integer_seconds * constants.cNanoSecondsPerSecond - ) - / constants.cNanoSecondsPerSecond - ) - return integer_seconds + fractional_seconds +def timedelta_2_seconds(time_delta: pd.Timedelta | pd.Series | pl.Series): + if isinstance(time_delta, pd.Timedelta): + return timedelta_2_seconds( + pl.Series([time_delta], dtype=pl.Duration(time_unit="ns")) + )[0] + if isinstance(time_delta, pd.Series): + return timedelta_2_seconds(pl.from_pandas(time_delta)).to_pandas() + assert isinstance(time_delta, pl.Series) + assert time_delta.dtype.time_unit == "ns" + return time_delta.dt.total_nanoseconds() / constants.cNanoSecondsPerSecond def timedelta_2_nanoseconds(time_delta: pd.Timedelta): From bbad83d109a23709a2c50ef90cb79e13a867f93d Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 4 Jun 2026 18:42:48 +0200 Subject: [PATCH 04/26] was truncating to microseconds before --- src/prx/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/prx/util.py b/src/prx/util.py index b914ef3..5fc149c 100644 --- a/src/prx/util.py +++ b/src/prx/util.py @@ -237,7 +237,7 @@ def week_and_seconds_2_timedelta(weeks, seconds): def timedelta_2_seconds(time_delta: pd.Timedelta | pd.Series | pl.Series): if isinstance(time_delta, pd.Timedelta): return timedelta_2_seconds( - pl.Series([time_delta], dtype=pl.Duration(time_unit="ns")) + pl.Series([time_delta.value], dtype=pl.Duration(time_unit="ns")) )[0] if isinstance(time_delta, pd.Series): return timedelta_2_seconds(pl.from_pandas(time_delta)).to_pandas() From c205851787bdb2a3b2882c202d40f3edfeeff67b Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 4 Jun 2026 18:42:53 +0200 Subject: [PATCH 05/26] format --- src/prx/test/test_util.py | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/src/prx/test/test_util.py b/src/prx/test/test_util.py index 6443529..ebad492 100644 --- a/src/prx/test/test_util.py +++ b/src/prx/test/test_util.py @@ -13,7 +13,26 @@ def test_timedelta_2_seconds(): - expected_timedelta_s = cSecondsPerDay + 1.23456789 - assert np.isclose(timedelta_2_seconds(pd.Timedelta(days=1, seconds=1.23456789)), expected_timedelta_s, atol=1e-9) - assert np.isclose(timedelta_2_seconds(pd.Series([pd.Timedelta(days=1, seconds=1.23456789)])).iloc[0], expected_timedelta_s, atol=1e-9) - assert np.isclose(timedelta_2_seconds(pl.Series(values=[timedelta(days=1, seconds=1.23456789)], dtype=pl.Duration(time_unit="ns")))[0], expected_timedelta_s, atol=1e-9) \ No newline at end of file + expected_timedelta_s = cSecondsPerDay + 1.23456789 + assert np.isclose( + timedelta_2_seconds(pd.Timedelta(days=1, seconds=1.23456789)), + expected_timedelta_s, + atol=1e-9, + ) + assert np.isclose( + timedelta_2_seconds(pd.Series([pd.Timedelta(days=1, seconds=1.23456789)])).iloc[ + 0 + ], + expected_timedelta_s, + atol=1e-9, + ) + assert np.isclose( + timedelta_2_seconds( + pl.Series( + values=[timedelta(days=1, seconds=1.23456789)], + dtype=pl.Duration(time_unit="ns"), + ) + )[0], + expected_timedelta_s, + atol=1e-9, + ) From 7d615075a9267fba4efebbe39d59da0c73101aca Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 4 Jun 2026 23:36:58 +0200 Subject: [PATCH 06/26] runs, fixing bugs --- src/prx/rinex_nav/evaluate.py | 137 ++++++++++++------------ src/prx/rinex_nav/test/test_evaluate.py | 17 +-- src/prx/util.py | 8 +- 3 files changed, 85 insertions(+), 77 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index b500640..2c5cf3c 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -575,10 +575,8 @@ def to_isagpst(time, timescale, gpst_utc_leapseconds): @timeit -def select_ephemerides(pandas_df: pd.DataFrame, pandas_query: pd.DataFrame): - df = pl.from_pandas(pandas_df) - query = pl.from_pandas(pandas_query) - df = df.drop_nulls(subset=["ephemeris_reference_time_isagpst"]) +def select_ephemerides(ephemerides: pl.DataFrame, query: pl.DataFrame): + ephemerides = ephemerides.drop_nulls(subset=["ephemeris_reference_time_isagpst"]) # Add fnav/inav indicator to query for to select the FNAV ephemeris for E5b signals, and INAV for other signals query = query.with_columns( pl.when( @@ -595,7 +593,7 @@ def select_ephemerides(pandas_df: pd.DataFrame, pandas_query: pd.DataFrame): .alias("fnav_or_inav") ) query = query.sort("query_time_isagpst").join_asof( - df.sort("ephemeris_reference_time_isagpst"), + ephemerides.sort("ephemeris_reference_time_isagpst"), left_on="query_time_isagpst", right_on="ephemeris_reference_time_isagpst", by=["sv", "fnav_or_inav"], @@ -615,12 +613,10 @@ def select_ephemerides(pandas_df: pd.DataFrame, pandas_query: pd.DataFrame): ).cast(bool), ) - result = query.to_pandas() - result["ephemeris_valid"] = result["ephemeris_valid"].astype(bool) - return result + return query -def extract_health_flag_from_query(query): +def extract_health_flag_from_query(query: pl.DataFrame) -> pl.DataFrame: """ Extracts the health flag for each row of a query from a `query` DataFrame containing ephemeris data. @@ -637,27 +633,25 @@ def extract_health_flag_from_query(query): "C" : "SatH1" """ - query = query.copy() - query["constellation"] = query["sv"].str[0] - - query["health_flag"] = query["health"] - if "C" in query["constellation"].unique(): - query.loc[query.constellation == "C", "health_flag"] = query.loc[ - query.constellation == "C", "SatH1" - ] - - return query["health_flag"] + query = query.with_columns( + pl.when(pl.col("sv").str.starts_with("C")) + .then(pl.col("SatH1")) + .otherwise(pl.col("health")) + .alias("health_flag") + ) + return query def compute_clock_offsets(df): - df["sat_clock_offset_m"] = constants.cGpsSpeedOfLight_mps * ( - df["SVclockBias"] - + df["SVclockDrift"] * df["query_time_wrt_clock_reference_time_s"] - + df["SVclockDriftRate"] * df["query_time_wrt_clock_reference_time_s"] ** 2 - ) - df["sat_clock_drift_mps"] = constants.cGpsSpeedOfLight_mps * ( - df["SVclockDrift"] - + 2 * df["SVclockDriftRate"] * df["query_time_wrt_clock_reference_time_s"] + df = df.with_columns( + sat_clock_offset_m=constants.cGpsSpeedOfLight_mps * pl.col("SVclockBias") + + pl.col("SVclockDrift") * pl.col("query_time_wrt_clock_reference_time_s") + + pl.col("SVclockDriftRate") + * pl.col("query_time_wrt_clock_reference_time_s") ** 2, + sat_clock_drift_mps=constants.cGpsSpeedOfLight_mps * pl.col("SVclockDrift") + + 2 + * pl.col("SVclockDriftRate") + * pl.col("query_time_wrt_clock_reference_time_s"), ) return df @@ -684,47 +678,47 @@ def compute_parallel( def compute( - rinex_nav_file_path, per_signal_query, is_query_corrected_by_sat_clock_offset=False + rinex_nav_file_path, + per_signal_query_pd, + is_query_corrected_by_sat_clock_offset=False, ): + per_signal_query = pl.from_pandas(per_signal_query_pd) query_columns = per_signal_query.columns - # per_signal_query is a pd.DataFrame with the following columns - # - time_of_reception_in_receiver_time - # - observation_value - # - signal - # - sv - # - query_time_isagpst rinex_nav_file_path = Path(rinex_nav_file_path) - ephemerides = parse_rinex_nav_file(rinex_nav_file_path) + ephemerides = pl.from_pandas(parse_rinex_nav_file(rinex_nav_file_path)) # Group delays and clock offsets can be signal-specific, so we need to match ephemerides to code signals, # not only to satellites # Example: Galileo transmits E5a clock and group delay parameters in the F/NAV message, but parameters for other # signals in the I/NAV message - per_signal_query = select_ephemerides(ephemerides, per_signal_query) + per_signal_query = select_ephemerides( + ephemerides=ephemerides, query=per_signal_query + ) # compute satellite clock bias if is_query_corrected_by_sat_clock_offset: per_signal_query = compute_clock_offsets(per_signal_query) else: # compute satellite clock offset iteratively - t = per_signal_query.query_time_wrt_clock_reference_time_s + t = per_signal_query["query_time_wrt_clock_reference_time_s"] for _ in range(2): per_signal_query = compute_clock_offsets(per_signal_query) - per_signal_query.query_time_wrt_clock_reference_time_s = ( - t - per_signal_query.sat_clock_offset_m / constants.cGpsSpeedOfLight_mps + per_signal_query = per_signal_query.with_columns( + query_time_wrt_clock_reference_time_s=( + t - pl.col("sat_clock_offset_m") / constants.cGpsSpeedOfLight_mps + ) ) # Apply sat clock correction to the query time for satellite position computation - per_signal_query.query_time_wrt_ephemeris_reference_time_s -= ( - per_signal_query.sat_clock_offset_m / constants.cGpsSpeedOfLight_mps + + per_signal_query = per_signal_query.with_columns( + query_time_wrt_ephemeris_reference_time_s=pl.col( + "query_time_wrt_ephemeris_reference_time_s" + ) + - pl.col("sat_clock_offset_m") / constants.cGpsSpeedOfLight_mps ) # Compute orbital states for each (satellite,ephemeris) pair only once: - per_sat_eph_query = ( - per_signal_query.groupby(["sv", "query_time_isagpst", "ephemeris_hash"]) - .first() - .reset_index() - ) - per_sat_eph_query = per_sat_eph_query.drop( - columns=["sat_clock_offset_m", "sat_clock_drift_mps"] - ) + per_sat_eph_query = per_signal_query.unique( + subset=["sv", "query_time_isagpst", "ephemeris_hash"] + ).drop(["sat_clock_offset_m", "sat_clock_drift_mps"]) def evaluate_orbit(sub_df): orbit_type = sub_df["orbit_type"].iloc[0] @@ -741,11 +735,12 @@ def evaluate_orbit(sub_df): sub_df[["x_m", "y_m", "z_m", "dx_mps", "dy_mps", "dz_mps"]] = np.nan return sub_df - per_sat_eph_query = per_sat_eph_query.groupby("orbit_type")[ - per_sat_eph_query.columns - ].apply(evaluate_orbit) - per_sat_eph_query = per_sat_eph_query.reset_index(drop=True) - per_sat_eph_query["health_flag"] = extract_health_flag_from_query(per_sat_eph_query) + per_sat_eph_query = pl.from_pandas( + per_sat_eph_query.to_pandas() + .groupby("orbit_type")[per_sat_eph_query.columns] + .apply(evaluate_orbit) + ) + per_sat_eph_query = extract_health_flag_from_query(per_sat_eph_query) columns_to_keep = [ "sv", "sat_pos_x_m", @@ -759,28 +754,36 @@ def evaluate_orbit(sub_df): "health_flag", "relativistic_clock_effect_m", ] - per_sat_eph_query = per_sat_eph_query[columns_to_keep] + per_sat_eph_query = per_sat_eph_query.select(columns_to_keep) # Merge the computed satellite states into the larger signal-specific query dataframe - per_signal_query = per_signal_query.merge( - per_sat_eph_query, on=["sv", "query_time_isagpst", "ephemeris_hash"] + per_signal_query = per_signal_query.join( + per_sat_eph_query, on=["sv", "query_time_isagpst", "ephemeris_hash"], how="left" ) - columns_to_keep = [ + columns_to_keep += [ "sat_clock_offset_m", "sat_clock_drift_mps", - ] + columns_to_keep - per_signal_query = compute_total_group_delays(per_signal_query) + ] + per_signal_query = pl.from_pandas( + compute_total_group_delays(per_signal_query.to_pandas()) + ) if "signal" in per_signal_query.columns: - columns_to_keep = ["signal", "sat_code_bias_m"] + columns_to_keep - columns_to_keep.append("frequency_slot") + columns_to_keep += ["signal", "sat_code_bias_m"] + columns_to_keep += ["frequency_slot"] computed_columns_to_keep = [ col for col in columns_to_keep if col not in query_columns ] - per_signal_query.loc[ - ~per_signal_query.ephemeris_valid, computed_columns_to_keep - ] = np.nan - per_signal_query = per_signal_query[columns_to_keep].reset_index(drop=True) - return per_signal_query + per_signal_query = per_signal_query.with_columns( + [ + pl.when(~pl.col("ephemeris_valid")) + .then(np.nan) + .otherwise(pl.col(col)) + .alias(col) + for col in computed_columns_to_keep + ] + ) + per_signal_query = per_signal_query.select(columns_to_keep) + return per_signal_query.to_pandas() def compute_total_group_delays( diff --git a/src/prx/rinex_nav/test/test_evaluate.py b/src/prx/rinex_nav/test/test_evaluate.py index e45f7b3..0cf1671 100644 --- a/src/prx/rinex_nav/test/test_evaluate.py +++ b/src/prx/rinex_nav/test/test_evaluate.py @@ -15,6 +15,7 @@ import shutil import pytest import itertools +import polars as pl # The following thresholds are the achieved maximum difference between broadcast and # MGEX precise orbit and clock solutions seen in this test. @@ -746,8 +747,8 @@ def test_select_ephemerides(): "TransTime": [9, 9, 100, 999], } ) - ephemerides = set_time_of_validity(ephemerides) - query = pd.DataFrame( + ephemerides = pl.from_pandas(set_time_of_validity(ephemerides)) + query = pl.from_pandas(pd.DataFrame( { "sv": ["E01", "G01", "G01"], "query_time_isagpst": [ @@ -757,15 +758,15 @@ def test_select_ephemerides(): ], "signal": ["C5X", "C1C", "C1C"], } - ) + )) query_with_ephemerides = select_ephemerides(ephemerides, query) - query_with_ephemerides = query_with_ephemerides.sort_values( + query_with_ephemerides = query_with_ephemerides.sort( by=["sv", "query_time_isagpst"] - ).reset_index(drop=True) - assert query_with_ephemerides.query_time_isagpst.equals( - pd.Series([pd.Timedelta("100s"), pd.Timedelta("50s"), pd.Timedelta("90s")]) ) - assert query_with_ephemerides.ephemeris_hash.equals(pd.Series([1, 2, 2])) + assert query_with_ephemerides["query_time_isagpst"].equals( + pl.Series([pd.Timedelta("100s").value, pd.Timedelta("50s").value, pd.Timedelta("90s").value]) + ) + assert query_with_ephemerides["ephemeris_hash"].equals(pl.Series([1, 2, 2])) def test_compute_health_flag(input_for_test_2): diff --git a/src/prx/util.py b/src/prx/util.py index 5fc149c..2a35a7f 100644 --- a/src/prx/util.py +++ b/src/prx/util.py @@ -234,7 +234,9 @@ def week_and_seconds_2_timedelta(weeks, seconds): return pd.Timedelta(weeks * constants.cSecondsPerWeek + seconds, "seconds") -def timedelta_2_seconds(time_delta: pd.Timedelta | pd.Series | pl.Series): +def timedelta_2_seconds( + time_delta: pd.Timedelta | pd.Series | pl.Series, +) -> float | pd.Series | pl.Series: if isinstance(time_delta, pd.Timedelta): return timedelta_2_seconds( pl.Series([time_delta.value], dtype=pl.Duration(time_unit="ns")) @@ -243,7 +245,9 @@ def timedelta_2_seconds(time_delta: pd.Timedelta | pd.Series | pl.Series): return timedelta_2_seconds(pl.from_pandas(time_delta)).to_pandas() assert isinstance(time_delta, pl.Series) assert time_delta.dtype.time_unit == "ns" - return time_delta.dt.total_nanoseconds() / constants.cNanoSecondsPerSecond + return ( + time_delta.dt.total_nanoseconds().cast(float) / constants.cNanoSecondsPerSecond + ) def timedelta_2_nanoseconds(time_delta: pd.Timedelta): From a3d25cd7c1c608b3c0c787eb4e4f3cdcbfd743ea Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 4 Jun 2026 23:54:34 +0200 Subject: [PATCH 07/26] fix --- src/prx/rinex_nav/evaluate.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 2c5cf3c..34dc09c 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -644,14 +644,14 @@ def extract_health_flag_from_query(query: pl.DataFrame) -> pl.DataFrame: def compute_clock_offsets(df): df = df.with_columns( - sat_clock_offset_m=constants.cGpsSpeedOfLight_mps * pl.col("SVclockBias") + sat_clock_offset_m=constants.cGpsSpeedOfLight_mps * (pl.col("SVclockBias") + pl.col("SVclockDrift") * pl.col("query_time_wrt_clock_reference_time_s") + pl.col("SVclockDriftRate") - * pl.col("query_time_wrt_clock_reference_time_s") ** 2, - sat_clock_drift_mps=constants.cGpsSpeedOfLight_mps * pl.col("SVclockDrift") + * pl.col("query_time_wrt_clock_reference_time_s") ** 2), + sat_clock_drift_mps=constants.cGpsSpeedOfLight_mps * (pl.col("SVclockDrift") + 2 * pl.col("SVclockDriftRate") - * pl.col("query_time_wrt_clock_reference_time_s"), + * pl.col("query_time_wrt_clock_reference_time_s")), ) return df @@ -759,16 +759,16 @@ def evaluate_orbit(sub_df): per_signal_query = per_signal_query.join( per_sat_eph_query, on=["sv", "query_time_isagpst", "ephemeris_hash"], how="left" ) - columns_to_keep += [ + columns_to_keep = [ "sat_clock_offset_m", "sat_clock_drift_mps", - ] + ] + columns_to_keep per_signal_query = pl.from_pandas( compute_total_group_delays(per_signal_query.to_pandas()) ) if "signal" in per_signal_query.columns: - columns_to_keep += ["signal", "sat_code_bias_m"] + columns_to_keep = ["signal", "sat_code_bias_m"] + columns_to_keep columns_to_keep += ["frequency_slot"] computed_columns_to_keep = [ col for col in columns_to_keep if col not in query_columns @@ -776,7 +776,7 @@ def evaluate_orbit(sub_df): per_signal_query = per_signal_query.with_columns( [ pl.when(~pl.col("ephemeris_valid")) - .then(np.nan) + .then(None) .otherwise(pl.col(col)) .alias(col) for col in computed_columns_to_keep From d8c418142ddee02a92351ba302da108e7f3ce2c3 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Fri, 5 Jun 2026 00:26:44 +0200 Subject: [PATCH 08/26] fixed --- src/prx/main.py | 13 +++++----- src/prx/rinex_nav/evaluate.py | 33 ++++++++++++++++--------- src/prx/rinex_nav/test/test_evaluate.py | 32 +++++++++++++++--------- 3 files changed, 48 insertions(+), 30 deletions(-) diff --git a/src/prx/main.py b/src/prx/main.py index 28e3a3d..6923f67 100644 --- a/src/prx/main.py +++ b/src/prx/main.py @@ -341,15 +341,14 @@ def build_records_levels_12( day_query["query_time_isagpst"] = day_query["query_time_isagpst"].astype( "datetime64[ns]" ) - sat_states_per_day.append( - pl.from_pandas( - rinex_evaluate.compute_parallel( - file, - day_query, - joblib_backend=joblib_backend, - ) + day_sat_states = pl.from_pandas( + rinex_evaluate.compute_parallel( + file, + day_query, + joblib_backend=joblib_backend, ) ) + sat_states_per_day.append(day_sat_states) if prx_level == 1: # drop sat group delay sat_states_per_day[-1] = sat_states_per_day[-1].drop(["sat_code_bias_m"]) sat_states = pl.concat(sat_states_per_day) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 34dc09c..4804765 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -527,6 +527,7 @@ def compute_ephemeris_and_clock_offset_reference_times(group): ) df = df.reset_index(drop=True) df = compute_gal_inav_fnav_indicators(df) + df["frequency_slot"] = int(1) if "R" in df.constellation.unique(): df["frequency_slot"] = df.FreqNum.where(df.sv.str[0] == "R", 1).astype(int) df.attrs["ionospheric_corr_GPS"] = nav_ds.ionospheric_corr_GPS @@ -644,14 +645,20 @@ def extract_health_flag_from_query(query: pl.DataFrame) -> pl.DataFrame: def compute_clock_offsets(df): df = df.with_columns( - sat_clock_offset_m=constants.cGpsSpeedOfLight_mps * (pl.col("SVclockBias") - + pl.col("SVclockDrift") * pl.col("query_time_wrt_clock_reference_time_s") - + pl.col("SVclockDriftRate") - * pl.col("query_time_wrt_clock_reference_time_s") ** 2), - sat_clock_drift_mps=constants.cGpsSpeedOfLight_mps * (pl.col("SVclockDrift") - + 2 - * pl.col("SVclockDriftRate") - * pl.col("query_time_wrt_clock_reference_time_s")), + sat_clock_offset_m=constants.cGpsSpeedOfLight_mps + * ( + pl.col("SVclockBias") + + pl.col("SVclockDrift") * pl.col("query_time_wrt_clock_reference_time_s") + + pl.col("SVclockDriftRate") + * pl.col("query_time_wrt_clock_reference_time_s") ** 2 + ), + sat_clock_drift_mps=constants.cGpsSpeedOfLight_mps + * ( + pl.col("SVclockDrift") + + 2 + * pl.col("SVclockDriftRate") + * pl.col("query_time_wrt_clock_reference_time_s") + ), ) return df @@ -674,7 +681,9 @@ def compute_parallel( ) for chunk in chunks ) - return pd.concat(processed_chunks) + result = pd.concat(processed_chunks) + result["frequency_slot"] = result["frequency_slot"].astype(float) + return result def compute( @@ -757,7 +766,9 @@ def evaluate_orbit(sub_df): per_sat_eph_query = per_sat_eph_query.select(columns_to_keep) # Merge the computed satellite states into the larger signal-specific query dataframe per_signal_query = per_signal_query.join( - per_sat_eph_query, on=["sv", "query_time_isagpst", "ephemeris_hash"], how="left" + per_sat_eph_query, + on=["sv", "query_time_isagpst", "ephemeris_hash"], + how="inner", ) columns_to_keep = [ "sat_clock_offset_m", @@ -783,7 +794,7 @@ def evaluate_orbit(sub_df): ] ) per_signal_query = per_signal_query.select(columns_to_keep) - return per_signal_query.to_pandas() + return per_signal_query.sort(by=["query_time_isagpst", "sv"]).to_pandas() def compute_total_group_delays( diff --git a/src/prx/rinex_nav/test/test_evaluate.py b/src/prx/rinex_nav/test/test_evaluate.py index 0cf1671..f931f1d 100644 --- a/src/prx/rinex_nav/test/test_evaluate.py +++ b/src/prx/rinex_nav/test/test_evaluate.py @@ -748,23 +748,31 @@ def test_select_ephemerides(): } ) ephemerides = pl.from_pandas(set_time_of_validity(ephemerides)) - query = pl.from_pandas(pd.DataFrame( - { - "sv": ["E01", "G01", "G01"], - "query_time_isagpst": [ - pd.Timedelta("100s"), - pd.Timedelta("50s"), - pd.Timedelta("90s"), - ], - "signal": ["C5X", "C1C", "C1C"], - } - )) + query = pl.from_pandas( + pd.DataFrame( + { + "sv": ["E01", "G01", "G01"], + "query_time_isagpst": [ + pd.Timedelta("100s"), + pd.Timedelta("50s"), + pd.Timedelta("90s"), + ], + "signal": ["C5X", "C1C", "C1C"], + } + ) + ) query_with_ephemerides = select_ephemerides(ephemerides, query) query_with_ephemerides = query_with_ephemerides.sort( by=["sv", "query_time_isagpst"] ) assert query_with_ephemerides["query_time_isagpst"].equals( - pl.Series([pd.Timedelta("100s").value, pd.Timedelta("50s").value, pd.Timedelta("90s").value]) + pl.Series( + [ + pd.Timedelta("100s").value, + pd.Timedelta("50s").value, + pd.Timedelta("90s").value, + ] + ) ) assert query_with_ephemerides["ephemeris_hash"].equals(pl.Series([1, 2, 2])) From 49d2f0e058429cc7ce3f31f1f876ddd16c3dff30 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Mon, 15 Jun 2026 18:53:59 +0200 Subject: [PATCH 09/26] compute array size --- pyproject.toml | 1 + src/prx/test/benchmark_copying.py | 34 +++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+) create mode 100644 src/prx/test/benchmark_copying.py diff --git a/pyproject.toml b/pyproject.toml index d5d2921..48683c4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,7 @@ dependencies = [ "polars>=1.31.0", "pyarrow>=20.0.0", "hatanaka>=2.8.1", + "psutil>=7.2.2", ] [project.scripts] diff --git a/src/prx/test/benchmark_copying.py b/src/prx/test/benchmark_copying.py new file mode 100644 index 0000000..24b29fd --- /dev/null +++ b/src/prx/test/benchmark_copying.py @@ -0,0 +1,34 @@ +import os +import time +import tracemalloc + +import numpy as np +import pandas as pd +import polars as pl +from string import ascii_uppercase +import psutil + +def main(): + array_bytes = 1e9 + n_columns = 10 + np_datatype = np.float64 + n_rows = int((array_bytes / n_columns) / 8) + time.sleep(1) + + pandas_df = pd.DataFrame({column_name: np.random.random(n_rows).astype(np_datatype) for column_name in ascii_uppercase[:n_columns]}) + assert all(dtype==np_datatype for dtype in pandas_df.dtypes) + time.sleep(1) + polars_df = pl.from_pandas(pandas_df) + polars_df = polars_df * 2 + #pandas_df = polars_df.to_pandas() + + time.sleep(1) + #input(f"Press Enter to exit process ({os.getpid()}) ...") + + +if __name__ == "__main__": + tracemalloc.start() + main() + current, peak = tracemalloc.get_traced_memory() + print(f"tracemalloc current : {current / 10**6} MB") + print(f"tracemalloc peak: {peak / 10**6} MB") From 9dd29394cbfeec6eec1deda81ccfdddce0558dbc Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Mon, 15 Jun 2026 19:17:34 +0200 Subject: [PATCH 10/26] exhibits huge peak RAM --- src/prx/test/benchmark_copying.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/prx/test/benchmark_copying.py b/src/prx/test/benchmark_copying.py index 24b29fd..7f0a481 100644 --- a/src/prx/test/benchmark_copying.py +++ b/src/prx/test/benchmark_copying.py @@ -6,29 +6,32 @@ import pandas as pd import polars as pl from string import ascii_uppercase -import psutil def main(): - array_bytes = 1e9 + array_bytes = 1e6 n_columns = 10 np_datatype = np.float64 n_rows = int((array_bytes / n_columns) / 8) - time.sleep(1) + time.sleep(1) pandas_df = pd.DataFrame({column_name: np.random.random(n_rows).astype(np_datatype) for column_name in ascii_uppercase[:n_columns]}) - assert all(dtype==np_datatype for dtype in pandas_df.dtypes) + time.sleep(1) + #pandas_df = wraps_polars(pandas_df) polars_df = pl.from_pandas(pandas_df) polars_df = polars_df * 2 - #pandas_df = polars_df.to_pandas() time.sleep(1) - #input(f"Press Enter to exit process ({os.getpid()}) ...") + pandas_df = polars_df.to_pandas() + + time.sleep(1) + input(f"Press Enter to exit process ({os.getpid()}) ...") if __name__ == "__main__": tracemalloc.start() main() + time.sleep(1) current, peak = tracemalloc.get_traced_memory() print(f"tracemalloc current : {current / 10**6} MB") print(f"tracemalloc peak: {peak / 10**6} MB") From 7ee48ea8faee975d27855d1e40aac679de008216 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Mon, 15 Jun 2026 21:03:19 +0200 Subject: [PATCH 11/26] trace allocation not going through pymalloc --- src/prx/test/benchmark.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/prx/test/benchmark.py b/src/prx/test/benchmark.py index a50b7e0..627e3c6 100644 --- a/src/prx/test/benchmark.py +++ b/src/prx/test/benchmark.py @@ -46,7 +46,7 @@ def run_case(case: dict, ram: bool, warm_parser_cache: bool) -> pd.DataFrame: memray_output.unlink(missing_ok=True) if not warm_parser_cache: disk_cache.clear() - with memray.Tracker(memray_output, follow_fork=True): + with memray.Tracker(memray_output, follow_fork=True, native_traces=True) as _: # Use multithreading here, memray does not track memory allocations in child processes with # joblib's "loky" backend. This likely makes prx slower, but we only care about memory allocation here. process( From 5dd689f1859086ae4fa9fb25445fadf6609f59f2 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Mon, 15 Jun 2026 21:03:49 +0200 Subject: [PATCH 12/26] hop --- src/prx/test/benchmark_copying.py | 41 +++++++++++++++++++------------ 1 file changed, 25 insertions(+), 16 deletions(-) diff --git a/src/prx/test/benchmark_copying.py b/src/prx/test/benchmark_copying.py index 7f0a481..8e340d0 100644 --- a/src/prx/test/benchmark_copying.py +++ b/src/prx/test/benchmark_copying.py @@ -1,4 +1,3 @@ -import os import time import tracemalloc @@ -7,31 +6,41 @@ import polars as pl from string import ascii_uppercase + +def wraps_polars(df: pd.DataFrame): + pl_df = pl.from_pandas(df) + pl_df = pl_df.select(pl.all() * 2) + return pl_df.to_pandas( + split_blocks=True, + self_destruct=True, + ) + + def main(): - array_bytes = 1e6 + array_bytes = 1e9 n_columns = 10 np_datatype = np.float64 n_rows = int((array_bytes / n_columns) / 8) - time.sleep(1) - pandas_df = pd.DataFrame({column_name: np.random.random(n_rows).astype(np_datatype) for column_name in ascii_uppercase[:n_columns]}) + pandas_df = pd.DataFrame( + { + column_name: np.random.random(n_rows).astype(np_datatype) + for column_name in ascii_uppercase[:n_columns] + } + ) time.sleep(1) - #pandas_df = wraps_polars(pandas_df) - polars_df = pl.from_pandas(pandas_df) - polars_df = polars_df * 2 + pandas_df = wraps_polars(pandas_df) time.sleep(1) - pandas_df = polars_df.to_pandas() - time.sleep(1) - input(f"Press Enter to exit process ({os.getpid()}) ...") + # input(f"Press Enter to exit process ({os.getpid()}) ...") if __name__ == "__main__": - tracemalloc.start() - main() - time.sleep(1) - current, peak = tracemalloc.get_traced_memory() - print(f"tracemalloc current : {current / 10**6} MB") - print(f"tracemalloc peak: {peak / 10**6} MB") + tracemalloc.start() + main() + time.sleep(1) + current, peak = tracemalloc.get_traced_memory() + print(f"tracemalloc current : {current / 10**6} MB") + print(f"tracemalloc peak: {peak / 10**6} MB") From e24bc88b8d379dad40782942000f124c5e622509 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 21:14:06 +0200 Subject: [PATCH 13/26] compare old and new PASSED [100%]2026-06-17 21:13:31,815 [DEBUG] [prx.util]: Function handle_bds_geos took 0 days 00:00:24.540819 to run. 2026-06-17 21:13:34,538 [DEBUG] [prx.util]: Function handle_bds_geos_faster took 0 days 00:00:02.708934 to run. --- src/prx/rinex_nav/evaluate.py | 60 +++++++++++++++++++++++-- src/prx/rinex_nav/test/test_evaluate.py | 20 ++++++++- src/prx/test/benchmark.py | 2 +- 3 files changed, 75 insertions(+), 7 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index b221cc2..2c15fa1 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -314,7 +314,7 @@ def orbital_plane_to_earth_centered_cartesian(eph): eph["dZ_k"] = eph.y_k * eph.di_k * np.cos(eph.i_k) + eph.dy_k * np.sin(eph.i_k) pass - +@timeit def handle_bds_geos(eph): # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see # Beidou_ICD_B3I_v1.0, Table 5-11 @@ -323,10 +323,10 @@ def handle_bds_geos(eph): return P_GK = np.reshape(geos[["X_k", "Y_k", "Z_k"]].to_numpy(), (-1, 1)) V_GK = np.reshape(geos[["dX_k", "dY_k", "dZ_k"]].to_numpy(), (-1, 1)) - z_angles = geos.OmegaEarthIcd_rps * geos.t_k + z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] rotation_matrices = [] + x_angle = util.deg_2_rad(-5.0) for i, z_angle in enumerate(z_angles): - x_angle = util.deg_2_rad(-5.0) Rx = np.array( [ [1, 0, 0], @@ -366,6 +366,58 @@ def frozen_to_rotating_bdcs(row): geos = geos.apply(frozen_to_rotating_bdcs, axis=1) eph[eph.is_bds_geo] = geos + return eph + + +@timeit +def handle_bds_geos_faster(eph): + # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see + # Beidou_ICD_B3I_v1.0, Table 5-11 + geos = eph[eph.is_bds_geo] + if geos.empty: + return + P_GK = np.reshape(geos[["X_k", "Y_k", "Z_k"]].to_numpy(), (-1, 1)) + V_GK = np.reshape(geos[["dX_k", "dY_k", "dZ_k"]].to_numpy(), (-1, 1)) + z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] + rotation_matrices = [] + x_angle = util.deg_2_rad(-5.0) + for i, z_angle in enumerate(z_angles): + Rx = np.array( + [ + [1, 0, 0], + [0, np.cos(x_angle), np.sin(x_angle)], + [0, -np.sin(x_angle), np.cos(x_angle)], + ] + ) + Rz = np.array( + [ + [np.cos(z_angle), np.sin(z_angle), 0], + [-np.sin(z_angle), np.cos(z_angle), 0], + [0, 0, 1], + ] + ) + rotation_matrices.append(np.matmul(Rz, Rx)) + R = scipy.sparse.block_diag(rotation_matrices) + P_K = R @ P_GK + P_K = np.reshape(P_K, (-1, 3)) + geos["X_k"] = P_K[:, 0] + geos["Y_k"] = P_K[:, 1] + geos["Z_k"] = P_K[:, 2] + # Velocity in inertial frame that coincides with BDCS at this time, ie a "frozen" ECEF frame + V_K_frozen = R @ V_GK + V_K_frozen = np.reshape(V_K_frozen, (-1, 3)) + geos["dX_k"] = V_K_frozen[:, 0] + geos["dY_k"] = V_K_frozen[:, 1] + geos["dZ_k"] = V_K_frozen[:, 2] + + # Add term due to ECEFs angular velocity w.r.t. the frozen frame + # Leverage the fact that there are only BDS GEOs + assert geos["OmegaEarthIcd_rps"].nunique() == 1 + OmegaEarthIcd_rps = geos["OmegaEarthIcd_rps"].iloc[0] + geos[["dX_k", "dY_k", "dZ_k"]] += np.cross(np.array([0, 0, - OmegaEarthIcd_rps]), geos[["X_k", "Y_k", "Z_k"]].to_numpy()) + + eph[eph.is_bds_geo] = geos + return eph # Adapted from gnss_lib_py's find_sat() @@ -392,7 +444,7 @@ def kepler_orbit_position_and_velocity(eph): ) position_in_orbital_plane(eph) orbital_plane_to_earth_centered_cartesian(eph) - handle_bds_geos(eph) + eph = handle_bds_geos(eph) eph = eph.rename( columns={ "X_k": "sat_pos_x_m", diff --git a/src/prx/rinex_nav/test/test_evaluate.py b/src/prx/rinex_nav/test/test_evaluate.py index e45f7b3..8babc17 100644 --- a/src/prx/rinex_nav/test/test_evaluate.py +++ b/src/prx/rinex_nav/test/test_evaluate.py @@ -5,13 +5,13 @@ from prx.rinex_nav.evaluate import ( select_ephemerides, set_time_of_validity, - parse_rinex_nav_file, + parse_rinex_nav_file, handle_bds_geos, handle_bds_geos_faster, ) from prx.rinex_obs.parser import parse_rinex_obs_file from prx.precise_corrections.sp3 import evaluate as sp3_evaluate from prx.rinex_nav import evaluate as rinex_nav_evaluate from prx import constants, converters -from prx.util import week_and_seconds_2_timedelta +from prx.util import week_and_seconds_2_timedelta, configure_logging import shutil import pytest import itertools @@ -833,3 +833,19 @@ def test_compute_health_flag(input_for_test_2): assert (values == test[2]).all() print("done") + + +def test_benchmark_geo_rotation(): + # GIVEN the following - not physically meaningful - satellite positions and velocities + df = pd.DataFrame(np.random.rand(int(1e6), 6)) + df.columns = ["X_k", "Y_k", "Z_k", "dX_k", "dY_k", "dZ_k"] + # With roughly half of them belonging to BDS GEO satellites + df["is_bds_geo"] = df["X_k"] > 0.5 + df["t_k"] = 1.23 + df["OmegaEarthIcd_rps"] = constants.cBdsOmegaDotEarth_rps + configure_logging("DEBUG") + reference = handle_bds_geos(df.copy()) + candidate = handle_bds_geos_faster(df.copy()) + assert reference.equals(candidate) + #assert len(result_df) == len(df) + pass \ No newline at end of file diff --git a/src/prx/test/benchmark.py b/src/prx/test/benchmark.py index f5483e4..a50b7e0 100644 --- a/src/prx/test/benchmark.py +++ b/src/prx/test/benchmark.py @@ -86,7 +86,7 @@ def main(ram: bool, obs_file: Path, warm_parser_cache: bool): configure_logging("DEBUG") cases = generate_inputs( - n_steps=10, + n_steps=6, obs_file=obs_file, root=obs_file.parent / "benchmark_datasets" if obs_file is not None else None, ) From b94bc408d7989cddb5de4d932bc31413862372b8 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 21:19:47 +0200 Subject: [PATCH 14/26] cleanup --- src/prx/rinex_nav/evaluate.py | 5 ++++- src/prx/rinex_nav/test/test_evaluate.py | 20 ++------------------ 2 files changed, 6 insertions(+), 19 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 2c15fa1..9ad3a10 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -314,6 +314,7 @@ def orbital_plane_to_earth_centered_cartesian(eph): eph["dZ_k"] = eph.y_k * eph.di_k * np.cos(eph.i_k) + eph.dy_k * np.sin(eph.i_k) pass + @timeit def handle_bds_geos(eph): # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see @@ -414,7 +415,9 @@ def handle_bds_geos_faster(eph): # Leverage the fact that there are only BDS GEOs assert geos["OmegaEarthIcd_rps"].nunique() == 1 OmegaEarthIcd_rps = geos["OmegaEarthIcd_rps"].iloc[0] - geos[["dX_k", "dY_k", "dZ_k"]] += np.cross(np.array([0, 0, - OmegaEarthIcd_rps]), geos[["X_k", "Y_k", "Z_k"]].to_numpy()) + geos[["dX_k", "dY_k", "dZ_k"]] += np.cross( + np.array([0, 0, -OmegaEarthIcd_rps]), geos[["X_k", "Y_k", "Z_k"]].to_numpy() + ) eph[eph.is_bds_geo] = geos return eph diff --git a/src/prx/rinex_nav/test/test_evaluate.py b/src/prx/rinex_nav/test/test_evaluate.py index 8babc17..e45f7b3 100644 --- a/src/prx/rinex_nav/test/test_evaluate.py +++ b/src/prx/rinex_nav/test/test_evaluate.py @@ -5,13 +5,13 @@ from prx.rinex_nav.evaluate import ( select_ephemerides, set_time_of_validity, - parse_rinex_nav_file, handle_bds_geos, handle_bds_geos_faster, + parse_rinex_nav_file, ) from prx.rinex_obs.parser import parse_rinex_obs_file from prx.precise_corrections.sp3 import evaluate as sp3_evaluate from prx.rinex_nav import evaluate as rinex_nav_evaluate from prx import constants, converters -from prx.util import week_and_seconds_2_timedelta, configure_logging +from prx.util import week_and_seconds_2_timedelta import shutil import pytest import itertools @@ -833,19 +833,3 @@ def test_compute_health_flag(input_for_test_2): assert (values == test[2]).all() print("done") - - -def test_benchmark_geo_rotation(): - # GIVEN the following - not physically meaningful - satellite positions and velocities - df = pd.DataFrame(np.random.rand(int(1e6), 6)) - df.columns = ["X_k", "Y_k", "Z_k", "dX_k", "dY_k", "dZ_k"] - # With roughly half of them belonging to BDS GEO satellites - df["is_bds_geo"] = df["X_k"] > 0.5 - df["t_k"] = 1.23 - df["OmegaEarthIcd_rps"] = constants.cBdsOmegaDotEarth_rps - configure_logging("DEBUG") - reference = handle_bds_geos(df.copy()) - candidate = handle_bds_geos_faster(df.copy()) - assert reference.equals(candidate) - #assert len(result_df) == len(df) - pass \ No newline at end of file From d57a8c8f6aaad7a5aa0b2deb9af3d3bbe1a6b04c Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 21:24:33 +0200 Subject: [PATCH 15/26] fix --- src/prx/rinex_nav/evaluate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 9ad3a10..4a0f248 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -371,12 +371,12 @@ def frozen_to_rotating_bdcs(row): @timeit -def handle_bds_geos_faster(eph): +def handle_bds_geos(eph): # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see # Beidou_ICD_B3I_v1.0, Table 5-11 geos = eph[eph.is_bds_geo] if geos.empty: - return + return eph P_GK = np.reshape(geos[["X_k", "Y_k", "Z_k"]].to_numpy(), (-1, 1)) V_GK = np.reshape(geos[["dX_k", "dY_k", "dZ_k"]].to_numpy(), (-1, 1)) z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] From 2c2f36c7fb912703cf10f87d410e33c1af73b9af Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 21:26:44 +0200 Subject: [PATCH 16/26] oups --- src/prx/rinex_nav/evaluate.py | 69 ++++------------------------------- 1 file changed, 7 insertions(+), 62 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 4a0f248..4034848 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -315,61 +315,6 @@ def orbital_plane_to_earth_centered_cartesian(eph): pass -@timeit -def handle_bds_geos(eph): - # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see - # Beidou_ICD_B3I_v1.0, Table 5-11 - geos = eph[eph.is_bds_geo] - if geos.empty: - return - P_GK = np.reshape(geos[["X_k", "Y_k", "Z_k"]].to_numpy(), (-1, 1)) - V_GK = np.reshape(geos[["dX_k", "dY_k", "dZ_k"]].to_numpy(), (-1, 1)) - z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] - rotation_matrices = [] - x_angle = util.deg_2_rad(-5.0) - for i, z_angle in enumerate(z_angles): - Rx = np.array( - [ - [1, 0, 0], - [0, np.cos(x_angle), np.sin(x_angle)], - [0, -np.sin(x_angle), np.cos(x_angle)], - ] - ) - Rz = np.array( - [ - [np.cos(z_angle), np.sin(z_angle), 0], - [-np.sin(z_angle), np.cos(z_angle), 0], - [0, 0, 1], - ] - ) - rotation_matrices.append(np.matmul(Rz, Rx)) - R = scipy.sparse.block_diag(rotation_matrices) - P_K = R @ P_GK - P_K = np.reshape(P_K, (-1, 3)) - geos["X_k"] = P_K[:, 0] - geos["Y_k"] = P_K[:, 1] - geos["Z_k"] = P_K[:, 2] - # Velocity in inertial frame that coincides with BDCS at this time, ie a "frozen" ECEF frame - V_K_frozen = R @ V_GK - V_K_frozen = np.reshape(V_K_frozen, (-1, 3)) - geos["dX_k"] = V_K_frozen[:, 0] - geos["dY_k"] = V_K_frozen[:, 1] - geos["dZ_k"] = V_K_frozen[:, 2] - - # Add term due to ECEFs angular velocity w.r.t. the frozen frame - - def frozen_to_rotating_bdcs(row): - p = np.array([row["X_k"], row["Y_k"], row["Z_k"]]) - v_frozen = np.array([row["dX_k"], row["dY_k"], row["dZ_k"]]) - v_rotating = v_frozen + np.cross(np.array([0, 0, -row.OmegaEarthIcd_rps]), p) - row[["dX_k", "dY_k", "dZ_k"]] = v_rotating - return row - - geos = geos.apply(frozen_to_rotating_bdcs, axis=1) - eph[eph.is_bds_geo] = geos - return eph - - @timeit def handle_bds_geos(eph): # Do special rotation from inertial to BDCS (ECEF) frame for Beidou GEO satellites, see @@ -382,14 +327,14 @@ def handle_bds_geos(eph): z_angles = geos["OmegaEarthIcd_rps"] * geos["t_k"] rotation_matrices = [] x_angle = util.deg_2_rad(-5.0) + Rx = np.array( + [ + [1, 0, 0], + [0, np.cos(x_angle), np.sin(x_angle)], + [0, -np.sin(x_angle), np.cos(x_angle)], + ] + ) for i, z_angle in enumerate(z_angles): - Rx = np.array( - [ - [1, 0, 0], - [0, np.cos(x_angle), np.sin(x_angle)], - [0, -np.sin(x_angle), np.cos(x_angle)], - ] - ) Rz = np.array( [ [np.cos(z_angle), np.sin(z_angle), 0], From 78d13134b7ce7742bd4207e616118d41cef8e335 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 22:17:07 +0200 Subject: [PATCH 17/26] remove row-wise applies --- src/prx/precise_corrections/sp3/evaluate.py | 5 ++- src/prx/rinex_nav/evaluate.py | 37 +++++++++------------ src/prx/rinex_nav/test/test_evaluate.py | 3 +- 3 files changed, 19 insertions(+), 26 deletions(-) diff --git a/src/prx/precise_corrections/sp3/evaluate.py b/src/prx/precise_corrections/sp3/evaluate.py index bc8338b..2794de6 100644 --- a/src/prx/precise_corrections/sp3/evaluate.py +++ b/src/prx/precise_corrections/sp3/evaluate.py @@ -7,6 +7,7 @@ from scipy.interpolate import KroghInterpolator from prx import constants, util from prx.precise_corrections.antex import antex_processing as atx_processing +from prx.util import timedelta_2_seconds log = logging.getLogger(__name__) @@ -32,9 +33,7 @@ def cached_load(file_path: Path, file_hash: str): if "position" not in col and "velocity" not in col: df.rename(columns={col: col.replace("_x", "")}, inplace=True) # Convert timestamps to seconds since GPST epoch - df["gpst_s"] = (df["time"] - constants.cGpstUtcEpoch).apply( - util.timedelta_2_seconds - ) + df["gpst_s"] = timedelta_2_seconds(df["time"] - constants.cGpstUtcEpoch) df.drop(columns=["time"], inplace=True) df["sat_clock_offset_m"] = ( constants.cGpsSpeedOfLight_mps * df["clock"] diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 4034848..1487491 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -9,7 +9,7 @@ import georinex from prx import util from prx import constants -from prx.util import timeit, try_repair_with_gfzrnx +from prx.util import timeit, try_repair_with_gfzrnx, timedelta_2_seconds log = logging.getLogger(__name__) @@ -40,7 +40,9 @@ def cached_load(rinex_file_path: Path, file_hash: str): return cached_load(rinex_file, file_content_hash) -def time_scale_integer_second_offset_wrt_gpst(time_scale, utc_gpst_leap_seconds=None): +def time_scale_integer_second_offset_wrt_gpst( + time_scale: str, utc_gpst_leap_seconds=None +): if time_scale in ["GPST", "SBAST", "QZSST", "IRNSST", "GST"]: return pd.Timedelta(seconds=0) if time_scale == "BDT": @@ -552,22 +554,13 @@ def compute_gal_inav_fnav_indicators(df): return df -def to_isagpst(time, timescale, gpst_utc_leapseconds): - if (isinstance(time, pd.Timedelta) or isinstance(time, pd.Series)) and isinstance( - timescale, str - ): - return time - time_scale_integer_second_offset_wrt_gpst( - timescale, gpst_utc_leapseconds - ) - if isinstance(time, pd.Series) and isinstance(timescale, pd.Series): - return time - timescale.apply( - lambda element: time_scale_integer_second_offset_wrt_gpst( - element, gpst_utc_leapseconds - ) - ) - - assert False, ( - f"Unexpected types: time is {type(time)}, timescale is {type(timescale)}" +def to_isagpst( + time: pd.Timedelta | pd.Series, + timescale: str | pd.Series, + gpst_utc_leapseconds: int, +) -> pd.Timedelta | pd.Series: + return time - time_scale_integer_second_offset_wrt_gpst( + timescale, gpst_utc_leapseconds ) @@ -593,12 +586,12 @@ def select_ephemerides(df, query): direction="backward", ) # Compute times w.r.t. orbit and clock reference times used by downstream computations - query["query_time_wrt_ephemeris_reference_time_s"] = ( + query["query_time_wrt_ephemeris_reference_time_s"] = timedelta_2_seconds( query["query_time_isagpst"] - query["ephemeris_reference_time_isagpst"] - ).apply(util.timedelta_2_seconds) - query["query_time_wrt_clock_reference_time_s"] = ( + ) + query["query_time_wrt_clock_reference_time_s"] = timedelta_2_seconds( query["query_time_isagpst"] - query["clock_reference_time_isagpst"] - ).apply(util.timedelta_2_seconds) + ) query["ephemeris_valid"] = (query["query_time_isagpst"] < query["validity_end"]) & ( query["query_time_isagpst"] > query["validity_start"] ) diff --git a/src/prx/rinex_nav/test/test_evaluate.py b/src/prx/rinex_nav/test/test_evaluate.py index e45f7b3..9283af4 100644 --- a/src/prx/rinex_nav/test/test_evaluate.py +++ b/src/prx/rinex_nav/test/test_evaluate.py @@ -11,7 +11,7 @@ from prx.precise_corrections.sp3 import evaluate as sp3_evaluate from prx.rinex_nav import evaluate as rinex_nav_evaluate from prx import constants, converters -from prx.util import week_and_seconds_2_timedelta +from prx.util import week_and_seconds_2_timedelta, disk_cache import shutil import pytest import itertools @@ -108,6 +108,7 @@ def input_for_test_2023(tmp_path_factory): def test_parse_nav_file(input_for_test): path_to_rnx3_nav_file = converters.anything_to_rinex_3(input_for_test["nav"]) + disk_cache.clear() df = parse_rinex_nav_file(path_to_rnx3_nav_file) assert not df.empty From 8956d60755061ccd52c95f4654c85339aad41b2a Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Wed, 17 Jun 2026 22:24:58 +0200 Subject: [PATCH 18/26] test it --- src/prx/test/test_util.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 src/prx/test/test_util.py diff --git a/src/prx/test/test_util.py b/src/prx/test/test_util.py new file mode 100644 index 0000000..ebad492 --- /dev/null +++ b/src/prx/test/test_util.py @@ -0,0 +1,38 @@ +import logging +from datetime import timedelta + +import numpy as np + +import pandas as pd +import polars as pl + +from prx.constants import cSecondsPerDay +from prx.util import timedelta_2_seconds + +log = logging.getLogger(__name__) + + +def test_timedelta_2_seconds(): + expected_timedelta_s = cSecondsPerDay + 1.23456789 + assert np.isclose( + timedelta_2_seconds(pd.Timedelta(days=1, seconds=1.23456789)), + expected_timedelta_s, + atol=1e-9, + ) + assert np.isclose( + timedelta_2_seconds(pd.Series([pd.Timedelta(days=1, seconds=1.23456789)])).iloc[ + 0 + ], + expected_timedelta_s, + atol=1e-9, + ) + assert np.isclose( + timedelta_2_seconds( + pl.Series( + values=[timedelta(days=1, seconds=1.23456789)], + dtype=pl.Duration(time_unit="ns"), + ) + )[0], + expected_timedelta_s, + atol=1e-9, + ) From c65ef4aa73fd4150cb2f9933f7276c065162ee94 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 18 Jun 2026 23:02:20 +0200 Subject: [PATCH 19/26] oups --- src/prx/util.py | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/prx/util.py b/src/prx/util.py index 57df626..06f3b5e 100644 --- a/src/prx/util.py +++ b/src/prx/util.py @@ -9,7 +9,7 @@ from pathlib import Path import importlib.metadata as md import git - +import polars as pl import georinex import joblib import numpy as np @@ -233,21 +233,20 @@ def week_and_seconds_2_timedelta(weeks, seconds): return pd.Timedelta(weeks * constants.cSecondsPerWeek + seconds, "seconds") -def timedelta_2_seconds(time_delta: pd.Timedelta): - if pd.isnull(time_delta): - return np.nan - assert isinstance(time_delta, pd.Timedelta), ( - "time_delta must be of type pd.Timedelta" - ) - integer_seconds = np.float64(round(time_delta.total_seconds())) - fractional_seconds = ( - np.float64( - timedelta_2_nanoseconds(time_delta) - - integer_seconds * constants.cNanoSecondsPerSecond - ) - / constants.cNanoSecondsPerSecond +def timedelta_2_seconds( + time_delta: pd.Timedelta | pd.Series | pl.Series, +) -> float | pd.Series | pl.Series: + if isinstance(time_delta, pd.Timedelta): + return timedelta_2_seconds( + pl.Series([time_delta.value], dtype=pl.Duration(time_unit="ns")) + )[0] + if isinstance(time_delta, pd.Series): + return timedelta_2_seconds(pl.from_pandas(time_delta)).to_pandas() + assert isinstance(time_delta, pl.Series) + assert time_delta.dtype.time_unit == "ns" + return ( + time_delta.dt.total_nanoseconds().cast(float) / constants.cNanoSecondsPerSecond ) - return integer_seconds + fractional_seconds def timedelta_2_nanoseconds(time_delta: pd.Timedelta): From bbcd9f907eaadf7ec02d946525a9dd7696102506 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Fri, 19 Jun 2026 10:02:57 +0200 Subject: [PATCH 20/26] fix --- src/prx/rinex_nav/evaluate.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 1487491..11b4be8 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -41,7 +41,7 @@ def cached_load(rinex_file_path: Path, file_hash: str): def time_scale_integer_second_offset_wrt_gpst( - time_scale: str, utc_gpst_leap_seconds=None + time_scale: str, utc_gpst_leap_seconds: int = None ): if time_scale in ["GPST", "SBAST", "QZSST", "IRNSST", "GST"]: return pd.Timedelta(seconds=0) @@ -556,8 +556,8 @@ def compute_gal_inav_fnav_indicators(df): def to_isagpst( time: pd.Timedelta | pd.Series, - timescale: str | pd.Series, - gpst_utc_leapseconds: int, + timescale: str, + gpst_utc_leapseconds: int | None, ) -> pd.Timedelta | pd.Series: return time - time_scale_integer_second_offset_wrt_gpst( timescale, gpst_utc_leapseconds From 389c8bf33aba80bf9b120c991ac31ae6623aa305 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 25 Jun 2026 08:47:17 +0200 Subject: [PATCH 21/26] switch --- src/prx/atmospheric_corrections.py | 101 +++++++++++++---------------- src/prx/main.py | 9 ++- src/prx/rinex_nav/evaluate.py | 14 ++-- src/prx/test/test_helpers.py | 2 +- src/prx/test/test_main.py | 2 +- src/prx/util.py | 26 +++++--- 6 files changed, 79 insertions(+), 75 deletions(-) diff --git a/src/prx/atmospheric_corrections.py b/src/prx/atmospheric_corrections.py index 727ffdf..23a16a6 100644 --- a/src/prx/atmospheric_corrections.py +++ b/src/prx/atmospheric_corrections.py @@ -2,9 +2,9 @@ import pandas as pd import georinex import logging - +import polars as pl from numpy.typing import NDArray - +from pathlib import Path from prx.util import deg_2_rad, ecef_2_geodetic, timedelta_2_weeks_and_seconds import prx.constants as constants @@ -85,8 +85,8 @@ def compute_l1_iono_delay_klobuchar( return iono_correction_l1_m -def add_iono_column( - flat_obs, rinex_3_ephemerides_files, approximate_receiver_ecef_position_m +def compute_iono_column( + flat_obs: pl.DataFrame, rinex_3_ephemerides_files: list[Path], approximate_receiver_ecef_position_m: NDArray[np.float64] ): # 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 @@ -95,62 +95,53 @@ def add_iono_column( file.name[12:19]: georinex.rinexheader(file) for file in rinex_3_ephemerides_files } - - idx_all_days = [] - iono_all_days = [] + [latitude_user_rad, longitude_user_rad, __] = ecef_2_geodetic( + approximate_receiver_ecef_position_m + ) + flat_obs = flat_obs.with_columns( gpsa=pl.lit(None), gpsb=pl.lit(None), iono_delay=pl.lit(None)) 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]) - - # Selection criteria: time of emission belonging to the day of the current NAV file - mask = ( - ( - flat_obs.time_of_emission_isagpst - >= pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy - 1) - ) - & ( - flat_obs.time_of_emission_isagpst - < pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy) - ) - & (flat_obs.observation_type.str.startswith("C")) - ) - mask_idx = mask.loc[mask].index - idx_all_days.append(mask_idx) - if "IONOSPHERIC CORR" in nav_header_dict[f"{year:03d}" + f"{doy:03d}"]: - logging.info(f"Computing iono delay for {year}-{doy:03d}") - time_of_emission_weeksecond_isagpst = timedelta_2_weeks_and_seconds( - flat_obs.loc[mask_idx, "time_of_emission_isagpst"] - - constants.system_time_scale_rinex_utc_epoch["GPST"] - )[1].to_numpy() - [latitude_user_rad, longitude_user_rad, __] = ecef_2_geodetic( - approximate_receiver_ecef_position_m - ) - iono_all_days.append( - compute_l1_iono_delay_klobuchar( - time_of_emission_weeksecond_isagpst, - nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"][ - "GPSA" - ], - nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"][ - "GPSB" - ], - flat_obs.loc[mask_idx, "elevation_rad"], - flat_obs.loc[mask_idx, "azimuth_rad"], - latitude_user_rad, - longitude_user_rad, - ) - * ( - constants.carrier_frequencies_hz()["G"]["L1"][1] ** 2 - / flat_obs.loc[mask_idx, "carrier_frequency_hz"] ** 2 - ) - ) - else: + if "IONOSPHERIC CORR" not in nav_header_dict[f"{year:03d}" + f"{doy:03d}"]: logging.warning(f"Missing iono model parameters for day {doy:03d}") - iono_all_days.append(np.full(mask_idx.shape, np.nan)) - delays = np.ones((len(flat_obs.index))) * np.nan - delays[np.concatenate(idx_all_days)] = np.concatenate(iono_all_days) - return delays + continue + # Assign iono model parameters to rows + # Assignment criteria: time of emission belonging to the day of the current NAV file + gpsa = nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"][ + "GPSA" + ] + gpsb = nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"][ + "GPSB" + ] + matching_rows = pl.when(pl.col("time_of_emission_isagpst") + >= pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy - 1), + pl.col("time_of_emission_isagpst") + < pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy), pl.col("observation_type").str.starts_with("C") + ) + flat_obs = flat_obs.with_columns(matching_rows.then(gpsa).otherwise(pl.col("gpsa")).alias("gpsa"), + matching_rows.then(gpsb).otherwise(pl.col("gpsb")).alias("gpsb") + ) + time_of_emission_weeksecond_isagpst = timedelta_2_weeks_and_seconds( + (flat_obs.select("time_of_emission_isagpst") - constants.system_time_scale_rinex_utc_epoch["GPST"]).to_series().cast(dtype=pl.Duration(time_unit="ns")) + )[1] + delay = compute_l1_iono_delay_klobuchar( + time_of_emission_weeksecond_isagpst, + nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"][ + "GPSA" + ], + nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"][ + "GPSB" + ], + flat_obs.loc[mask_idx, "elevation_rad"], + flat_obs.loc[mask_idx, "azimuth_rad"], + latitude_user_rad, + longitude_user_rad, + ) * ( + constants.carrier_frequencies_hz()["G"]["L1"][1] ** 2 + / flat_obs.loc[mask_idx, "carrier_frequency_hz"] ** 2 + ) + return delay def compute_tropo_delay_saastamoinen(height, el, lat, humi=0.7): diff --git a/src/prx/main.py b/src/prx/main.py index 6923f67..5baa70b 100644 --- a/src/prx/main.py +++ b/src/prx/main.py @@ -216,7 +216,10 @@ def assign_carrier_frequencies(flat_obs): return flat_obs +from line_profiler import profile + @util.timeit +@profile def build_records_levels_12( rinex_3_obs_file, rinex_3_ephemerides_files, @@ -423,8 +426,8 @@ def build_records_levels_12( if prx_level == 2: # add iono correction - iono_delay = atmo.add_iono_column( - flat_obs.to_pandas(), + iono_delay = atmo.compute_iono_column( + flat_obs, rinex_3_ephemerides_files, approximate_receiver_ecef_position_m, ) @@ -638,7 +641,7 @@ def assign_carrier_frequencies(flat_obs): rnx3_nav_files = nav_file_discovery.discover_or_download_auxiliary_files( rinex_3_obs_file )["broadcast_ephemerides"] - iono_delay = atmo.add_iono_column( + iono_delay = atmo.compute_iono_column( flat_obs, rnx3_nav_files, approximate_receiver_ecef_position_m ) flat_obs["iono_delay_m"] = iono_delay diff --git a/src/prx/rinex_nav/evaluate.py b/src/prx/rinex_nav/evaluate.py index 68c9b04..358ca7c 100644 --- a/src/prx/rinex_nav/evaluate.py +++ b/src/prx/rinex_nav/evaluate.py @@ -668,12 +668,15 @@ def compute_parallel( # split dataframe into `n_chunks` smaller dataframes n_chunks = min(len(per_signal_query.index), 4) chunks = np.array_split(per_signal_query, n_chunks) - processed_chunks = parallel( - delayed(compute)( - rinex_nav_file_path, chunk, is_query_corrected_by_sat_clock_offset + if joblib_backend == "sequential": + processed_chunks = [compute(rinex_nav_file_path, chunk, is_query_corrected_by_sat_clock_offset) for chunk in chunks] + else: + processed_chunks = parallel( + delayed(compute)( + rinex_nav_file_path, chunk, is_query_corrected_by_sat_clock_offset + ) + for chunk in chunks ) - for chunk in chunks - ) result = pd.concat(processed_chunks) result["frequency_slot"] = result["frequency_slot"].astype(float) return result @@ -709,7 +712,6 @@ def compute( ) ) # Apply sat clock correction to the query time for satellite position computation - per_signal_query = per_signal_query.with_columns( query_time_wrt_ephemeris_reference_time_s=pl.col( "query_time_wrt_ephemeris_reference_time_s" diff --git a/src/prx/test/test_helpers.py b/src/prx/test/test_helpers.py index 63d37bf..9e1a1e4 100644 --- a/src/prx/test/test_helpers.py +++ b/src/prx/test/test_helpers.py @@ -230,7 +230,7 @@ def test_timedelta_2_weeks_and_seconds(): seconds_of_week_expected = [280800, 281400, 317400, 302400, np.nan] np.testing.assert_array_equal(week_computed, week_expected) - np.testing.assert_array_equal(seconds_of_week_computed, seconds_of_week_expected) + np.testing.assert_allclose(seconds_of_week_computed, seconds_of_week_expected, atol=1e-15, rtol=0) # We also expect the function to work for Series of timestamps week_series_computed, seconds_of_week_series_computed = ( diff --git a/src/prx/test/test_main.py b/src/prx/test/test_main.py index 55e5d30..1ec4578 100644 --- a/src/prx/test/test_main.py +++ b/src/prx/test/test_main.py @@ -178,7 +178,7 @@ def test_prx_command_line_call(input_for_test_tlse): def test_prx_function_call(input_for_test_tlse): test_file = input_for_test_tlse - main.process(observation_file_path=test_file, prx_level=2) + main.process(observation_file_path=test_file, prx_level=2, joblib_backend="sequential") expected_prx_file = Path(str(test_file).replace("crx.gz", "csv")) assert expected_prx_file.exists() df = pd.read_csv(expected_prx_file, comment="#") diff --git a/src/prx/util.py b/src/prx/util.py index 06f3b5e..5433d70 100644 --- a/src/prx/util.py +++ b/src/prx/util.py @@ -219,14 +219,19 @@ def timestamp_to_mid_day(ts): ) -def timedelta_2_weeks_and_seconds(time_delta: pd.Timedelta | pd.Series): - if isinstance(time_delta, pd.Timedelta): - wn_series, tow_series = timedelta_2_weeks_and_seconds(pd.Series([time_delta])) - return wn_series.iloc[0], tow_series.iloc[0] - in_nanoseconds = time_delta / pd.Timedelta(1, "ns") - weeks = np.floor(in_nanoseconds / constants.cNanoSecondsPerWeek) - week_nanoseconds = in_nanoseconds - weeks * constants.cNanoSecondsPerWeek - return weeks, week_nanoseconds.astype(np.float64) / constants.cNanoSecondsPerSecond +def timedelta_2_weeks_and_seconds(time_delta: pd.Timedelta | pd.Series | pl.Series): + if time_delta is pd.NaT: + return np.nan, np.nan + if isinstance(time_delta, pd.Timedelta) : + w, s = timedelta_2_weeks_and_seconds(pl.Series([time_delta.value], dtype=pl.Duration(time_unit="ns"))) + return w[0], s[0] + if isinstance(time_delta, pd.Series): + w, s = timedelta_2_weeks_and_seconds(pl.from_pandas(time_delta)) + return w, s + seconds = timedelta_2_seconds(time_delta) + weeks = (seconds / constants.cSecondsPerWeek).floor() + week_seconds = seconds - weeks * constants.cSecondsPerWeek + return weeks, week_seconds def week_and_seconds_2_timedelta(weeks, seconds): @@ -244,8 +249,11 @@ def timedelta_2_seconds( return timedelta_2_seconds(pl.from_pandas(time_delta)).to_pandas() assert isinstance(time_delta, pl.Series) assert time_delta.dtype.time_unit == "ns" + # Don't lose float resolution unnecessarily here + nanoseconds = time_delta.dt.total_nanoseconds() + integer_seconds = (nanoseconds / constants.cNanoSecondsPerSecond).floor().cast(int) return ( - time_delta.dt.total_nanoseconds().cast(float) / constants.cNanoSecondsPerSecond + integer_seconds + (nanoseconds - integer_seconds*constants.cNanoSecondsPerSecond).cast(float) / constants.cNanoSecondsPerSecond ) From f6e1d2dbc91cbe2fe6d03a76ab312601e7c6df46 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 25 Jun 2026 23:34:17 +0200 Subject: [PATCH 22/26] works --- src/prx/atmospheric_corrections.py | 65 +++++++++++++----------------- src/prx/main.py | 7 ++-- 2 files changed, 32 insertions(+), 40 deletions(-) diff --git a/src/prx/atmospheric_corrections.py b/src/prx/atmospheric_corrections.py index de0f53e..60b0ed4 100644 --- a/src/prx/atmospheric_corrections.py +++ b/src/prx/atmospheric_corrections.py @@ -100,9 +100,7 @@ def compute_iono_column( [latitude_user_rad, longitude_user_rad, __] = ecef_2_geodetic( approximate_receiver_ecef_position_m ) - flat_obs = flat_obs.with_columns( - gpsa=pl.lit(None), gpsb=pl.lit(None), iono_delay=pl.lit(None) - ) + flat_obs = flat_obs.with_row_index() for file in rinex_3_ephemerides_files: # get year and doy from NAV filename year = int(file.name[12:16]) @@ -112,40 +110,35 @@ def compute_iono_column( continue # Assign iono model parameters to rows # Assignment criteria: time of emission belonging to the day of the current NAV file - gpsa = nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"]["GPSA"] - gpsb = nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"]["GPSB"] - matching_rows = pl.when( - pl.col("time_of_emission_isagpst") - >= pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy - 1), - pl.col("time_of_emission_isagpst") - < pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy), - pl.col("observation_type").str.starts_with("C"), - ) - flat_obs = flat_obs.with_columns( - matching_rows.then(gpsa).otherwise(pl.col("gpsa")).alias("gpsa"), - matching_rows.then(gpsb).otherwise(pl.col("gpsb")).alias("gpsb"), - ) - time_of_emission_weeksecond_isagpst = timedelta_2_weeks_and_seconds( - ( - flat_obs.select("time_of_emission_isagpst") - - constants.system_time_scale_rinex_utc_epoch["GPST"] + matching_rows = ((pl.col("time_of_emission_isagpst") >= pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy - 1)) & + (pl.col("time_of_emission_isagpst") < pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy)) & + pl.col("observation_type").str.starts_with("C")) + day_df = flat_obs.filter(matching_rows) + time_of_emission_weeksecond_isagpst = timedelta_2_weeks_and_seconds( + ( + day_df.select("time_of_emission_isagpst") + - constants.system_time_scale_rinex_utc_epoch["GPST"] + ) + .to_series() + .cast(dtype=pl.Duration(time_unit="ns")) + )[1] + delay = compute_l1_iono_delay_klobuchar( + time_of_emission_weeksecond_isagpst.to_numpy().reshape((-1, 1)), + nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"]["GPSA"], + nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"]["GPSB"], + day_df.select(pl.col("elevation_rad")).to_numpy().reshape((-1, 1)), + day_df.select(pl.col("azimuth_rad")).to_numpy().reshape((-1, 1)), + latitude_user_rad, + longitude_user_rad, + ) * ( + constants.carrier_frequencies_hz()["G"]["L1"][1] ** 2 + / day_df.select(pl.col("carrier_frequency_hz")).to_numpy().reshape((-1, 1)) ** 2 + ) + day_df = day_df.select(["index"]).with_columns( + pl.Series("iono_delay_m", delay.flatten()) ) - .to_series() - .cast(dtype=pl.Duration(time_unit="ns")) - )[1] - delay = compute_l1_iono_delay_klobuchar( - time_of_emission_weeksecond_isagpst, - nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"]["GPSA"], - nav_header_dict[f"{year:03d}" + f"{doy:03d}"]["IONOSPHERIC CORR"]["GPSB"], - flat_obs.loc[mask_idx, "elevation_rad"], - flat_obs.loc[mask_idx, "azimuth_rad"], - latitude_user_rad, - longitude_user_rad, - ) * ( - constants.carrier_frequencies_hz()["G"]["L1"][1] ** 2 - / flat_obs.loc[mask_idx, "carrier_frequency_hz"] ** 2 - ) - return delay + flat_obs = flat_obs.join(day_df, on="index", how="left") + return flat_obs.drop("index") def compute_tropo_delay_saastamoinen(height, el, lat, humi=0.7): diff --git a/src/prx/main.py b/src/prx/main.py index d736179..f5ec833 100644 --- a/src/prx/main.py +++ b/src/prx/main.py @@ -253,7 +253,7 @@ def build_records_levels_12( log.info("Computing times of emission in satellite time") per_sat = flat_obs.pivot( index=["time_of_reception_in_receiver_time", "satellite"], - columns=["observation_type"], + on=["observation_type"], values="observation_value", ) per_sat = per_sat.with_columns( @@ -427,12 +427,11 @@ def build_records_levels_12( if prx_level == 2: # add iono correction - iono_delay = atmo.compute_iono_column( + flat_obs = atmo.compute_iono_column( flat_obs, rinex_3_ephemerides_files, approximate_receiver_ecef_position_m, ) - flat_obs = flat_obs.with_columns(iono_delay_m=iono_delay) return flat_obs.to_pandas() @@ -470,7 +469,7 @@ def build_records_level_3( log.info("Computing times of emission in satellite time") per_sat = flat_obs.pivot( index=["time_of_reception_in_receiver_time", "satellite"], - columns=["observation_type"], + on=["observation_type"], values="observation_value", ).reset_index() per_sat["time_scale"] = ( From 9394be78a79a6b881129cfd5460a57aee4a75d6d Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 25 Jun 2026 23:37:35 +0200 Subject: [PATCH 23/26] drop --- src/prx/atmospheric_corrections.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/prx/atmospheric_corrections.py b/src/prx/atmospheric_corrections.py index 60b0ed4..ea59303 100644 --- a/src/prx/atmospheric_corrections.py +++ b/src/prx/atmospheric_corrections.py @@ -113,7 +113,7 @@ def compute_iono_column( matching_rows = ((pl.col("time_of_emission_isagpst") >= pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy - 1)) & (pl.col("time_of_emission_isagpst") < pd.Timestamp(year=year, month=1, day=1) + pd.Timedelta(days=doy)) & pl.col("observation_type").str.starts_with("C")) - day_df = flat_obs.filter(matching_rows) + day_df = flat_obs.filter(matching_rows).select("index", "time_of_emission_isagpst", "elevation_rad", "azimuth_rad", "carrier_frequency_hz") time_of_emission_weeksecond_isagpst = timedelta_2_weeks_and_seconds( ( day_df.select("time_of_emission_isagpst") From 13d10f43f2701f33d6cd4a3f9c30e73d89d81e14 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 25 Jun 2026 23:42:51 +0200 Subject: [PATCH 24/26] about 10% faster From b86d228d74c37008ffe9555173e9b7a802579ab9 Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Thu, 25 Jun 2026 23:51:08 +0200 Subject: [PATCH 25/26] drop unnessesary conversion --- src/prx/main.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/prx/main.py b/src/prx/main.py index f5ec833..af536f2 100644 --- a/src/prx/main.py +++ b/src/prx/main.py @@ -27,11 +27,10 @@ @util.timeit def write_prx_file( prx_header: dict, - prx_records_pd: pd.DataFrame, + prx_records: pl.DataFrame, file_name_without_extension: Path, ): output_file = Path(f"{str(file_name_without_extension)}.csv") - prx_records = pl.from_pandas(prx_records_pd) prx_records = prx_records.with_columns( (pl.col("elevation_rad") * cDegPerRad).alias("sat_elevation_deg"), (pl.col("azimuth_rad") * cDegPerRad).alias("sat_azimuth_deg"), @@ -433,7 +432,7 @@ def build_records_levels_12( approximate_receiver_ecef_position_m, ) - return flat_obs.to_pandas() + return flat_obs def build_records_level_3( @@ -702,13 +701,13 @@ def process( metadata["processing_start_time"] = t0 # build record - records = build_records_level_3( + records = pl.from_pandas(build_records_level_3( rinex_3_obs_file, aux_files["sp3_orb"], aux_files["atx"], metadata["approximate_receiver_ecef_position_m"], model_tropo, - ) + )) metadata["processing_time"] = str( pd.Timestamp.now() - metadata["processing_start_time"] ) From 53d8bef1c39f2cd3b95b987756f76d878d5a514b Mon Sep 17 00:00:00 2001 From: Jan Bolting Date: Fri, 26 Jun 2026 00:41:31 +0200 Subject: [PATCH 26/26] bof --- src/prx/main.py | 14 -------------- src/prx/test/benchmark.py | 9 +++++---- 2 files changed, 5 insertions(+), 18 deletions(-) diff --git a/src/prx/main.py b/src/prx/main.py index af536f2..6f8b2f3 100644 --- a/src/prx/main.py +++ b/src/prx/main.py @@ -619,20 +619,6 @@ def build_records_level_3( # set frequency slot to 1 for non-GLONASS satellites flat_obs.loc[flat_obs.satellite.str[0] != "R", "frequency_slot"] = int(1) - def assign_carrier_frequencies(flat_obs): - freq_dict = pd.json_normalize(carrier_frequencies_hz(), sep="_").to_dict( - orient="records" - )[0] - assignable = flat_obs.frequency_slot.notna() - keys = ( - flat_obs.satellite[assignable].str[0] - + "_L" - + flat_obs["observation_type"][assignable].str[1] - + "_" - + flat_obs.frequency_slot[assignable].astype(int).astype(str) - ) - flat_obs.loc[:, "carrier_frequency_hz"] = keys.map(freq_dict) - return flat_obs flat_obs = assign_carrier_frequencies(flat_obs).drop(columns=["frequency_slot"]) diff --git a/src/prx/test/benchmark.py b/src/prx/test/benchmark.py index 627e3c6..e8997e7 100644 --- a/src/prx/test/benchmark.py +++ b/src/prx/test/benchmark.py @@ -35,7 +35,9 @@ def run_case(case: dict, ram: bool, warm_parser_cache: bool) -> pd.DataFrame: p = cProfile.Profile() p.enable() # cProfile does not profile subprocesses with joblib's loky backend, use threading + t0 = pd.Timestamp.now() process(observation_file_path=obs_file, joblib_backend="threading") + t_exec = pd.Timestamp.now() - t0 p.disable() # RAM @@ -55,7 +57,9 @@ def run_case(case: dict, ram: bool, warm_parser_cache: bool) -> pd.DataFrame: reader = memray.FileReader(memray_output) metadata = reader.metadata peak_ram_mb = metadata.peak_memory / 1024 / 1024 - + logger.info( + f"Processed {obs_file.name} in {t_exec}: {case['epochs'] / t_exec.seconds} epochs/s, peak RAM [Mb]: {peak_ram_mb}" + ) # Run time stats_file = Path("benchmark_prx.prof").resolve() p.dump_stats(stats_file) @@ -67,9 +71,6 @@ def run_case(case: dict, ram: bool, warm_parser_cache: bool) -> pd.DataFrame: .sort_values(by="tottime", ascending=False) .reset_index(drop=True) ) - logger.info( - f"Processed {obs_file.name} in {df.iloc[0, :]['tottime']} seconds: {case['epochs'] / df.iloc[0, :]['tottime']} epochs/s, peak RAM [Mb]: {peak_ram_mb}" - ) df = df[["func", "tottime"]] df["function"] = df["func"].apply(lambda x: getattr(x, "co_name", None)) df["file"] = df["func"].apply(lambda x: getattr(x, "co_filename", None))