geodezyx.gnss_edu package

Created on 12/08/2025 14:58:49

@author: psakic

Submodules

geodezyx.gnss_edu.gnss_edu module

Created on Tue Feb 11 16:39:34 2025

@author: snahmani

geodezyx.gnss_edu.gnss_edu.Sagnac_rotate_around_z(row)

Calcule la correction Sagnac pour les coordonnées satellites. Retourne un pd.Series avec les nouvelles colonnes X_sat_corr, Y_sat_corr, Z_sat_corr.

geodezyx.gnss_edu.gnss_edu.best_prn_for_interval(df: DataFrame, t0: Timestamp, t1: Timestamp, snr_col: str, snr_min: float, sampling: Timedelta, pool_prns: list[str]) → str | None

Pick a PRN from pool_prns that fully covers [t0, t1] with SNR>=snr_min. Tie-breaker: highest mean SNR over the interval.

Returns None if no PRN can fully cover the interval.

geodezyx.gnss_edu.gnss_edu.build_active_pivot_schedule(df: DataFrame, selected_prns: list[str], snr_col: str = 'S1', snr_min: float = 40.0, sampling: Timedelta = Timedelta('0 days 00:00:30'), switch_margin_db: float = 2.0) → Series

Build a stable active pivot schedule from a pre-selected PRN set. Switch only if the best candidate is better than current by switch_margin_db.

Returns:

active_pivot

Return type:

pd.Series indexed by expected epochs, values are PRN or None

geodezyx.gnss_edu.gnss_edu.build_residual_dataframe(df_obs_used: DataFrame, A: ndarray, B: ndarray, dP_est: ndarray, predicted_col: str = 'predicted_m', residual_col: str = 'residual_m') → DataFrame

Build a residual DataFrame aligned with the observation DataFrame actually used in the least-squares adjustment.

Parameters:

• df_obs_used (pd.DataFrame) – Observation DataFrame used to build A and B. Its row order must match the row order of A and B.

• A (np.ndarray) – Design matrix.

• B (np.ndarray) – Observation vector.

• dP_est (np.ndarray) – Estimated parameter vector.

• predicted_col (str, default="predicted_m") – Name of the column storing the predicted values A @ dP_est.

• residual_col (str, default="residual_m") – Name of the column storing the residuals B - A @ dP_est.

Returns:

Copy of df_obs_used enriched with predicted values and residuals.

Return type:

pd.DataFrame

geodezyx.gnss_edu.gnss_edu.check_full_coverage_from_active_pivot(active_pivot: Series) → dict

Simple proof helper: - counts None epochs - returns first/last None epoch if any

geodezyx.gnss_edu.gnss_edu.detect_intra_arc_holes(df, sampling=Timedelta('0 days 00:00:30'), gap=Timedelta('0 days 00:30:00'))

Detect missing observations inside satellite arcs.

A hole is defined as:

sampling*2 < dt < gap

Parameters:

• df (pandas.DataFrame) – GNSS dataframe indexed by (epoch, prn)

• sampling (Timedelta) – Expected RINEX sampling interval

• gap (Timedelta) – Arc segmentation threshold

Returns:

epoch_prev : last valid observation epoch : first observation after the hole dt : time difference n_missing : number of missing epochs

Return type:

DataFrame with columns

geodezyx.gnss_edu.gnss_edu.enforce_min_dwell(df: DataFrame, active_pivot: Series, snr_col: str, snr_min: float, sampling: Timedelta, min_dwell: Timedelta) → Series

Stability post-filter: If a segment is shorter than min_dwell, try to absorb it into the previous pivot if the previous pivot can cover the short interval.

This reduces last-minute tiny segments (e.g. 25 min) without harming coverage.

geodezyx.gnss_edu.gnss_edu.expected_epochs(t0: Timestamp, t1: Timestamp, sampling: Timedelta) → DatetimeIndex

Return the expected epoch grid (inclusive endpoints) for a given sampling.

geodezyx.gnss_edu.gnss_edu.fix_short_pivot_segments(df: DataFrame, active_pivot: Series, snr_col: str, snr_min: float, sampling: Timedelta, min_duration: Timedelta, pool: str = 'selected_only', selected_prns: list[str] | None = None) → Series

