PyTsfit

PyTsfit fits position time series of GNSS/GPS stations. It estimates secular velocities, seasonal (annual / semi-annual) signals, coseismic offsets, non-earthquake breaks and postseismic transients from PBO .pos or GAMIT/GLOBK .neu time series, and provides robust uncertainty estimates via realistic-sigma error scaling, iterative outlier editing and MCMC sampling.

Features

• Trajectory model — constant offset, linear (secular velocity), annual and semi-annual seasonal terms, earthquake offsets, non-earthquake breaks and postseismic transients.
• Input formats — PBO .pos files (posData) and GAMIT/GLOBK .neu files (neuData).
• Prior constraints — optional prior secular velocity, coseismic displacement, non-earthquake break and seasonal (e.g. GRACE-derived) values are subtracted before fitting via the correction class.
• Quality control — two mechanisms ported from GAMIT/GLOBK tsfit:
  • realistic-sigma error scaling (real_stats → realistic_sigma): Herring bin-averaging plus FOGM extrapolation, giving realistic uncertainties that account for time-correlated noise;
  • iterative outlier editing (edit_ns → flag_outliers): post-fit n-sigma rejection with MAD-based robust scale and hysteresis.
• Outputs — secular velocities, coseismic offsets, breaks, postseismic displacements and postseismic time series, plus observed/modeled figures and parameter tables. All outputs are GMT-compatible (comment-header columns).
• MCMC uncertainty — emcee-based sampling of the trajectory parameters (coseismic and postseismic variants), including MPI-parallel chains via schwimmbad, and corner plots of the posterior.

Installation

pip install -e .

Dependencies (installed automatically): numpy, scipy, matplotlib, pandas, pyproj, scikit-learn, emcee, schwimmbad, corner, plotly, pyyaml.

Requires Python >= 3.7.

Quick start

The main entry point is the do_pytsfit console script, driven by a YAML configuration file:

do_pytsfit --cfgfile config.yaml

A full example configuration lives at src/pytsfit/scripts/config.yaml. Command-line options override individual YAML parameters; the YAML file remains the single source of defaults. Options not given on the command line fall back to the values in the config file, e.g.:

# Override a few parameters for a single run
do_pytsfit --cfgfile config.yaml --sitelist BJFS XIAJ --tsdir ../pos/ \
           --timespan 2000 2024 --annual true --outlier true --nsigma 3.0

--sitelist takes one or more site names and overrides the sitefile in the configuration.

Interactive UI (Streamlit)

An interactive web interface wraps the same Python code as the CLI. It lets you browse stations on a map, inspect the raw N/E/U time series, adjust the trajectory model and fit options, run the least-squares fit and view the observed/modeled plots, residuals, parameter tables and the same output tables the CLI writes (velocity, eq offsets, breaks), all in the browser:

streamlit run src/pytsfit/ui/app.py

The UI is a thin shell over the existing API (pytsfit/ui/api.py); no fitting or plotting logic is re-implemented. All charts are interactive Plotly figures (no Mapbox token required). The sidebar mirrors the YAML config — dict_input, dict_param, dict_fit, dict_plot and the output toggles. Fit results are kept per session so you can switch parameters without re-fitting everything.

Configuration

The YAML configuration has five sections:

dict_input:                # input files and time series
    eqfile           : './eq_rename'
    prior_velfile    : ''              # prior secular velocity
    prior_offsetfile : ''              # prior coseismic offsets
    prior_periodfile : ''              # prior seasonal terms
    sitefile         : 'cmonoc.cmnc'   # site list (one name per line)
    timespan         : [1998, 2024]    # [start, end] in decimal years
    tsdir            : '../pos/'       # directory holding the time series
    tsformat         : 'pos'           # 'pos' or 'neu'
dict_param:              # which terms to estimate
    constant    : True
    linear      : True                 # secular velocity
    annual      : True
    semiannual  : True
    break       : True                 # non-earthquake breaks from eqfile
    eqoffset_ne : True                 # coseismic offsets, N/E
    eqoffset_up : True                 # coseismic offsets, U
    eqpost_ne   : False                # postseismic transients, N/E
    eqpost_up   : False                # postseismic transients, U
dict_plot:               # plotting options
    detrend    : True
    debreak    : True
    deeqoffset : False
    depost     : False
    deseason   : False
    showfig    : False
    figformat  : 'jpg'
dict_fit:                # fitting / quality-control options
    sigma_scale    : 'nrms'   # none | nrms | realistic
    min_rsig       : 30       # minimum used points to attempt FOGM
    min_sigscale   : 10       # minimum dof below which sig_scale = 1
    outlier        : False    # enable iterative outlier editing
    nsigma         : 4.0      # rejection threshold
    outlier_scale  : 'mad'    # 'mad' (robust) or 'sigma'
    max_iter       : 10       # maximum edit/refit passes
    restore_factor : 0.9      # hysteresis band for restoring points
    max_sigma      : 0.0      # pre-fit sigma screen (mm); 0 = off
