Skip to content

add saastamoinen tropo delay computation #185

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jtec merged 10 commits into from
Dec 2, 2025
Merged

add saastamoinen tropo delay computation #185

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
6f4ed33
- add saastamoinen tropo delay computation
plutonheaven Jul 7, 2025
f1905cd
Merge branch 'main' into 202507_tropo_saastamoinen
plutonheaven Jul 9, 2025
5920368
Merge branch 'main' into 202507_tropo_saastamoinen
plutonheaven Jul 23, 2025
1cb614a
add choice of tropo model to command line arg
plutonheaven Jul 23, 2025
06d946d
add test with call to specific tropo model
plutonheaven Jul 23, 2025
4508f9d
add documentation about this feature
plutonheaven Jul 23, 2025
c623fc2
Merge remote-tracking branch 'origin/main' into 202507_tropo_saastamo…
jtec Dec 2, 2025
e15cac2
Merge remote-tracking branch 'origin/main' into 202507_tropo_saastamo…
jtec Dec 2, 2025
515e621
fix ephemeris discovery
jtec Dec 2, 2025
12a7d86
Merge branch 'main' into 202507_tropo_saastamoinen
jtec Dec 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,8 @@ uv run python src/prx/main.py --observation_file_path <path_to_rinex_file>
You can specify `--prx_level` according to the type of computation you need. `--prx_level 1` is adapted for DGNSS or RTK
processing, `--prx_level 2` is adapted for SPP processing and is the default value.

There is an optional argument to select the tropospheric delay model, by adding `--tropo saastamoinen` (default) or `--tropo unb3m`

You might have to add `<path to prx root>/src/prx` to your `PYTHONPATH` environment variable if you run
into import errors.

Expand Down
68 changes: 36 additions & 32 deletions documents/dev_status.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,11 @@
## Features of PRX

Different information are required depending oh the positioning algorithm chosen. Three types of positioning algorithms have been identified, each requiring specific data that are organized into different levels of processing in PRX.
Different information are required depending oh the positioning algorithm chosen. Three types of positioning algorithms
have been identified, each requiring specific data that are organized into different levels of processing in PRX.
The following tables gather the different information needed according to the level of processing.
- **Level 1** : suitable for DGNSS/RTK processing.
- **Level 2** : designed for Single Point Positioning (SPP).
- **Level 3** : designed for Precise Point Positioning (PPP).
- **Level 1** : suitable for DGNSS/RTK processing.
- **Level 2** : designed for Single Point Positioning (SPP).
- **Level 3** : designed for Precise Point Positioning (PPP).

### Level 1 – DGNSS / RTK