Remove segments shorter than min_duration.

Priority

1. merge into previous pivot if it covers the short interval

2. merge into next pivot if it covers the short interval

3. fallback: pick best PRN that covers the short interval (optional pool)

returns:

active_pivot_fixed

rtype:

pd.Series

geodezyx.gnss_edu.gnss_edu.get_approx_position(fichier)

Recherche la ligne contenant ‘APPROX POSITION XYZ’ dans le fichier et retourne les trois premières valeurs sous forme de numpy.array.

geodezyx.gnss_edu.gnss_edu.greedy_pivot_set_cover(df: DataFrame, snr_col: str, snr_min: float, sampling: Timedelta = Timedelta('0 days 00:00:30'), require_full_coverage: bool = True, return_diagnostics: bool = True)

Greedy set-cover selection of PRNs to cover ALL expected epochs with SNR >= snr_min.

Guarantee: - If require_full_coverage=True: raise RuntimeError if full coverage is impossible.

Returns:

• selected_prns (list[str])

• diagnostics (dict (optional))

geodezyx.gnss_edu.gnss_edu.grep_file(pattern, filename)

Recherche un motif dans un fichier et retourne les lignes correspondantes.

geodezyx.gnss_edu.gnss_edu.hatch_carrier_smoothing_by_satellite(df: DataFrame, code_col: str, phase_col: str, window: int = 100, sampling_seconds: float = 30.0, max_gap_factor: float = 1.5, slip_threshold_m: float | None = 10.0) → Series

Apply Hatch carrier smoothing to a code observable independently for each satellite.

This function implements a simple code-phase smoothing strategy in which the code observable provides the absolute range level, while the carrier phase provides a much less noisy short-term range evolution.

The smoothing is performed separately for each satellite track, assuming that the input DataFrame uses a MultiIndex containing at least:

• ‘epoch’ : observation epoch

• ‘prn’ : satellite identifier

Parameters:

• df (pd.DataFrame) – Input observation DataFrame indexed by a MultiIndex containing at least (‘epoch’, ‘prn’).

• code_col (str) – Name of the code observable column, in meters.

• phase_col (str) – Name of the carrier-phase observable column, in meters. Important: this must already be expressed in meters, not in cycles.

• window (int, default=100) –

  Maximum Hatch smoothing length, in number of epochs.

  The effective smoothing duration is approximately:

  window * sampling_seconds

  With the default values used in this TP:

  100 * 30 s = 3000 s ≈ 50 min

  A larger window gives stronger smoothing but requires longer continuity of the phase observable.

• sampling_seconds (float, default=30.0) –

  Nominal time interval, in seconds, between two consecutive epochs for a given satellite.

  This is used to detect abnormal time gaps. For example:

  • 30.0 for a 30-second RINEX

  • 1.0 for a 1-second RINEX

• max_gap_factor (float, default=1.5) –

  Time-gap tolerance factor used to decide whether the smoothing filter must be reset.

  The filter is reset when the time interval between two consecutive epochs exceeds:

  max_gap_factor * sampling_seconds

  With the default values:

  1.5 * 30 s = 45 s

  Therefore, the filter is reset if a time gap larger than 45 seconds is detected.

• slip_threshold_m (float | None, default=10.0) –

  Optional crude consistency threshold, in meters, used as a simple protection against cycle slips or other discontinuities.

  At each epoch, the function compares:

  • the code increment : P(k) - P(k-1)

  • the phase increment : L(k) - L(k-1)

  If:

  abs( [P(k)-P(k-1)] - [L(k)-L(k-1)] ) > slip_threshold_m

  the filter is reset.

  This is only a simple pedagogical safeguard, not a rigorous cycle-slip detection method.

  If set to None, this additional test is disabled.

Returns:

Smoothed code observable, in meters, aligned with df.index.

Return type:

pd.Series

Notes

The recursive Hatch formula is:

P_tilde(k) = (1/n) * P(k)