dict_output:             # what to write and where
    tsfig       : False       # observed/modeled figure per site
    param       : False       # parameter table
    obsmod      : False       # observed and modeled series file
    eqpostts    : False       # postseismic time series files (obs/mod)
    velfile     : ''          # velocity output file
    eqoffset    : ''          # coseismic offset output file
    break       : ''          # break output file
    eqpostdisp  : ''          # postseismic displacement output file
    eqpost_tspan: [2015, 2023]

Fit options (dict_fit)

sigma_scale controls how the covariance returned by the least-squares fit is scaled:

• none — use the raw covariance, no scaling;
• nrms — scale by chi2/dof (scipy default behaviour);
• realistic — Herring bin-averaging plus FOGM extrapolation (real_stats), giving uncertainties that account for time-correlated noise.

outlier enables iterative post-fit n-sigma editing with the robust MAD scale. The exact semantics of each option are documented in pytsfit.qualitycontrol.DEFAULT_FIT_OPTS.

Input files

Time series

• PBO .pos — PBO-style position time series, columns decyr N E U SN SE SU (read by posData).
• GAMIT/GLOBK .neu — same columns (read by neuData).

Files are matched by glob.glob('{tsdir}/{site}*.{tsformat}'), so several time series per site are allowed.

Earthquake / break file (eqfile)

A catalog of earthquakes (and/or breaks) listing each event's code, location, epoch and the affected stations. Used by eqcatalog (coseismic offsets), breakcatalog (non-earthquake breaks) and eqPostList (postseismic events).

Prior information files

Optional priors are subtracted before fitting and given to the correction class:

• prior_velfile — prior secular velocities: Lon, Lat, Ve, Vn, Sig_ve, Sig_vn, Cor_en, Site, Vu, Sig_vu (units mm/yr).
• prior_offsetfile — prior coseismic displacements: E, N, U, Site, decimal-year (units mm).
• prior_periodfile — prior seasonal terms (e.g. GRACE-derived): E_sa, E_ca, E_ssa, E_csa, N_sa, N_ca, N_ssa, N_csa, U_sa, U_ca, U_ssa, U_csa, Site where sa/ca are the annual sine/cosine amplitudes and ssa/csa the semi-annual ones (units mm). Sites missing from this file are assigned zeros.

Outputs

Output files are opened in append mode and written as comment-header columns so they can be consumed directly by GMT or np.genfromtxt:

• velocity — Lon, Lat, Ve, Vn, Se, Sn, Corr, Vu, Su, Site, plus time-span statistics.
• coseismic offsets — Lon, Lat, E/N/U offsets and uncertainties, EQ code.
• breaks — non-earthquake break offsets.
• postseismic displacement — E/N/U displacement and WRMS over eqpost_tspan.

MCMC scripts

In addition to the fast least-squares fit, the package provides emcee-based sampling of the trajectory parameters:

Script	Purpose
do_pytsfit_mcmc_coseismic.py	MCMC sampling of coseismic offsets / velocities / seasonal terms
do_pytsfit_mcmc_postseismic.py	MCMC sampling of postseismic transients; writes postseismic.gmtvec
tsfit_mcmc.py	MCMC fitting (single process)
tsfit_mcmc_mpi.py / tsfit_mcmc_mpi_new.py	MPI-parallel MCMC fitting via schwimmbad
plot_corner.py	corner plot of the chains saved in chain.npz

These live in src/pytsfit/scripts/ and read the same config.yaml layout. plot_corner.py can be run directly on the chain.npz produced by the MCMC scripts.

Package layout

The code is organized by responsibility in src/pytsfit/:

Module	Responsibility
data.py	time-series readers (neuData, posData)
models.py	events and priors (earthquake, eqcatalog, eqPost, eqPostList, offset, breakcatalog, correction)
tsfitting.py	the fitting engine
qualitycontrol.py	realistic-sigma error scaling and iterative outlier editing
output.py	output_* writers and plot_obs_mod
geotools.py, GPSTime.py	geodetic coordinate and time utilities
PyTsfit.py	compatibility shim re-exporting the public API (incl. build_param_dict)
ui/	Streamlit front-end (app.py) and thin wrappers (api.py)
scripts/	do_pytsfit entry point and auxiliary/MCMC scripts

The public API is re-exported from pytsfit.PyTsfit, e.g.:

from pytsfit.PyTsfit import posData, tsfitting, output_velo

Testing

pytest

The suite covers the CLI/YAML override logic, quality-control functions (realistic sigma, outlier editing) and characterization tests with golden baselines.

License

MIT (see pyproject.toml).

Author

Bin Zhao — zhaobin@cgps.ac.cn, Institute of Seismology, CEA.