Expand All @@ -14,50 +15,53 @@ Those parameters are computed from the broadcast navigation message (`rinex nav`
uv run python src/prx/main.py --observation_file_path <path_to_rinex_file> --prx_level 1
```

| Parameters | Name in PRX file | Status |
|----------------------------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|-----------|
| GNSS observations | `time_of_reception_in_receiver_time`<br/>`C_obs_m`, `D_obs_hz`, `L_obs_cycles`, `S_obs_dBHz`,<br/>`rnx_obs_identifier`, `constellation`, `prn`, `frequency_slot` (for GLONASS) | ✅ |
| Loss of Lock indicator | `LLI` | ✅ |
| Satellite health flag | `health_flag` | ✅ |
| Satellite position and velocity | `sat_pos_x_m`, `sat_pos_y_m`, `sat_pos_z_m`,<br> `sat_vel_x_mps`, `sat_vel_y_mps`, `sat_vel_z_mps` | ✅ |
| Satellite clock offset and drift | `sat_clock_offset_m`, `sat_clock_drift_mps`, `relativistic_clock_effect_m` | ✅ |
| Satellite elevation and azimuth | `sat_elevation_deg`, `sat_azimuth_deg` | ✅ |
| Ephemerides dataset identifier | `ephemeris_hash` | ✅ |
| Parameters | Name in PRX file | Status |
|----------------------------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|--------|
| GNSS observations | `time_of_reception_in_receiver_time`<br/>`C_obs_m`, `D_obs_hz`, `L_obs_cycles`, `S_obs_dBHz`,<br/>`rnx_obs_identifier`, `constellation`, `prn`, `frequency_slot` (for GLONASS) | ✅ |
| Loss of Lock indicator | `LLI` | ✅ |
| Satellite health flag | `health_flag` | ✅ |
| Satellite position and velocity | `sat_pos_x_m`, `sat_pos_y_m`, `sat_pos_z_m`,<br> `sat_vel_x_mps`, `sat_vel_y_mps`, `sat_vel_z_mps` | ✅ |
| Satellite clock offset and drift | `sat_clock_offset_m`, `sat_clock_drift_mps`, `relativistic_clock_effect_m` | ✅ |
| Satellite elevation and azimuth | `sat_elevation_deg`, `sat_azimuth_deg` | ✅ |
| Ephemerides dataset identifier | `ephemeris_hash` | ✅ |

### Level 2 - SPP

### Level 2 - SPP
Same as Level 1, with the following additional parameters. Those parameters are still computed from the broadcast navigation message.
Same as Level 1, with the following additional parameters. Those parameters are still computed from the broadcast
navigation message.

```
uv run python src/prx/main.py --observation_file_path <path_to_rinex_file> --prx_level 2
uv run python src/prx/main.py --observation_file_path <path_to_rinex_file> --prx_level 2 \
[--tropo {saastamoinen | unb3m}]
```

| Parameters | Name in PRX file | Status |
|----------------------------------|--------------------------------|-----------|
| Sagnac effect | `sagnac_effect_m` | ✅ |
| Tropospheric delay | `tropo_delay_m` | ✅ |
| Satellite code bias | `sat_code_bias_m` | ✅ |
| Ionospheric delay | `iono_delay_m` | ✅ |
| Parameters | Name in PRX file | Status |
|--------------------------------------------------|-------------------|--------|
| Sagnac effect | `sagnac_effect_m` | ✅ |
| Tropospheric delay (Saastamoinen or UNB3m model) | `tropo_delay_m` | ✅ |
| Satellite code bias | `sat_code_bias_m` | ✅ |
| Ionospheric delay | `iono_delay_m` | ✅ |

---

### Level 3 - PPP
### Level 3 - PPP

Same as level 1 and 2, with the following additional parameters. The parameters are computed using IGS products.

```
uv run python src/prx/main.py --observation_file_path <path_to_rinex_file> --prx_level 3
```

| Parameters | Name in PRX file | Status |
|-------------------------------------------------------------------|---------------------------------------------------------------------------------------------------|--------|
| Satellite position and velocity (computed using `sp3` and `clk` files) | `sat_pos_x_m`, `sat_pos_y_m`, `sat_pos_z_m`,<br>`sat_vel_x_mps`, `sat_vel_y_mps`, `sat_vel_z_mps` | ❌ |
| Parameters | Name in PRX file | Status |
|---------------------------------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------|--------|
| Satellite position and velocity (computed using `sp3` and `clk` files) | `sat_pos_x_m`, `sat_pos_y_m`, `sat_pos_z_m`,<br>`sat_vel_x_mps`, `sat_vel_y_mps`, `sat_vel_z_mps` | ❌ |
| Satellite clock offset and drift (including relativistic effect) (computed using `sp3` and `clk` files) | `sat_clock_offset_m`, `sat_clock_drift_mps` | ❌ |
| Tropospheric delay (computed using `tropex` files) | `tropo_delay_m` | ❌ |
| Tropospheric mapping function | to be completed | ❌ |
| Ionospheric delay (computed using `ionex` files) | `iono_delay_m` | ❌ |
| Satellite code & phase bias (computed using `bia` files) | `sat_code_bias_m` | ❌ |
| Satellite & receiver phase center offset and variation (computed using `antex` files) | to be completed | ❌ |
| Solid Earth Tide | to be completed | ❌ |
| Tropospheric delay (computed using `tropex` files) | `tropo_delay_m` | ❌ |
| Tropospheric mapping function | to be completed | ❌ |
| Ionospheric delay (computed using `ionex` files) | `iono_delay_m` | ❌ |
| Satellite code & phase bias (computed using `bia` files) | `sat_code_bias_m` | ❌ |
| Satellite & receiver phase center offset and variation (computed using `antex` files) | to be completed | ❌ |
| Solid Earth Tide | to be completed | ❌ |



Expand Down
120 changes: 89 additions & 31 deletions src/prx/atmospheric_corrections.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,52 @@ def add_iono_column(
return np.concatenate(idx_all_days), np.concatenate(iono_all_days)


def compute_tropo_delay_saastamoinen(height, el, lat, humi=0.7):
"""
code source: rtklib v2.4.3 b33, rtkcmn.c, tropmodel

INPUTS:
- height: geodetic height in m. np.ndarray of shape (n,)
- el in rad. np.ndarray of shape (n,)
- lat in rad. np.ndarray of shape (n,)
- humi in ratio. scalar

OUTPUTS:
- tropo_tot, tropo_h, tropo_w in m. np.ndarray of shape (n,)
"""
assert len(height) == len(el) == len(lat)
assert (-100 <= height).all()
assert (height <= 1e4).all()
assert (
el[~np.isnan(el)] >= -np.deg2rad(5)
).all() # consider slightly negative elevation...

# standard atmosphere model
height_new = height.copy() # to avoid changing input
height_new[height_new < 0] = 0
height = height_new

p = 1013.25 * (1 - 2.2557e-5 * height) ** 5.2568
t = 15 - 6.5e-3 * height + 273.16
e_wat = 6.108 * np.exp((17.15 * t - 4684) / (t - 38.45)) * humi

# zenith angle in rad
z = np.pi / 2 - el

# hydrostatic tropospheric delay
tropo_h = (
0.0022768 * p / (1 - 0.00266 * np.cos(2 * lat) - 0.28e-6 * height) / np.cos(z)
)

# wet tropospheric delay
tropo_w = 0.002277 / np.cos(z) * (1255 / t + 0.05) * e_wat

# total tropospheric delay (projected on the slant path)
tropo_tot = tropo_h + tropo_w

return tropo_tot, tropo_h, tropo_w


def compute_tropo_delay_unb3m(
latitude_user_rad, height_user_m, day_of_year, elevation_sat_rad
):
Expand Down Expand Up @@ -372,36 +418,48 @@ def compute_tropo_delay_unb3m(
)


def add_tropo_column(sat_states, flat_obs, receiver_ecef_position_m):
# We need Timestamps to compute tropo delays
sat_states = sat_states.merge(
flat_obs[
[
"satellite",
"time_of_emission_isagpst",
"time_of_reception_in_receiver_time",
]
].drop_duplicates(),
on=["satellite", "time_of_emission_isagpst"],
how="left",
)
def compute_tropo_delay(sat_states, flat_obs, receiver_ecef_position_m, model):
"""
model is either "saastamoinen" or "unb3m"
"""
[latitude_user_rad, __, height_user_m] = ecef_2_geodetic(receiver_ecef_position_m)
days_of_year = np.array(
sat_states["time_of_reception_in_receiver_time"]
.apply(lambda element: element.timetuple().tm_yday)
.to_numpy()
)
(
tropo_delay_m,
__,
__,
__,
__,
) = compute_tropo_delay_unb3m(
latitude_user_rad * np.ones(days_of_year.shape),
height_user_m * np.ones(days_of_year.shape),
days_of_year,
sat_states.elevation_rad.to_numpy(),
)
# sat_states = sat_states.drop(columns=["time_of_reception_in_receiver_time"])
match model:
case "saastamoinen":
tropo_delay_m, _, _ = compute_tropo_delay_saastamoinen(
height_user_m * np.ones((len(sat_states),)),
sat_states.elevation_rad.to_numpy(),
latitude_user_rad * np.ones((len(sat_states),)),
)
case "unb3m":
# We need Timestamps to compute tropo delays
sat_states = sat_states.merge(
flat_obs[
[
"satellite",
"time_of_emission_isagpst",
"time_of_reception_in_receiver_time",
]
].drop_duplicates(),
on=["satellite", "time_of_emission_isagpst"],
how="left",
)
days_of_year = np.array(
sat_states["time_of_reception_in_receiver_time"]
.apply(lambda element: element.timetuple().tm_yday)
.to_numpy()
)
(
tropo_delay_m,
__,
__,
__,
__,
) = compute_tropo_delay_unb3m(
latitude_user_rad * np.ones((len(sat_states),)),
height_user_m * np.ones((len(sat_states),)),
days_of_year,
sat_states.elevation_rad.to_numpy(),
)
case _:
assert False, f"tropospheric model not recognized: {model}"
return tropo_delay_m
16 changes: 12 additions & 4 deletions src/prx/main.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -183,6 +183,7 @@ def build_records_levels_12(
rinex_3_ephemerides_files,
approximate_receiver_ecef_position_m,
prx_level,
model_tropo,
):
"""
Creates a flat_obs dataframe including columns for prx processing levels 1 and 2.
Expand Down Expand Up @@ -322,8 +323,8 @@ def build_records_levels_12(
sat_states[["sat_pos_x_m", "sat_pos_y_m", "sat_pos_z_m"]].to_numpy(),
approximate_receiver_ecef_position_m,
)
sat_states["tropo_delay_m"] = atmo.add_tropo_column(
sat_states, flat_obs, approximate_receiver_ecef_position_m
sat_states["tropo_delay_m"] = atmo.compute_tropo_delay(
sat_states, flat_obs, approximate_receiver_ecef_position_m, model_tropo
)

# Merge sat states into observation dataframe. Due to Galileo's FNAV/INAV ephemerides
Expand Down Expand Up @@ -378,7 +379,7 @@ def assign_carrier_frequencies(flat_obs):


@prx.util.timeit
def process(observation_file_path: Path, prx_level=2):
def process(observation_file_path: Path, prx_level=2, model_tropo="saastamoinen"):
t0 = pd.Timestamp.now()
# We expect a Path, but might get a string here:
observation_file_path = Path(observation_file_path)
Expand Down Expand Up @@ -406,6 +407,7 @@ def process(observation_file_path: Path, prx_level=2):
aux_files["broadcast_ephemerides"],
metadata["approximate_receiver_ecef_position_m"],
prx_level,
model_tropo,
)
case 3:
assert False, (
Expand Down Expand Up @@ -439,11 +441,17 @@ def process(observation_file_path: Path, prx_level=2):
choices=[1, 2, 3],
default=2,
)
parser.add_argument(
"--tropo",
type=str,
choices=["saastamoinen", "unb3m"],
default="saastamoinen",
)
args = parser.parse_args()
if args.observation_file_path is None:
log.error("No observation file path provided.")
sys.exit(1)
if not Path(args.observation_file_path).exists():
log.error(f"Observation file {args.observation_file_path} does not exist.")
sys.exit(1)
process(Path(args.observation_file_path), args.prx_level)
process(Path(args.observation_file_path), args.prx_level, args.tropo)
13 changes: 9 additions & 4 deletions src/prx/rinex_nav/nav_file_discovery.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,12 +157,17 @@ def discover_or_download_ephemerides(
t_start: pd.Timestamp, t_end: pd.Timestamp, folder, constellations
):
# If there are any navigation files provided by the user, use them, otherwise use IGS files.
candidates = [anything_to_rinex_3(file) for file in folder.rglob("*")]
candidates = [candidate for candidate in candidates if candidate is not None]
user_provided_nav_files = [
anything_to_rinex_3(f)
for f in folder.rglob("*")
if is_rinex_3_nav_file(anything_to_rinex_3(f))
and is_rinex_3_mixed_mgex_broadcast_ephemerides_file(anything_to_rinex_3(f))
f
for f in candidates
if (
is_rinex_3_nav_file(f)
and is_rinex_3_mixed_mgex_broadcast_ephemerides_file(f)
)
]
user_provided_nav_files = list(set(user_provided_nav_files))
if len(user_provided_nav_files) > 0:
return user_provided_nav_files
# Ephemeris files cover at least a day, so first round time stamps to midday here
Expand Down
Loading