• ((n - 1)/n) * [P_tilde(k-1) + L(k) - L(k-1)]

where:

• P_tilde(k) is the smoothed code,

• P(k) is the raw code,

• L(k) is the carrier phase in meters,

• n is the current smoothing length, increasing progressively up to the maximum value ‘window’.

The filter is reset when:

• the code or phase is missing,

• the time gap is too large,

• the optional code-phase increment consistency test fails.

In this TP, the default values are chosen for a simple and robust introductory use case based on 30-second GNSS data.

geodezyx.gnss_edu.gnss_edu.hatch_carrier_smoothing_by_satellite_fb(df: DataFrame, code_col: str, phase_col: str, window: int = 100, sampling_seconds: float = 30.0, max_gap_factor: float = 1.5, slip_threshold_m: float | None = 10.0) → DataFrame

Apply forward-backward Hatch carrier smoothing independently for each satellite.

The function returns three aligned series:

• forward Hatch smoothing

• backward Hatch smoothing

• forward-backward average

Parameters:

• df (pd.DataFrame) – Input observation DataFrame indexed by a MultiIndex containing at least (‘epoch’, ‘prn’).

• code_col (str) – Name of the code observable column, in meters.

• phase_col (str) – Name of the carrier-phase observable column, in meters. Important: this phase must already be expressed in meters.

• window (int, default=100) – Maximum Hatch smoothing length, in number of epochs.

• sampling_seconds (float, default=30.0) – Nominal sampling interval in seconds.

• max_gap_factor (float, default=1.5) – Reset the filter when the time gap exceeds max_gap_factor * sampling_seconds.

• slip_threshold_m (float | None, default=10.0) – Optional crude consistency threshold, in meters, applied to the difference between code and phase increments. If exceeded, the filter is reset. If None, no such test is applied.

Returns:

DataFrame with three columns:

• f”{code_col}_smooth_fwd”

• f”{code_col}_smooth_bwd”

• f”{code_col}_smooth_fb”

Return type:

pd.DataFrame

Notes

The forward-backward average is a simple way to reduce the edge effects of a purely forward causal smoother. It remains pedagogical and easy to inspect, but it is not an optimal state-space smoother.

geodezyx.gnss_edu.gnss_edu.load_and_clean_rinex(path)

Charge et nettoie un fichier RINEX d’observation.

Paramètres :

• path : chemin vers le fichier RINEX.

Retourne :

• dfDataFrame nettoyé, indexé par [‘epoch’, ‘prn’],

  et contenant une colonne ‘ind_ligne’ indiquant le numéro de ligne.

geodezyx.gnss_edu.gnss_edu.pivot_schedule_to_segments(active_pivot: Series, sampling: Timedelta, drop_none: bool = True) → DataFrame

Convert an active pivot schedule to segments: columns = prn, start, end, duration, n_epochs.

geodezyx.gnss_edu.gnss_edu.plot_gnss_sd_by_prn(df_SD: DataFrame, observable: str = 'SD_L1', **kwargs)

Convenience wrapper for SD plots (keeps the TP API stable).

geodezyx.gnss_edu.gnss_edu.plot_gnss_timeseries_by_prn(df: DataFrame, y: str, gap: Timedelta = Timedelta('0 days 00:30:00'), label_arcs: bool = True, arc_label_fontsize: int = 8, arc_label_xytext: tuple[int, int] = (3, 3), show_legend: bool = False, legend_outside: bool = True, legend_ncol: int = 4, title: str | None = None, xlabel: str = 'Time (epoch)', ylabel: str | None = None, figsize: tuple[int, int] = (10, 6))

Plot a GNSS time series by satellite PRN with gap handling.

This function is generic: it works for any DataFrame that:

• is indexed by a MultiIndex with levels (‘epoch’, ‘prn’)

• contains a numeric column y

Parameters:

• df (pandas.DataFrame) – MultiIndex (epoch, prn).

• y (str) – Column name to plot.

• gap (pandas.Timedelta) – Time gap threshold used to split continuous segments (“arcs”).

• label_arcs (bool) – If True, write PRN text at the start of each arc.

• arc_label_fontsize (int) – Font size for PRN labels.

• arc_label_xytext (tuple) – (dx, dy) offset in points for the arc labels.

• show_legend (bool) – If True, show legend (often unnecessary if label_arcs=True).

• legend_outside (bool) – If True, place legend outside the axes (right side).

• legend_ncol (int) – Legend columns.

• title (str or None) – Plot title. If None, a default title is generated.

• xlabel (str) – X axis label.

• ylabel (str or None) – Y axis label. If None, uses y.

• figsize (tuple) – Figure size.

Return type:

fig, ax

geodezyx.gnss_edu.gnss_edu.plot_residual_analysis(*args, df_residuals: DataFrame | None = None, figure_title: str | None = None, save_path=None, P_est=None, P_rnx_header=None, residual_col: str = 'residual_m', predicted_col: str = 'predicted_m', top_left_plot: str = 'timeseries', prn_gap: Timedelta = Timedelta('0 days 00:30:00'), label_arcs: bool = True, **kwargs)

Public residual-analysis function with backward-compatible routing.

Preferred API

plot_residual_analysis(

df_residuals=df_residuals, figure_title=…, save_path=…, P_est=…, P_rnx_header=…, top_left_plot=”timeseries” or “by_prn”,

)

Legacy API (deprecated)

plot_residual_analysis(

A, B, dP_est, figure_title=…, save_path=…, P_est=…, P_rnx_header=…,

)

geodezyx.gnss_edu.gnss_edu.plot_sd_derivative_by_prn(df_SD: DataFrame, obs: str, gap: Timedelta = Timedelta('0 days 00:30:00'), label_arcs: bool = True, normalize_by_dt: bool = True, elev_col: str | None = None, cutoff_deg: float | None = None, title: str | None = None, phase_to_meters: bool = False, wavelength_m: float | None = None)

Plot temporal derivative of a single-difference observable by PRN.

Notes (pedagogical)

• If obs is a carrier-phase SD_L* (stored in cycles), you may set phase_to_meters=True to convert cycles -> meters using the wavelength. This is ONLY for comparison with code (meters); it does not remove ambiguities.

param df_SD:

MultiIndex (epoch, prn). Must contain column obs.

type df_SD:

DataFrame

param obs:

Column to differentiate (e.g., “SD_L1” cycles, “SD_C1” meters).

type obs:

str

param phase_to_meters:

If True and obs starts with “SD_L”, convert cycles to meters before differencing.

type phase_to_meters:

bool

param wavelength_m:

Override wavelength (meters) used for cycles->meters conversion. If None and obs is SD_L1/SD_L2/SD_L5, use conv.L*_WAVELENGTH.

type wavelength_m:

float or None

geodezyx.gnss_edu.gnss_edu.plot_series(df, col1, col2=None, coeff1=1.0, coeff2=1.0, seuil=3600, renderer='browser')

Affiche les séries temporelles pour chaque satellite. Si col2 est fourni, affiche la série : coeff1 * col1 - coeff2 * col2. Sinon, affiche la série de la colonne col1 directement. La série est découpée en segments lorsqu’un “trou” (écart > seuil) est détecté.

Paramètres

dfDataFrame

DataFrame avec un index multi-niveaux contenant au moins le niveau ‘prn’ et les colonnes col1 (et éventuellement col2).

col1str

Nom de la première colonne.

col2str, optional

Nom de la seconde colonne (optionnel). Si None, on affiche col1.

coeff1float

Coefficient multiplicateur pour la première colonne (défaut 1.0).

coeff2float

Coefficient multiplicateur pour la seconde colonne (défaut 1.0).

seuilfloat

Seuil en secondes pour considérer un “trou” dans la série (défaut 3600).

rendererstr

Renderer Plotly (ex : “browser” ou “iframe”).

Retourne

figFigure

Figure Plotly contenant les courbes tracées.

geodezyx.gnss_edu.gnss_edu.plot_tracking_timeline(df: DataFrame, sampling: Timedelta = Timedelta('0 days 00:00:30'), snr_col: str | None = None, snr_min: float | None = None, title: str | None = None)

Timeline heatmap of tracking availability by PRN.

Black = epoch exists for this PRN (and SNR >= snr_min if provided) White = missing epoch (or SNR < snr_min)

Parameters:

• df (DataFrame) – Must be indexed by MultiIndex (‘epoch’,’prn’).

• sampling (Timedelta) – Expected observation interval (e.g. 30 s).

• snr_col (str or None) – Column name for SNR (e.g., ‘S1’). If None, only presence/absence is shown.

• snr_min (float or None) – If provided with snr_col, only epochs with SNR >= snr_min are considered “present”.

geodezyx.gnss_edu.gnss_edu.plot_tracking_timeline_with_pivots(df: DataFrame, selected_prns: list[str], sampling: Timedelta = Timedelta('0 days 00:00:30'), snr_col: str = 'S1', snr_min: float = 40.0, active_pivot: Series | None = None, title: str | None = None)

Tracking timeline with SNR threshold + pivot satellites highlighted.

Color meaning

0 = white : missing epoch OR SNR < snr_min 1 = black : usable (SNR >= snr_min) but NOT selected as pivot 2 = green : usable (SNR >= snr_min) AND selected as pivot 3 = dark green : active pivot at this epoch (provided by active_pivot)

param df:

MultiIndex (epoch, prn). Must contain snr_col.

type df:

DataFrame

param selected_prns:

PRNs selected by the greedy algorithm (candidate pivots).

type selected_prns:

list[str]

param sampling:

Expected sampling interval (e.g. 30 s).

type sampling:

Timedelta

param snr_col:

SNR column name (e.g., “S1”).

type snr_col:

str

param snr_min:

Minimum SNR threshold to consider a satellite usable.

type snr_min:

float

param active_pivot:

Series indexed by expected_epochs, values are PRN (string) or None/NaN. If provided, it is used to mark state=3 (dark green).

type active_pivot:

pd.Series | None

param title:

Plot title override.

type title:

str | None

rtype:

fig, ax, info

geodezyx.gnss_edu.gnss_edu.prn_covers_interval(df: DataFrame, prn: str, t0: Timestamp, t1: Timestamp, snr_col: str, snr_min: float, sampling: Timedelta) → bool

True if PRN has non-NaN SNR >= snr_min at ALL expected epochs in [t0, t1].

geodezyx.gnss_edu.gnss_edu.prn_snr_on_grid(df: DataFrame, prn: str, snr_col: str, epochs: DatetimeIndex) → Series

Return SNR series for one PRN reindexed on epochs. Missing epochs -> NaN.

geodezyx.gnss_edu.gpt3 module

Created on Fri Feb 9 07:43:47 2024

@author: snahmani

geodezyx.gnss_edu.gpt3.gpt3_5_fast(mjd=None, lat=None, lon=None, h_ell=None, it=None, grid=None)

geodezyx.gnss_edu.gpt3.gpt3_5_fast_readGrid(filename='./gpt3_5.grd')

geodezyx.gnss_edu.klobuchar module

Created on Tue Feb 6 10:30:57 2024

Filière ING3 - PPMD - Traitement de la mesure de phase

@author: Samuel Nahmani (1,2) https://www.ipgp.fr/annuaire/nahmani/) contact : nahmani@ipgp.fr ou samuel.nahmani@ign.fr (1) Université Paris Cité, Institut de physique du globe de Paris, CNRS, IGN, F-75005 Paris, France. (2) Univ Gustave Eiffel, ENSG, IGN, F-77455 Marne-la-Vallée, France.

Version: 1.0 Dépendances: numpy, geodezyx, datetime

geodezyx.gnss_edu.klobuchar.klobuchar(phi, lambda_, elev, azimuth, tow, alpha, beta)

Compute ionospheric range correction for GPS L1 frequency using Klobuchar model.

Translation to Python of Meysam Mahooti (2024). Klobuchar Ionospheric Delay Model https://www.mathworks.com/matlabcentral/fileexchange/59530-klobuchar-ionospheric-delay-model MATLAB Central File Exchange. Retrieved February 7, 2024.

Parameters:

• phi (float) – Geodetic latitude of receiver (degrees).

• lambda (float) – Geodetic longitude of receiver (degrees).

• elev (float) – Elevation angle of satellite (degrees).

• azimuth (float) – Geodetic azimuth of satellite (degrees).

• tow (float) – Time of Week (seconds).

• alpha (array_like) – The coefficients of a cubic equation representing the amplitude of the vertical delay (4 coefficients - 8 bits each).

• beta (array_like) – The coefficients of a cubic equation representing the period of the model (4 coefficients - 8 bits each).

Returns:

d_ion1 – Ionospheric slant range correction for the L1 frequency (metres).

Return type:

float

References

Klobuchar, J.A., (1996) “Ionospheric Effects on GPS”, in Parkinson, Spilker (ed), “Global Positioning System Theory and Applications, pp.513-514.

ICD-GPS-200, Rev. C, (1997), pp. 125-128

NATO, (1991), “Technical Characteristics of the NAVSTAR GPS”, pp. A-6-31 - A-6-33

geodezyx.gnss_edu.read_vmf1_grid module

read_vmf1_grid.py — autonome (sans jdutil) Lecture + interpolation des coefficients VMF1 (ah, aw) aux coordonnées station et à des dates MJD arbitraires, à partir des fichiers VMFG_YYYYMMDD.HHH.

• Interpolation spatiale:

  • tente scipy.interpolate.griddata (linear)

  • sinon repli nearest neighbor (NumPy) si SciPy indisponible

• Interpolation temporelle: linéaire entre pas 6h

• Chemins: ./Mapping_Fcn/vmf1/VMFG_YYYYMMDD.HHH

Format attendu des fichiers (après les lignes commençant par ‘!’):

lat lon ah aw zhd zwd (au minimum: lat, lon, ah, aw)

Auteur: adapté pour être autonome (remplacement de jdutil)

class geodezyx.gnss_edu.read_vmf1_grid.VMFPoint(lat: 'float', lon: 'float', ah: 'float', aw: 'float')

Bases: object

ah: float

aw: float

lat: float

lon: float

geodezyx.gnss_edu.read_vmf1_grid.add_hours(dt, hours: int)

geodezyx.gnss_edu.read_vmf1_grid.datetime_to_mjd(dt) → float

Convertit un datetime UTC naïf/aware en MJD (UTC).

geodezyx.gnss_edu.read_vmf1_grid.floor_to_6h(dt)

Ramène dt à l’heure multiple de 6h immédiatement inférieure (UTC).

geodezyx.gnss_edu.read_vmf1_grid.jd_to_calendar(jd: float) → Tuple[int, int, int, float]

Convertit un Julian Day (float) en (year, month, day, hour_float). Implémentation autonome (remplace jdutil.jd_to_date).

geodezyx.gnss_edu.read_vmf1_grid.mjd_to_datetime(mjd: float)

Conversion MJD -> datetime UTC naïf (sans tzinfo).

geodezyx.gnss_edu.read_vmf1_grid.mjd_to_jd(mjd: float) → float

geodezyx.gnss_edu.read_vmf1_grid.read_vmf1_grid(mjd: Iterable[float], ell: Iterable[float]) → Tuple[ndarray, ndarray]

Retourne (ah_vec, aw_vec) pour chaque MJD donné et une station ell=[lat_deg, lon_deg, h_m].

Paramètres

mjdIterable[float]

Liste/array de dates en Modified Julian Day (UTC).

ell[lat_deg, lon_deg, h_m]

Latitude (deg), Longitude (deg), Hauteur ellipsoïdale (m). (La hauteur n’est pas utilisée ici pour ah/aw, mais on garde la signature d’interface.)

Renvoie

ah_vec : np.ndarray shape (N,) aw_vec : np.ndarray shape (N,)

Notes

• Les fichiers requis sont : ./Mapping_Fcn/vmf1/VMFG_YYYYMMDD.H{00,06,12,18} pour les 6h entourant chaque MJD.

• Interpolation temporelle: linéaire entre t0=⌊MJD⌋_6h et t1=t0+6h.