pyacs.gts.Gts module

Individual geodetic time series (Gts) class and helpers.

The Gts class holds NEU (or XYZ) time series with metadata. Key attributes: data (2D: dec_year, N, E, U, S_N, S_E, S_U, S_NE, S_NU, S_EU), code, lon/lat/h, X0,Y0,Z0,t0, data_xyz, outliers, offsets_*, velocity, annual, semi_annual, ifile, log, metadata. Units: dates in decimal year, coordinates in m, phases in radians.

class pyacs.gts.Gts.Gts(code=None, lat=None, lon=None, h=None, X0=None, Y0=None, Z0=None, t0=None, data=None, data_xyz=None, data_corr_neu=None, data_corr_xyz=None, offsets_dates=[], offsets_values=None, outliers=[], annual=None, semi_annual=None, velocity=None, ifile=None, log=None, metadata=None)[source]

Bases: object

add_obs(date, NEUSNSESUCNECNUCEU, in_place=False, check=True, verbose=False)

Add observation(s) as DN, DE, DU to a time series.

Parameters:

date (float or list or ndarray) – Date(s) in decimal year.
NEUSNSESUCNECNUCEU (list or ndarray) – Values to add: at least NEU (North, East, Up). Optional: SN, SE, SU, CNE, CNU, CEU (standard deviations and correlations). If not provided, SN=SE=SU=0.001, CNE=CNU=CEU=0.
in_place (bool, optional) – If True, add to current Gts; if False, return a new Gts.
check (bool, optional) – Check time order and duplicate dates.
verbose (bool, optional) – Verbose mode.

Returns:

New Gts or the modified Gts if in_place.

Return type:

Notes

If .data_xyz exists, it will be set to None for consistency.

add_obs_xyz(date, XYZSXSYSZCXYCXZCYZ, in_place=False, check=True, neu=True, verbose=False)

Add observation(s) as XYZ to a time series.

Parameters:

date (float or list or ndarray) – Date(s) in decimal year.
XYZSXSYSZCXYCXZCYZ (list or ndarray) – Values to add: at least X, Y, Z. Optional: SX, SY, SZ, CXY, CXZ, CYZ (standard deviations and correlations). If not provided, SX=SY=SZ=0.001, CXY=CXZ=CYZ=0.
in_place (bool, optional) – If True, add to current Gts; if False, return a new Gts.
check (bool, optional) – Check time order, duplicate dates, and regenerate NEU (.data).
neu (bool, optional) – Regenerate .data from updated .data_xyz.
verbose (bool, optional) – Verbose mode.

Returns:

New Gts or the modified Gts if in_place.

Return type:

Notes

By default .data is updated from .data_xyz and X0, Y0, Z0 are updated.

add_offsets_dates(offsets_dates, in_place=False): add_offsets_dates to a time series if in_place = True then replace the current time series

add_vel_sigma(b_fn=4, verbose=True)

Estimate velocity uncertainties from residuals (white + flicker noise).

Velocity and residual time series must exist: run detrend(), detrend_annual(), or detrend_seasonal() first. Uncertainty on velocity components (N, E, U) is computed using white noise from the residual scatter and flicker noise (Williams, 2003, eq. 19 and 23). Returns a new Gts with velocity sigmas set.

Parameters:

b_fn (float, optional) – Flicker noise parameter in mm/yr^0.25. Default is 4.
verbose (bool, optional) – If True, print flicker noise variance. Default is True.

Returns:

New Gts with velocity uncertainties set in velocity[3:6] (m/yr), or None if self.velocity is not set.

Return type:

Gts or None

Notes

White noise is estimated from the first differences of the residual time series. The combined white + flicker sigma is scaled by component RMS so that the minimum component (N, E, or U) keeps the nominal flicker level.

References

Williams, S. D. P. (2003). The effect of coloured noise on the uncertainties of rates from geodetic time series. Journal of Geodesy, 76(9-10), 483-494.

apply_offsets(np_offset, opposite=False, in_place=False, verbose=False)

Apply given offsets to the time series.

Parameters:

np_offset (numpy.ndarray or list) – 1D or 2D array (or list of lists): offset_dates, north, east, up, s_north, s_east, s_up.
opposite (bool, optional) – If True, apply opposite of offsets. Default is False.
in_place (bool, optional) – If True, modify self; else return new Gts. Default is False.
verbose (bool, optional) – Verbose mode. Default is False.

Returns:

self (if in_place) or new Gts.

Return type:

cdata(data=False, data_xyz=False, tol=0.001, verbose=False)

Check data and/or data_xyz attributes for consistency and validity.

Parameters:

data (bool, optional) – If True, check the data attribute. Default is False.
data_xyz (bool, optional) – If True, check the data_xyz attribute. Default is False.
tol (float, optional) – Tolerance in days for two dates to be considered the same. Default is 0.001.
verbose (bool, optional) – If True, print details. Default is False.

Returns:

True if checks pass, False otherwise.

Return type:

bool

Notes

In future, this routine may also check consistency between .data and .data_xyz.

check_l1_trend(l1ts, component='ENU', min_samples_per_segment=4, threshold_bias_res_detect=40, threshold_vel=8, plot=False)

Inspect the result from l1trend model of a time series. Returns a list periods ([sdate,edate]) where bad modelling is suspected.

Parameters:

ts (pyacs.gts.Gts.Gts) – Raw time series
l1ts (pyacs.gts.Gts.Gts) – L1-trend time series
component (str) – Components to be analyzed. Default: ‘EN’
min_samples_per_segment (int) – Minimum number of samples for a segment to be inspected (default: 6)
threshold_bias_res_detect (float) – Threshold to detect bias residuals (default: 40)
threshold_vel (float) – Instantaneous velocity threshold for a segment to be considered in detection (default: 8 mm/yr)
plot (bool) – Plot the results (default: False)

Returns:

(H_period, H_cp, H_cp_pb) - H_period: Dictionary of suspicious periods for each component - H_cp: Dictionary of breakpoint indices for each component - H_cp_pb: Dictionary of problematic breakpoint indices for each component

Return type:

tuple

clean_l1trend(raw_gts, threshold='auto')

Cleans breakpoints from an already L1-trend filtered time series. Removes breakpoints that are too close in value (<0.5 mm/yr). Returns a new Gts object interpolated on the dates from l1trend_gts.

Parameters:

raw_gts (pyacs.gts.Gts.Gts) – Original (position) time series.
threshold (float or 'auto') – Threshold for cleaning. If ‘auto’, computed from median difference.

Returns:

Cleaned/interpolated Gts object.

Return type:

pyacs.gts.Gts.Gts

copy(data=True, data_xyz=True, loutliers=True)

Return a deep copy of the time series.

By default copies .data, .data_xyz, outliers, etc. Behaviour can be overridden per attribute below.

Parameters:

data (bool or numpy.ndarray, optional) – True = copy .data; False = set to None; or (n,10) array. Default is True.
data_xyz (bool or numpy.ndarray, optional) – True = copy .data_xyz; False = None; or (n,10) array. Default is True.
loutliers (bool, optional) – If False, do not copy outliers. Default is True.

Returns:

New Gts instance.

Return type:

correct_duplicated_dates(action='correct', tol=0.1, in_place=False, verbose=False)

Check or remove duplicated dates in a time series.

Parameters:

action (str, optional) – ‘correct’ (default) or ‘check’.
tol (float, optional) – Tolerance in days for two dates to be considered the same (default 0.1 day).
in_place (bool, optional) – If True, modify in place; otherwise return a new Gts.
verbose (bool, optional) – Verbose mode.

Returns:

Current Gts if in_place, else new Gts with duplicates removed (when action=’correct’).

Return type:

decimate(time_step=30.0, dates=[], method='median', verbose=False)

Decimate a time series to a regular time step.

Parameters:

time_step (float, optional) – Time step in days.
dates (list, optional) – List of dates where points are forced to be written regardless of time_step.
method (str, optional) – Method to compute position in each bin: ‘median’, ‘mean’, or ‘exact’.
verbose (bool, optional) – Verbose mode.

Returns:

New decimated Gts.

Return type:

decyear2days(ref_date='', in_place=False)

Convert time series dates from decimal year to days after a reference date.

Parameters:

ref_date (float or str, optional) – Reference date (decimal year or parseable by guess_date). Default first data date.
in_place (bool, optional) – If True, modify in place; otherwise return a new Gts.

Returns:

Time series with dates in days (or self if in_place).

Return type:

delete_small_offsets(offsets, del_by_pricise=False)

Estimate offsets with clean data, then remove offsets that are too small.

Parameters:

offsets (list) – List of offset dates (decimal year).
del_by_pricise (bool, optional) – If True, also require offset > sigma. Default is False.

Returns:

Filtered list of offset dates, or None.

Return type:

list or None

detrend(method='L2', periods=[], exclude_periods=[], log_dir=None)

Detrend time series and store velocity (and optional offsets) in attributes.

Parameters:

method (str, optional) – Estimation method (e.g. ‘L2’). Default is ‘L2’.
periods (list, optional) – Periods used for velocity estimate. Default is [].
exclude_periods (list, optional) – Periods excluded from velocity estimate. Default is [].
log_dir (str, optional) – Directory for Gts-specific log file {self.code}.log. Default is None.

Returns:

Detrended time series (velocity attribute set).

Return type:

Notes

Outliers in Gts.outliers are omitted; offsets in Gts.offsets_dates are estimated simultaneously.

detrend_annual(method='L2', periods=None, exclude_periods=None)

Estimate trend + annual terms and remove them from the time series.

Velocity and annual are saved in Gts.velocity and Gts.annual.

Parameters:

method (str, optional) – Estimation method (e.g. ‘L2’).
periods (list, optional) – Periods used for estimation.
exclude_periods (list, optional) – Periods to exclude from estimation.

Returns:

Detrended time series (new instance).

Return type:

Notes

Outliers (Gts.outliers) are omitted in the estimation; offsets (Gts.offsets_dates) are estimated simultaneously.

detrend_hectorp(component='NEU', **kwargs)

Detrend a time series using the Hector model (HECTOR from Machiel Bos).

References: Bos et al. (2013), J. Geodesy 87(4), 351-360. See https://gitlab.com/machielsimonbos/hectorp and estimatetrend wiki.

Parameters:

component (str, optional) – Component(s) to detrend (‘N’, ‘E’, ‘U’ or ‘NEU’). Default is ‘NEU’.
**kwargs (dict, optional) – Additional arguments for the Hector control file (estimatetrend).

Returns:

Residual time series with velocity (and sigmas), offsets attributes.

Return type:

detrend_median(delta_day=None, periods=[], exclude_periods=[], verbose=False, auto=False, log_dir=None)

Compute velocity from median of 1-year pair displacements (MIDAS-style).

If the time series is shorter than one year, it is left unchanged.

Parameters:

delta_day (int or None, optional) – Time delta in days for pairs. None = one year, 0 = relax mode for campaign data. Default is None.
periods (list, optional) – Periods to include for trend. Default is [].
exclude_periods (list, optional) – Periods to exclude. Default is [].
verbose (bool, optional) – If True, print progress. Default is False.
auto (bool, optional) – If True, try delta_day=None then 0, else fallback to regular detrend. Default is False.
log_dir (str, optional) – Directory for log file {self.code}_detrend_median.log. Default is None.

Returns:

Detrended Gts, or None if series shorter than 1 year.

Return type:

Gts or None

detrend_pytrf(noise=['wn', 'fn'], method='Nelder-Mead', fixed_correlated_noise_value=[None, None], log_dir=None)

Detrend using PyTRF; estimate velocity and offsets with realistic sigmas.

Wrapper around pytrf (https://github.com/prebischung/pytrf, P. Rebischung). Simplified trajectory: constant velocity, offsets, annual/semi-annual, white + power-law noise. For more complex models use pytrf directly. Power-law noise can be fixed (useful for campaign data).

Parameters:

noise (list, optional) – Noise model: ‘wn’ (mandatory) plus ‘fn’, ‘rw’, or ‘pl’. Default is [‘wn’,’fn’].
method (str, optional) – Optimization: ‘Nelder-Mead’, ‘BFGS’, ‘CG’, ‘Newton’, ‘Powell’. Default is ‘Nelder-Mead’.
fixed_correlated_noise_value (list, optional) – [s2, a] for power-law; None to estimate. Default is [None, None].
log_dir (str, optional) – Directory for pytrf yaml log; None = current dir. Default is None.

Returns:

Residual time series with velocity (and sigmas), offsets.

Return type:

Notes

‘Nelder-Mead’ is default and often more robust than ‘BFGS’ despite being slower.

detrend_seasonal(method='L2', periods=None, exclude_periods=None)

Estimate trend + annual + semi-annual terms and remove them from the time series.

Velocity, annual and semi_annual are saved in Gts.velocity, Gts.annual, Gts.semi_annual.

Parameters:

method (str, optional) – Estimation method (e.g. ‘L2’).
periods (list, optional) – Periods used for estimation.
exclude_periods (list, optional) – Periods to exclude from estimation.

Returns:

Detrended time series (new instance).

Return type:

Notes

Outliers (Gts.outliers) are omitted; offsets (Gts.offsets_dates) are estimated simultaneously.

detrend_seasonal_median(wl=11, verbose=False): Calculates a velocity using the median of pair of displacements exactly separated by one year, inspired from MIDAS and then removes repeating yearly signal If the time series has less than three years of data, then the time series is kept untouched.

differentiate()

Differentiate the current time series (first difference along time).

Returns:: Differentiated time series as a new Gts object.
Return type:: Gts

Notes

Differentiation is performed on .data. .data_xyz is set to None.

disp2vel(alpha='auto', ndays_mf=1)

Compute time series derivative via L1 trend filter (BIC-optimal parameter).

Uses https://pypi.org/project/trendfilter/

Parameters:

alpha (str or float, optional) – Filter parameter; ‘auto’ for BIC-optimal. Default is ‘auto’.
ndays_mf (int, optional) – Median filter window (days) before L1-trend. Default is 1.

Returns:

New Gts instance; units m/yr.

Return type:

displacement(sdate=None, edate=None, window=None, method='median', speriod=[], eperiod=[], rounding='day', verbose=True)

Calculate displacement between two dates or two periods.

Parameters:

sdate (float, optional) – Start date in decimal year.
edate (float, optional) – End date in decimal year.
window (float, optional) – Time window in days for searching available dates.
method (str, optional) – Method to calculate position: ‘median’ or ‘mean’. Default is ‘median’.
speriod (list, optional) – Period [start, end] for calculating the start position.
eperiod (list, optional) – Period [start, end] for calculating the end position.
rounding (str, optional) – Rounding for dates: ‘day’, ‘hour’, ‘minute’, or ‘second’. Default is ‘day’.
verbose (bool, optional) – Verbose mode.

Returns:

Displacement as np.array([dn, de, du, sdn, sde, sdu]).

Return type:

ndarray

edge_filter(alpha='auto', ndays_mf=1)

L1 edge filter.

Uses https://pypi.org/project/trendfilter/

Parameters:

alpha (str or float, optional) – Filter parameter; ‘auto’ for BIC-optimal. Default is ‘auto’.
ndays_mf (int, optional) – Median filter window (days) applied before L1 edge filter. Default is 1 (no median filter).

Returns:

L1 edge filter as a new Gts instance.

Return type:

estimate_local_offset(window_length=4, in_place=False)

Estimate local offset amplitudes using window_length positions before and after each offset.

Returns:: New Gts (or self if in_place) with offsets_values set.
Return type:: Gts

exclude_periods(lperiod, in_place=False, verbose=False)

Exclude periods of a Gts.

Parameters:

lperiod (list) – A list [start_date,end_date] or a list of periods e.g. [[2000.1,2003.5],[2009.3,2010.8]].
in_place (bool, optional) – If True, will make change in place; if False, returns a new time series.
verbose (bool, optional) – Verbose mode.

Notes

X0, Y0, Z0 attributes will be changed if necessary.
Handles both .data and .data_xyz.

extract_dates(dates, tol=0.05, rounding=None, in_place=False, verbose=True)

Return a time series extracted for a given list of dates.

Parameters:

dates (list or ndarray) – Dates as a list or 1D numpy array of decimal dates.
tol (float, optional) – Date tolerance in days to assert that two dates are equal (default 0.05 day).
rounding (str, optional) – Rounding for date comparison (‘second’, ‘minute’, ‘hour’, ‘day’). Inferred from tol if None.
in_place (bool, optional) – If True, will make change in place; if False, returns a new time series.
verbose (bool, optional) – Verbose mode.

extract_ndates_after_date(date, n, verbose=False)

Extract n values after a given date.

If n values are not available, returns all available values after date. .data is set to None if no value at all is available.

Parameters:

date (float) – Date in decimal year.
n (int) – Number of observations to extract.
verbose (bool, optional) – Verbose mode.

Returns:

New Gts with at most n observations after the given date.

Return type:

extract_ndates_around_date(date, n)

Extract n values before and n values after a given date.

If n values are not available, returns all available values. .data is set to None if no value at all is available.

Parameters:

date (float) – Date in decimal year.
n (int) – Number of observations to extract on each side.

Returns:

New Gts with at most n observations before and n after the given date.

Return type:

extract_ndates_before_date(date, n, verbose=False)

Extract n values before a given date.

If n values are not available, returns all available values before date. .data is set to None if no value at all is available.

Parameters:

date (float) – Date in decimal year.
n (int) – Number of observations to extract.
verbose (bool, optional) – Verbose mode.

Returns:

New Gts with at most n observations before the given date.

Return type:

extract_periods(lperiod, in_place=False, verbose=False, no_reset=False, ignore_data_xyz=False)

Extract periods of a Gts.

Parameters:

lperiod (list) – A list [start_date,end_date] or a list of periods e.g. [[2000.1,2003.5],[2009.3,2010.8]].
in_place (bool, optional) – If True, will make change in place; if False, returns a new time series.
verbose (bool, optional) – Verbose mode.
no_reset (bool, optional) – If True, do not reset X0, Y0, Z0 to first epoch of extracted data.
ignore_data_xyz (bool, optional) – If True, work on .data only and ignore .data_xyz.

Notes

X0, Y0, Z0 attributes will be changed if necessary.
Handles both .data and .data_xyz.

extrapolate(idates, kind='linear')

Extrapolate the time series at dates outside the observation span.

Parameters:

idates (list or ndarray) – Dates in decimal year where extrapolation will be performed.
kind (str, optional) – Interpolation kind for scipy.interpolate.interp1d (e.g. ‘linear’, ‘nearest’, ‘cubic’).

Returns:

New Gts with extrapolated values at the given dates.

Return type:

find_large_uncertainty(sigma_threshold=20, verbose=True, lcomponent='NE', log_dir=None)

Flag dates with large uncertainty as outliers.

Parameters:

sigma_threshold (float, optional) – Threshold in mm for flagging. Default is 20.
verbose (bool, optional) – If True, print progress. Default is True.
lcomponent (str or list, optional) – Components to check (‘N’, ‘E’, ‘U’). Default is ‘NE’.
log_dir (str, optional) – Directory for Gts-specific log file {self.code}.log. Default is None.

Returns:

self (outliers updated).

Return type:

find_offsets(threshold=3, n_max_offsets=9, conf_level=95, lcomponent='NE', verbose=True, in_place=False)

Simple empirical procedure to find offsets.

Parameters:

threshold (float, optional) – Threshold for preliminary offset detection. Default is 3.
n_max_offsets (int, optional) – Maximum number of offsets to detect. Default is 9.
conf_level (float, optional) – Confidence level (percent) to accept offset. Default is 95.
lcomponent (str, optional) – Components for detection (‘N’,’E’,’U’). Default is ‘NE’.
verbose (bool, optional) – Verbose mode. Default is True.
in_place (bool, optional) – If True, modify self. Default is False.

Returns:

New Gts (or self if in_place) with offsets_dates and outliers set.

Return type:

find_offsets_edge_filter(threshold=0.6, search_lbda=[3, 5, 7, 10, 20, 50, 100, 200, 300], delta_day=100, in_place=False, lcomponent='NE', verbose=True, debug=True, log=False, eq_file=None)

find_offsets_ivel(ivelts=None, lcomponent='EN', threshold_offset=100, offsets_file=None)

Find offsets using l1trend filtering and instantaneous velocity

parameters:

ivelts: Sgts: ivel time series
lcomponent: string: component to use for the analysis
threshold_offset: float: threshold for the offset in mm/yr
offsets_file: string: name of the file to append the offsets. If None, no file is written

notes:

The function uses the l1trend of the time series. It then detects large ivel values and check whether they last for more than one day. If so, an offset is estimated by integrating ivel over two successive dates and offset is applied to the time series. Format of offsets_file is: code date (datetime.isoformat()) lon lat offset_e offset_n offset_u std_offset_e std_offset_n std_offset_u comment

find_offsets_t_scan(threshold=0.8, window=250, in_place=False, lcomponent='NE', verbose=True, debug=True): Find suspected offsets using t-scan step detection.

find_outlier_around_date(date, conf_level=95, n=3, lcomponent='NE', verbose=True)

Find an outlier around a given date using F-ratio test.

Returns the index of the outlier, or [] if none found.

Parameters:

date (float) – Date in decimal year.
conf_level (float, optional) – Confidence level for F-ratio test (default 95).
n (int, optional) – Number of dates on each side of date (default 3).
lcomponent (str, optional) – Components to test: ‘N’, ‘E’, ‘U’, ‘NE’, ‘NEU’ (default ‘NE’).
verbose (bool, optional) – Verbose mode.

Returns:

self with outlier flagged, or [] if no significant outlier.

Return type:

Gts or list

find_outliers_l1trend(lam, threshold, period=None, gap=10, components='NE', plot=False, verbose=False, in_place=False)

Find outliers using L1 trend filtering and residual threshold.

Parameters:

lam (float) – Lambda parameter for L1 trend filtering.
threshold (float) – Residuals above threshold * (scale) are flagged as outliers.
period (list or None, optional) – Period(s) for searching outliers.
gap (int or float, optional) – Gap in days to split series (default 10).
components (str, optional) – Components for detection (default ‘NE’).
plot (bool, optional) – If True, plot filter result and flagged outliers.
verbose (bool, optional) – Verbose mode.
in_place (bool, optional) – If True, modify current Gts; otherwise return a new Gts.

Returns:

Gts with outliers set (new instance or self if in_place).

Return type:

find_outliers_percentage(percentage=0.03, in_place=False, verbose=False, component='NEU', periods=None, excluded_periods=None): detrend a time series and ranks the residuals by increasing absolute value populate the outliers with the x % largest ones on each component

find_outliers_simple(threshold=100, window_length=10, in_place=False, verbose=False, component='NEU', periods=None, excluded_periods=None)

find_outliers_sliding_window(threshold=3, in_place=False, verbose=True, periods=[[]], excluded_periods=[[]], component='NE', window_len=15, automatic=True): Find outliers using sliding windows

find_outliers_vondrak(threshold=10, fc=2.0, in_place=False, verbose=True, periods=[[]], excluded_periods=[[]], component='NE'): Find outliers using a Vondrak filter

find_time_offsets(option=None, ndays=7, th_detection_rms=3, th_detection_offset=3): Find the time of suspected offsets by RMS time series calculated over ndays.

flag_outliers_using_l1trend(l1trend, threshold=5)

Flag outliers in the time series based on a user provided l1trend filter representation of the time series.

This function identifies outliers by comparing the original time series to its L1-trend representation. Outliers are flagged based on their deviation from the trend using the Median Absolute Deviation (MAD) method.

Parameters:

self (Gts object) – The time series to flag outliers in.
l1trend (Gts object) – The l1trend representation of the time series.
threshold (float, optional) – The threshold for flagging outliers. Default is 5 for 5 times the median absolute deviation (MAD) from the median.

Returns:

The original time series object with outliers flagged in self.outliers

Return type:

Notes

The method uses the Median Absolute Deviation (MAD) which is more robust to outliers than standard deviation. Outliers are identified as points that deviate more than threshold * MAD from the median of the absolute differences between the original time series and the L1-trend representation.

force_daily(in_place=False): force a time series to be daily with dates at 12:00:00 of every day

frame(frame=None, verbose=False): Rotates a time series according to an Euler pole. Returns a new Gts instance.

get_coseismic(eq_date, window_days=5, sample_after=1, method='median', in_place=False)

Get coseismic displacement at a given date. Coseismic displacement is estimated as the position difference between the median of window_days before the earthquake date and the median of sample_after samples after the earthquake date.

note: only median method implemented

get_coseismic_l1trend(date)

Get coseismic offsets using pyacs customized l1 trend filter

param: date: date in decimal year, or datetime.datetime, or (mday,month,year)

get_unr(site, verbose=False)

Fetch a time series from UNR (IGS20 txyz) and return as Gts.

Parameters:

site (str) – 4-letter site code.
verbose (bool, optional) – If True, print progress. Default is False.

Returns:

Loaded time series or None on failure.

Return type:

Gts or None

get_unr_loading(site, verbose=False)

Get NTAL, NTOL, HYDL, MASC, RESE loading predictions from UNR.

See http://geodesy.unr.edu/gps_timeseries/README_tenv3load.txt for format.

Parameters:

site (str) – Station 4-letter code.
verbose (bool, optional) – Verbose mode.

Returns:

(NTAL, NTOL, HYDL, MASC, RESE) loading time series.

Return type:

tuple of 5 Gts

Examples

NTAL_CSEC, NTOL_CSEC, HYDL_CSEC, MASC_CSEC, RESE_CSEC = Sgts().get_unr_loading(‘CSEC’)

Notes

UNR .tenv3 files include load prediction triplets (NTAL, NTOL, HYDL, MASC, RESE, WUS, WUSA). Columns 21-23: NTAL; 24-26: NTOL; 27-29: HYDL; 30-32: MASC; 33-35: RESE. Prior to 2002 GRACE mascon data are unavailable (extrapolation used).

get_values_at_date(date, data_type='data_xyz')

return the values of the time series at a given date. date can be provided as a datetime, string, decimal year, modified julian day, or (doy,year) data_type can be ‘data’ or ‘data_xyz’

parameters:

date: datetime, string, float, tuple: the date to look for
data_type: string: ‘data’ or ‘data_xyz’ depending on the type of data to return

info(info=2): Print various informations about the time series

insert_gts_data(gts, in_place=False, verbose=False)

Insert data (and/or .data_xyz) of another Gts into the current Gts.

Parameters:

gts (Gts) – Time series to be inserted.
in_place (bool, optional) – If True, modify current Gts; if False, return a new Gts.
verbose (bool, optional) – Verbose mode.

Returns:

New Gts or the modified Gts if in_place.

Return type:

insert_ts(ts, rounding='day', data='xyz', overlap=True)

Insert another Gts into this one (merge or overwrite over its period).

Parameters:

ts (Gts) – Gts to be inserted.
rounding (str, optional) – Date rounding for deciding whether an entry is replaced: ‘second’, ‘minute’, ‘hour’, ‘day’.
data (str or None, optional) – Attribute to update: ‘xyz’ for .data_xyz or None for .data.
overlap (bool, optional) – If True, update only on matching dates; if False, ts overwrites current Gts over ts period.

Returns:

New Gts. .data or .data_xyz is set to None according to the data argument.

Return type:

interpolate(date='day', kind='linear', gap=10, in_place=False, verbose=False)

Interpolate the time series at regular or given dates.

Parameters:

date (str or ndarray, optional) – ‘day’ performs daily interpolation; alternatively a 1D numpy array with datetime or decimal year.
kind (str, optional) – scipy.interpolate.interp1d kind argument (e.g. ‘linear’).
gap (int or float, optional) – Maximum gap in days for interpolation; series is split at larger gaps.
in_place (bool, optional) – If True, modify in place; otherwise return a new Gts.
verbose (bool, optional) – Verbose mode.

Returns:

Interpolated time series (new instance unless in_place).

Return type:

ivel()

Compute instantaneous velocity (time series should be filtered first).

Returns:: New Gts instance with velocity; units m/yr.
Return type:: Gts

Notes

Output dates are shifted by 1 day so the first ivel date is the second date of the input time series, and so on.

l1trend(alpha='auto', ndays_mf=1, criterion='BIC', lcomponent='ENU', algo='golden', pre_process=4500, bounds=[-2, 1], tol=0.01, refine=True, min_samples_per_segment=5, threshold_bias_res_detect=5, threshold_vel=3, min_sample_segment_refinement=1, norm='L2', verbose=False, log_dir=None, component=None)

Performs l1 trend filter on GPS time series. For each component, optimal filter parameter is searched using specified criterion. Uses https://pypi.org/project/trendfilter/

Parameters:

alpha (float or str) – Either hyperparameter as float or ‘auto’ to find optimal filtering using criterion
ndays_mf (int) – Parameter for a median filter to be applied prior to l1-trend filtering
criterion (str) – Criterion to select automatic optimal l1-trend filtering. Options: ‘BIC’ (default), ‘AICc’, ‘Cp’
lcomponent (str) – List of components to be filtered. Other components will remain unchanged
algo (str) – Algorithm to find optimal alpha. Options: ‘golden’ (default) or ‘custom’
pre_process (float) – Pre-process time series to avoid convergence problems in l1trend. It corrects for unrealistically large velocity using the threshold provided (mm/yr). Correction is then added again after l1_trend. Default is 3600 (3.6 m/yr)
bounds (list) – Bounds for golden section search algorithm. Default is [-2, 1]
tol (float) – Tolerance for golden section search algorithm. Default is 0.01
refine (bool) – If True, will run refinement after l1 trend filtering
min_samples_per_segment (int) – Minimum number of samples per segment for refinement
threshold_bias_res_detect (float) – Threshold for bias detection in refinement
threshold_vel (float) – Threshold for velocity detection in refinement
min_sample_segment_refinement (int) – Minimum number of samples per segment in refinement
norm (str) – Norm to use for refinement. Options: ‘L2’ (default)
verbose (bool) – Enable verbose output for debugging
log_dir (str or None) – If provided, directory where the log file will be created. The log file will be named ‘SITE_l1trendi.log’ where SITE is the site code. If None, no file logging is performed.
component (str or None) – Components to process (‘N’, ‘E’, ‘U’, ‘NE’, ‘EN’, ‘NU’, ‘UN’, ‘EU’, ‘UE’, ‘NEU’, etc.). If None, uses lcomponent parameter. Used for component-specific statistics calculation.

Returns:

l1-trend filtered time series as new Gts instance

Return type:

Raises:

ValueError – If invalid parameters are provided
RuntimeError – If processing fails

l1trend_optimal_workflow(criterion='AICc', component='NEU', log_dir=None)

Optimal workflow for l1 trend analysis.

This workflow combines multiple L1-trend processing steps for optimal results: 1. Initial L1-trend filtering with AICc criterion 2. Breakpoint simplification to remove redundant breakpoints 3. Local refinement using original data for better fit 4. Final simplification to clean up the refined result

Parameters:

criterion (str, optional) – Criterion for L1-trend optimization (‘AICc’, ‘BIC’, ‘Cp’)
component (str, optional) – Components to process (‘N’, ‘E’, ‘U’, ‘NE’, ‘NU’, ‘EU’, ‘NEU’). Default is ‘NEU’
log_dir (str, optional) – Directory for Gts-specific logging. If provided, creates a log file named {self.code}.log

Returns:

Final optimized L1-trend time series

Return type:

l1trend_to_breakpoints(tol='auto', threshold=[1.0, 1.0, 5.0])

Convert a Gts resulting from a L1-trend-filtering to a dictionary of breakpoints. The breakpoints are computed by looking for significant changes in the slope of the time series.

Parameters:

tol (float or 'auto') – Tolerance for detecting breakpoints in mm/yr. If ‘auto’, finds the greatest tol such that the max difference between interpolated breakpoints and the time series is below a threshold for each component.
threshold (list of floats) – List of thresholds in mm for each component (‘N’, ‘E’, ‘U’) to determine the best tolerance if tol is ‘auto’. Default is [1., 1., 5.] for ‘N’, ‘E’, and ‘U’ respectively. If tol is a float, this parameter is ignored.

Returns:

bp – Dictionary with keys as component names (‘E’, ‘N’, ‘U’) and dates/values as bp[component][0], H_bp[component][1]

Return type:

l1trendi(alpha='auto', ndays_mf=1, criterion='BIC', lcomponent='ENU', algo='golden', pre_process=4500, bounds=[-2, 1], tol=0.01, refine=True, min_samples_per_segment=5, threshold_bias_res_detect=5, threshold_vel=3, min_sample_segment_refinement=1, norm='L2', verbose=False, log_dir=None, component=None)

Performs l1 trend filter on GPS time series. For each component, optimal filter parameter is searched using specified criterion. Uses https://pypi.org/project/trendfilter/

Parameters:

alpha (float or str) – Either hyperparameter as float or ‘auto’ to find optimal filtering using criterion
ndays_mf (int) – Parameter for a median filter to be applied prior to l1-trend filtering
criterion (str) – Criterion to select automatic optimal l1-trend filtering. Options: ‘BIC’ (default), ‘AICc’, ‘Cp’
lcomponent (str) – List of components to be filtered. Other components will remain unchanged
algo (str) – Algorithm to find optimal alpha. Options: ‘golden’ (default) or ‘custom’
pre_process (float) – Pre-process time series to avoid convergence problems in l1trend. It corrects for unrealistically large velocity using the threshold provided (mm/yr). Correction is then added again after l1_trend. Default is 3600 (3.6 m/yr)
bounds (list) – Bounds for golden section search algorithm. Default is [-2, 1]
tol (float) – Tolerance for golden section search algorithm. Default is 0.01
refine (bool) – If True, will run refinement after l1 trend filtering
min_samples_per_segment (int) – Minimum number of samples per segment for refinement
threshold_bias_res_detect (float) – Threshold for bias detection in refinement
threshold_vel (float) – Threshold for velocity detection in refinement
min_sample_segment_refinement (int) – Minimum number of samples per segment in refinement
norm (str) – Norm to use for refinement. Options: ‘L2’ (default)
verbose (bool) – Enable verbose output for debugging
log_dir (str or None) – If provided, directory where the log file will be created. The log file will be named ‘SITE_l1trendi.log’ where SITE is the site code. If None, no file logging is performed.
component (str or None) – Components to process (‘N’, ‘E’, ‘U’, ‘NE’, ‘EN’, ‘NU’, ‘UN’, ‘EU’, ‘UE’, ‘NEU’, etc.). If None, uses lcomponent parameter. Used for component-specific statistics calculation.

Returns:

l1-trend filtered time series as new Gts instance

Return type:

Raises:

ValueError – If invalid parameters are provided
RuntimeError – If processing fails

local_offset_robust(date, n, verbose=False, debug=False)

Estimate a local offset (no velocity) with a robust method.

Parameters:

date (float) – Offset date in decimal year.
n (int) – Number of dates before/after used in estimation.
verbose (bool, optional) – Verbose mode. Default is False.
debug (bool, optional) – Debug output. Default is False.

Returns:

1D array [date, north, east, up, s_north, s_east, s_up].

Return type:

numpy.ndarray

make_dynamic_apr(apr, time_step=30.0, pos_tol=0.03, dates=[], gap=20.0, verbose=False)

Create a dynamic apr file for GAMIT (coordinates at different times, no velocity).

Parameters:

apr (str) – Output apr file path (globk format).
time_step (float, optional) – Time step for writing dates in days. Default is 30.
pos_tol (float, optional) – Position tolerance (m); exceedance triggers a new date. Default is 0.03.
dates (list, optional) – Decimal-year dates where writing is forced. Default is [].
gap (float, optional) – Gap in days; if no data for longer than gap, observation is forced and tested against pos_tol. Default is 20.
verbose (bool, optional) – If True, print progress. Default is False.

Returns:

self.

Return type:

make_model(option='detrend', method='L2', loutlier=None): Estimate linear model parameters using least squares. input: data: Gts format option are: ‘detrend’/’detrend_annual’/’detrend_seasonal’ output: new Gts object: time series is now the residuals wrt to the model and its associated values (vel, annual, semi-annual etc).

median_filter(n, in_place=False, verbose=True)

Filter time series with scipy.signal.medfilt (median filter).

Parameters:

n (int) – Window size (must be odd).
in_place (bool, optional) – If True, replace self. Default is False.
verbose (bool, optional) – Verbose mode. Default is True.

Returns:

Filtered time series.

Return type:

Notes

Applied to .data; .data_xyz is set to None.

minimum_component(mask_period=[], p=1, fcut=None, Q=None, in_place=False, verbose=True)

Minimum component filtering for Gts (background signal in the presence of spikes).

Parameters:

mask_period (list, optional) – Period(s) to ignore for smoothing (list or list of lists).
p (int, optional) – Polynomial degree for the fit (default 1).
fcut (float, optional) – Cutoff frequency for the low-pass filter. Default f_nyq / sqrt(N).
Q (float, optional) – Strength of the low-pass filter; larger Q = steeper cutoff. Default 0.1 * fcut.
in_place (bool, optional) – If True, replace the current time series; otherwise return a new Gts.
verbose (bool, optional) – Verbose mode.

Returns:

Filtered time series.

Return type:

Notes

Follows “Practical Statistics for Astronomers” (Wall & Jenkins) and Wall, J, A&A 122:371, 1997.

mmodel(): Generates a modeled time series from the parameters read in self

n_obs(date)

Return the number of available observations for the given date or period.

Parameters:: date (str or list) – Date or period. Formats: decimal year or pandas-style (e.g. ‘2018-01-01’, [‘2018-01-01’,’2018-02-01’], or [2019., 2019.5] for decimal year).
Returns:: Number of observations in the current time series for the given date/period.
Return type:: int

Examples

A specific day: ts.RIOP.n_obs(‘2018-01-01’) A period: ts.RIOP.n_obs([‘2018-01-01’,’2018-02-01’]) A period in decimal year: ts.RIOP.n_obs([2019., 2019.5])

neu2xyz(corr=False, verbose=False)

Populate .data_xyz from .data (requires X0,Y0,Z0 or lon, lat, h).

Parameters:

corr (bool, optional) – If True, compute standard deviations and correlations. Default is False.
verbose (bool, optional) – Verbose mode. Default is False.

Returns:

self.

Return type:

np_datetime_2_eq_time(leap_sec=0.0, eq_time=0.0)

Convert dict of datetime to seconds relative to eq_time.

If input is in GPS time, leap_sec corrects for GPS-UTC (e.g. 17 s on 2016-02-13).

Parameters:

data (dict) – Mapping index -> datetime.datetime.
leap_sec (float, optional) – GPS time minus UTC in seconds (e.g. 17).
eq_time (datetime, optional) – Earthquake time in UTC (reference for seconds).

Returns:

Seconds relative to eq_time (and leap_sec correction if applied).

Return type:

ndarray

np_yyyy_mm_dd_hh_mm_ss_2_datetime()

Convert numpy array [year, month, mday, hour, minute, sec] to datetime objects.

Parameters:: data (ndarray) – Array of shape (n, 6).
Returns:: Mapping index -> datetime.datetime.
Return type:: dict

np_yyyy_mm_dd_hh_mm_ss_2_decyear()

Convert numpy array [year, month, mday, hour, minute, sec] to decimal year.

Parameters:: data (ndarray) – Array of shape (n, 6) with year, month, day, hour, minute, sec.
Returns:: 1D array of decimal years.
Return type:: ndarray

plot(title=None, loffset=True, loutliers=True, verbose=False, date=[], yaxis=None, min_yaxis=None, yupaxis=None, xticks_minor_locator=1, lcomponent=['N', 'E', 'U'], error_scale=1.0, lperiod=[[]], lvline=[], save_dir_plots='.', save=None, show=True, unit='mm', date_unit='cal', date_ref=0.0, center=True, superimposed=None, lcolor=['r', 'g', 'c', 'm', 'y', 'k', 'b'], label=None, legend=False, set_zero_at_date=None, grid=True, plot_size=None, info=[], xlabel_fmt=None, **kwargs)

Plot North-East-Up time series with offsets/outliers using Matplotlib.

Coordinates in meters; default display unit is mm (use unit=’m’ for meters).

Parameters:

title (str, optional) – Plot title (added to site name). Default is None.
loffset (bool, optional) – Draw vertical lines at offset dates. Default is True.
loutliers (bool, optional) – Plot outliers. Default is True.
verbose (bool, optional) – Verbose mode. Default is False.
date (list, optional) – [sdate, edate] for plot range (decimal years or days per date_unit). Default is [].
yaxis (list, optional) – [min_y, max_y] for y-axis; auto if not provided.
min_yaxis (list, optional) – Alternative y-axis limits.
yupaxis (list, optional) – Same as yaxis for Up component only.
xticks_minor_locator (float or str, optional) – Minor tick spacing; float for ‘decyear’/’days’, str e.g. ‘%Y’ for ‘cal’. Default is 1.
lcomponent (list, optional) – Components to plot. Default is [‘N’,’E’,’U’].
error_scale (float, optional) – Error bar scale (0 = no bars). Default is 1.0.
lperiod (list or dict, optional) – Background periods (light salmon); or dict per component. Default is [[]].
lvline (list, optional) – Dates for vertical lines (green). Default is [].
save_dir_plots (str, optional) – Directory for saving. Default is ‘.’.
save (bool or str, optional) – True = auto name, or filename. Default is None.
show (bool, optional) – If True, show plot. Default is True.
unit (str, optional) – ‘m’, ‘cm’, ‘mm’. Default is ‘mm’.
date_unit (str, optional) – ‘decyear’, ‘cal’, or ‘days’. Default is ‘cal’.
date_ref (float, optional) – Reference date. Default is 0.0.
center (bool, optional) – Center y-axis on mean. Default is True.
superimposed (Gts, optional) – Gts to overlay. Default is None.
lcolor (list, optional) – Colors for superimposed. Default is [‘r’,’g’,’c’,’m’,’y’,’k’,’b’].
label (str, optional) – Legend label for superimposed. Default is None.
legend (bool, optional) – Show legend. Default is False.
set_zero_at_date (float or list, optional) – Date (or [date, off_n, off_e, off_u]) to align master and superimposed. Default is None.
grid (bool, optional) – Draw grid. Default is True.
plot_size (tuple, optional) – Figure size. Default is best guess.
info (list, optional) – Titles for subplots. Default is [].
xlabel_fmt (str, optional) – Format for x-axis labels. Default is None.
**kwargs (dict, optional) – Passed to matplotlib.pyplot.errorbar (linewidth, marker, markersize, etc.).

Return type:

None

classmethod read(tsfile, fmt=None, verbose=False)[source]

Read a time series from file (format auto-detected or set by fmt).

Parameters:

tsfile (str) – Path to the time series file.
fmt (str, optional) – Format name (e.g. ‘pride’, ‘pride_pos’, ‘mb_file’). If None, formats are tried.
verbose (bool, optional) – If True, print progress. Default is False.

Returns:

Loaded time series or None on failure.

Return type:

Gts or None

read_cats_file(idir='.', ifile=None, gmt=True, verbose=False)

Read a CATS file and load the time series into this Gts.

Parameters:

idir (str, optional) – Directory for CATS files.
ifile (str, optional) – Path to CATS file. If None, uses idir/cats_<code>.dat.
gmt (bool or str, optional) – If True, read lon/lat from ../maps_en_velo.gmt; if str, path to GMT file.
verbose (bool, optional) – Verbose mode.

Returns:

self (data loaded from file).

Return type:

read_eq_rename(eq_rename, in_place=False, verbose=False)

Read current site info from a globk eq_rename file.

Populates outliers and offsets_dates; excluded periods from the file are added to outliers.

Parameters:

eq_rename (str) – Path to eq_rename file (globk format).
in_place (bool, optional) – If True, modify self; if False, return a new Gts. Default is False.
verbose (bool, optional) – If True, print progress. Default is False.

Returns:

self (if in_place) or new Gts instance.

Return type:

read_kenv(ifile, date_type='jd'): Read kenv file (magnet) format for time series

read_lon_lat(gmt_file, verbose=False)

Read a GMT psvelo file and set Gts.lon and Gts.lat.

Parameters:

gmt_file (str) – Path to GMT psvelo file.
verbose (bool, optional) – If True, print progress. Default is False.

Returns:

self (lon, lat populated).

Return type:

read_mb_file(idir='.', ifile=None, gmt=True, verbose=False): Read GAMIT/GLOBK mb_files in a directory and actually loads the time series

read_offset_dates(offset_file)

Read an offset file and set offsets_dates (pyacs format).

File format: code and dates; dates parsed by pyacs.guess_date.

Parameters:: offset_file (str) – Path to offset file.
Returns:: self (offsets_dates populated).
Return type:: Gts

read_pos(tsdir='.', tsfile=None, xyz=True, verbose=False)

Read GAMIT/GLOBK PBO pos file and load the time series into this Gts.

Parameters:

tsdir (str, optional) – Directory containing pos file(s). Default is ‘.’.
tsfile (str, optional) – Pos file to load. If None, a file CODE*.pos is sought. Default is None.
xyz (bool, optional) – If True, read XYZ and sx,sy,sz, corr columns. Default is True.
verbose (bool, optional) – If True, print progress. Default is False.

Returns:

self (data, code, X0, Y0, Z0, t0 populated).

Return type:

Notes

A pos file contains (almost) all needed info. If tsfile is None, read_pos looks for a file named CODE*.pos.

read_pride(tsdir='.', tsfile=None, xyz=True, verbose=False)

Read PRIDE-PPPAR kinematic result file.

Parameters:

tsdir (str, optional) – Directory of pride-pppar kinematic files.
tsfile (str, optional) – Pride-pppar kinematic file to load. If None, looks for kin_*code.
xyz (bool, optional) – If True, keep data in XYZ; otherwise convert to NEU.
verbose (bool, optional) – Verbose mode.

Notes

If tsfile is None, read_pride looks for a file named kin_*code.

read_pride_pos(tsdir='.', tsfile=None, verbose=False)

Read PRIDE-PPPAR static result file(s).

Parameters:

tsdir (str, optional) – Directory containing pride-pppar pos static files. Default ‘.’.
tsfile (str, optional) – Specific file to load. If None, searches for pos_*code in tsdir.
verbose (bool, optional) – If True, enable verbose output.

Notes

If tsfile is None, looks for files named pos_*<code> in tsdir.

read_series(tsdir='.', tsfile=None, xyz=True, verbose=False)

Read GipsyX series time series file.

Parameters:

tsdir (str, optional) – Directory containing the series file.
tsfile (str, optional) – Series file to load. If None, looks for CODE*.series.
xyz (bool, optional) – If True, keep/store XYZ; otherwise NEU.
verbose (bool, optional) – Verbose mode.

Notes

This format does not include the absolute position; relative only.

read_sol(tsdir='.', ifile=None, verbose=False)

Read SGC sol file and load the time series into this Gts instance.

Parameters:

tsdir (str, optional) – Directory of sol/pos files. Default is ‘.’.
ifile (str, optional) – Path to the sol file to load. If None, a file CODE*.pos may be used.
verbose (bool, optional) – If True, print progress. Default is False.

Returns:

self (data, code, X0, Y0, Z0, t0 are populated from file).

Return type:

Notes

A pos/sol file contains (almost) all needed info. If ifile is None, read_pos looks for a file named CODE*.pos.

read_tdp(idir='.', ifile=None, gmt=True, verbose=False): Read tdp (Gipsy kinematics provided by Cedric Twardzik 17/04/2018) format for time series

read_track_NEU(tsdir='.', tsfile=None, leap_sec=0.0, eq_time=None, verbose=False): Read a GAMIT/GLOBK Track output file generated with the option out_type NEU in this case dates are seconds by default the seconds are with respect to the first epoch of measurements If option leap_sec is provided with a value > 0.0, then GPS time is corrected for the difference between GPTS time and UTC If eq_time is provided, it is assumed to be UTC. Expected format is YYYY:MM:MD:HH:MM:SS.S

realistic_sigma(option='tsfit', in_place=False, verbose=False)

Calculate realistic velocity uncertainties (GLOBK tsfit or CATS).

Parameters:

option (str, optional) – ‘tsfit’: GLOBK T. Herring realistic sigma; ‘cats_pl’: CATS with noise type (–model=pl:); ‘cats_seasonal_pl’: CATS with seasonal (–model=pl: –sinusoid=1y1); ‘cats_flicker’: CATS flicker (–model=pl:k-1); ‘cats_seasonal_flicker’: CATS seasonal + flicker.
in_place (bool, optional) – If True, update Gts in place.
verbose (bool, optional) – Verbose mode.

Returns:

self (possibly with updated velocity uncertainties) or result of CATS/tsfit.

Return type:

Gts or None

refine_l1trend(rawts, lcomponent='ENU', min_samples_per_segment=10, threshold_bias_res_detect=10, threshold_vel=8, min_sample_segment_refinement=1, norm='L2', output='ts')

Refine results from l1trend by optimizing the date of breakpoints for suspected periods. The algorithm first identifies periods where the l1trend model might be improved. For every candidate period, the algorithm checks different metrics to decide if the model should be refined. Finally, for the selected periods, it optimizes the date of breakpoints.

Parameters:

rawts (pyacs.gts.Gts.Gts) – Raw time series, originally processed by l1trend
lcomponent (str) – Component to refine (default: ‘ENU’)
min_samples_per_segment (int) – Minimum number of samples for a segment to be considered as a candidate for improvement (default: 10)
threshold_bias_res_detect (float) – Threshold to detect bias residuals (default: 10)
threshold_vel (float) – Threshold to detect high velocities (default: 8)
min_sample_segment_refinement (int) – Minimum number of samples for a segment in the refined model (default: 1)
norm (str) – Norm to be minimized for improvement (default: ‘L2’)
output (str) – Output type ‘ts’ for the refined Gts, ‘info’ for the periods suspected, ‘both’ for both (default: ‘ts’)

Returns:

Refined Gts object or information about suspected periods

Return type:

pyacs.gts.Gts.Gts or tuple

remove_outliers(periods=None, in_place=False, log_dir=None)

Remove outliers listed in self.outliers.

Parameters:

periods (list, optional) – If provided, restrict to these periods. Default is None.
in_place (bool, optional) – If True, modify self; otherwise return a new Gts. Default is False.
log_dir (str, optional) – Directory for Gts-specific log file {self.code}.log. Default is None.

Returns:

New Gts without outliers, or self if in_place=True.

Return type:

Gts or None

remove_pole(pole, pole_type='euler', verbose=True): Remove velocity predicted by an Euler pole or a rotation rate vector from a time series. pole is a 1D array with 3 values. Requires self.lon & self.lat attributes to have been filled before. Returns a new Gts instance.

remove_velocity(vel_neu, in_place=False): remove velocity from a time series vel_neu is a 1D array of any arbitrary length, but with the velocities (NEU) to be removed in the first 3 columns if in_place = True then replace the current time series

reorder(verbose=False)

Reorder .data and/or .data_xyz by increasing dates (always in place).

Parameters:: verbose (bool, optional) – Verbose mode.
Returns:: self (modified in place).
Return type:: Gts

rotate(angle, in_place=False)

Rotate the horizontal axes (E, N) by an angle.

Parameters:

angle (float) – Angle in decimal degrees, clockwise.
in_place (bool, optional) – If True, replace the current time series; otherwise return a new Gts.

Returns:

Rotated time series (self if in_place, else new Gts). .data_xyz set to None.

Return type:

save_apr(apr, epoch=None, verbose=False, excluded_periods=None)

Save Gts analysis results to a globk-format apr file.

Parameters:

apr (str) – Output apr file path (globk format).
epoch (float, optional) – Epoch in decimal year for coordinates. Default is None (use period mid).
verbose (bool, optional) – If True, print progress. Default is False.
excluded_periods (list, optional) – Periods to exclude. Default is None.

Returns:

self.

Return type:

Notes

Following globk convention, site names are XXXX_1PS, XXXX_2PS, etc. between offset dates.

save_eq_rename(eq_rename, verbose=False, excluded_periods=None)

Save Gts analysis results to a globk-format eq_rename file.

Parameters:

eq_rename (str) – Output eq_rename file path (globk format).
verbose (bool, optional) – If True, print progress. Default is False.
excluded_periods (list, optional) – Periods to exclude. Default is None.

Returns:

self.

Return type:

save_offsets(ofile, verbose=True, comment='', up=False, info=False)

Append offset values to a text file (GMT psvelo format).

Parameters:

ofile (str) – Output offset file path.
verbose (bool, optional) – If True, print progress. Default is True.
comment (str, optional) – Comment; ‘#’ is prepended if not present. Default is ‘’.
up (bool, optional) – If True, Ve/SVe/SVen set to 0 and Vu as 4th/6th fields. Default is False.
info (bool, optional) – If True, write extra info. Default is False.

Returns:

self.

Return type:

save_velocity(gmt_file, verbose=True, comment=None, up=False)

Append velocity estimates (with uncertainties) to a GMT psvelo file.

Parameters:

gmt_file (str) – Output GMT psvelo file (appended if it exists).
verbose (bool, optional) – If True, print progress. Default is True.
comment (str, optional) – Comment line; ‘#’ is prepended if not present.
up (bool, optional) – If True, Ve/SVe/SVen set to 0 and Vu written as 4th/6th fields. Default is False.

Returns:

self.

Return type:

savitzky_golay(in_place=False, verbose=True, window_length=15, polyorder=3, deriv=0, delta=1.0, mode='interp', cval=0.0)

Return a filtered time series using scipy.signal.savgol_filter.

Parameters:

in_place (bool, optional) – If True, replace the current time series; otherwise return a new Gts.
verbose (bool, optional) – Verbose mode.
window_length (int, optional) – Length of the filter window (must be odd).
polyorder (int, optional) – Order of the polynomial.
deriv (int, optional) – Order of derivative (0 = smoothing).
delta (float, optional) – Sample spacing.
mode (str, optional) – Extension mode (‘interp’, ‘mirror’, etc.).
cval (float, optional) – Value for mode=’constant’.

Returns:

Filtered time series.

Return type:

See also

scipy.signal.savgol_filter

set_zero_at_date(date, offset=None)

Translate the time series so that values are zero at a given date.

If the provided date does not exist in the series, the next available date is used.

Parameters:

date (float) – Date in decimal year.
offset (float or list or ndarray, optional) – Offset in mm to add. Float (same for N,E,U) or 3-element list/array for N,E,U.

Returns:

New Gts with translation applied; .data_xyz set to None.

Return type:

sigma_cats(in_place=False, verbose=False, k='k-1', seasonal='')

Run CATS to obtain realistic velocity uncertainties.

Parameters:

in_place (bool, optional) – If True, update Gts in place.
verbose (bool, optional) – Verbose mode.
k (str, optional) – CATS noise model (e.g. ‘k-1’ for flicker).
seasonal (str, optional) – Seasonal option for CATS.

Returns:

self or result from CATS run.

Return type:

Gts or None

sigma_vel_tsfit(in_place=False, verbose=False): runs tsfit for getting realistic sigma

simplify_l1trend(tolerance=0.5, components='ENU')

Remove unnecessary breakpoints from an L1-trend filtered time series.

This function iteratively removes breakpoints and tests if the simplified model still fits the original time series within a specified tolerance.

Parameters:

tolerance (float) – Maximum allowed difference (in mm) between original and simplified model. Default is 0.5 mm.
components (str) – Components to process. Default is ‘ENU’.

Returns:

Simplified Gts object with unnecessary breakpoints removed.

Return type:

pyacs.gts.Gts.Gts

simplify_l1trend_with_fisher_test(rawts, components='ENU', alpha=0.05)

Simplify an L1-trend filtered time series by comparing to the unfiltered time series using Fisher-Snedecor tests.

This function iteratively removes breakpoints and tests their usefulness by comparing the fit quality with and without each breakpoint using Fisher-Snedecor tests.

Parameters:

rawts (pyacs.gts.Gts.Gts) – The unfiltered (raw) time series to compare against.
components (str) – Components to process. Default is ‘ENU’.
alpha (float) – Significance level for the Fisher-Snedecor test. Default is 0.05.

Returns:

Simplified Gts object with unnecessary breakpoints removed.

Return type:

pyacs.gts.Gts.Gts

smooth(window_len=11, window='hanning', in_place=False, verbose=False, component='NEU'): smooth a time series

spline(smoothing=1, degree=5, date=None)

Model the time series with a smoothing spline (UnivariateSpline).

Parameters:

smoothing (float, optional) – Positive smoothing factor; number of knots increased until sum((w*(y-spl(x)))**2) <= s (s = smoothing * 1e-3).
degree (int, optional) – Degree of the spline (<= 5). Default 5.
date (ndarray or str or None, optional) – Interpolation dates: 1D array (decimal year), ‘day’ for daily, or None for data dates only.

Returns:

New Gts instance with spline-smoothed NEU.

Return type:

split(split_dates, verbose=True)

Split time series at given dates.

Dates <= split date define the left segment. Dates are converted to integer seconds via at.decyear2seconds.

Parameters:

split_dates (list or numpy.ndarray) – Decimal-year dates for split points.
verbose (bool, optional) – Verbose mode. Default is True.

Returns:

Gts objects from the split.

Return type:

list of Gts

split_gap(gap=10, verbose=False)

Split the time series at gaps larger than a given number of days.

Parameters:

gap (int or float, optional) – Gap in days above which the series is split.
verbose (bool, optional) – Verbose mode.

Returns:

List of Gts split from the original.

Return type:

list of Gts

substract_ts(ts, tol=0.05, verbose=True)

Subtract the provided time series from the current one (at common dates).

Parameters:

ts (Gts) – Time series to subtract (Gts instance).
tol (float, optional) – Date tolerance in days to consider two dates identical (default 0.05 day).
verbose (bool, optional) – Verbose mode.

Returns:

New Gts with code self.code+’_’+ts.code; .data_xyz set to None.

Return type:

substract_ts_daily(ts, verbose=True)

Subtract the given time series from the current one (sample-by-sample).

Parameters:

ts (Gts) – Time series to subtract (Gts instance).
verbose (bool, optional) – If True, print progress. Default is True.

Returns:

New Gts instance (current minus ts).

Return type:

Notes

Assumes daily time series; dates are matched.

sum_l1_trend()

Summarize L1-trend filtered time series with comprehensive statistics.

This function analyzes a previously L1-trend filtered time series and provides detailed statistics including: - Time series period - Number of breakpoints per component - Velocity statistics (mean, median, maximum) - Median duration between successive breakpoints

Parameters:: self (Gts object) – The L1-trend filtered time series to analyze.
Returns:: The original time series object (unchanged)
Return type:: Gts

Notes

This method is designed to work on time series that have already been processed with L1-trend filtering. It extracts breakpoints and computes velocity statistics for each component (N, E, U).

suspect_offsets(threshold=3, verbose=True, lcomponent='NE', n_max_offsets=10, in_place=False): Try to find offsets in a time series from day-to-day differences.

suspect_offsets_mf(threshold=3, verbose=True, lcomponent='NE', n_max_offsets=5, in_place=False, debug=False): Try to find offsets in a time series using a median filter.

test_offset_significance(date, conf_level=95, lcomponent='NE', verbose=True, debug=False, mode='local')

Test whether an offset is statistically significant.

Parameters:

date (float) – Offset date in decimal year.
conf_level (float, optional) – Confidence level in percent. Default is 95.
lcomponent (str, optional) – Components to test (‘N’,’E’,’U’). Default is ‘NE’.
verbose (bool, optional) – Verbose mode. Default is True.
debug (bool, optional) – Debug output. Default is False.
mode (str, optional) – ‘local’, ‘detrend’, or ‘detrend_seasonal’. Default is ‘local’.

Returns:

True if significant, else False.

Return type:

bool

test_offsets(verbose=False, debug=True, window_length=None): Test offsets: delete small, F-ratio test, re-check small offsets.

to_pandas_df(data_xyz=False, uncertainty=False, round=False)

Convert a pyacs Gts to a pandas DataFrame.

Parameters:

data_xyz (bool, optional) – If True, use .data_xyz (X,Y,Z); otherwise use .data (N,E,U).
uncertainty (bool, optional) – If True, include uncertainty columns.
round (bool, optional) – If True, round index (unused).

Returns:

Index is datetime; columns are position and optionally uncertainties.

Return type:

pandas.DataFrame

Notes

Uncertainties are included only when uncertainty=True.

to_pytrf()

Convert pyacs Gts to pytrf ts object.

Uses the pytrf library from https://github.com/prebischung/pytrf.

Returns:: Time series in pytrf format (East, North, Up; MJD; covariances).
Return type:: pytrf.ts.ts

trajectory(model_type, offset_dates=[], eq_dates=[], H_fix={}, H_constraints={}, H_bounds={}, component='NEU', verbose=False)

Estimate parameters of a (non-linear) trajectory model for the time series.

Model: y(t) = trend + annual + semi-annual + offsets + post-seismic (psd_log).

Parameters:

model_type (str) – String of key-word parameters to estimate: ‘trend’, ‘annual’, ‘semi-annual’, ‘seasonal’, ‘offset’, ‘psd_log’. E.g. ‘trend-seasonal-offset-psd_log’ for full model.
offset_dates (list, optional) – List of offset dates in decimal year.
eq_dates (list, optional) – List of earthquake dates for post-seismic (psd_log) estimation.
H_fix (dict, optional) – Parameters to fix, e.g. {‘psd_log_offset_00’: [10., 15., 0.], ‘psd_log_tau_00’: [100., 100., 100.]}.
H_constraints (dict, optional) – Parameters to constrain (center, sigma), e.g. {‘psd_log_tau_01’: [[500., 50], [500., 50], [500., 50]]}.
H_bounds (dict, optional) – Bounds, e.g. {‘psd_log_tau_02’: [[2*365., 3*365.], …]}.
component (str, optional) – Components to estimate (‘NEU’ or subset).
verbose (bool, optional) – Verbose mode.

Returns:

(results_dict, model_Gts, residual_Gts, daily_predictions_Gts). Unlike most pyacs.gts functions, trajectory returns these 4 elements.

Return type:

tuple

vondrak(fc, in_place=False, verbose=True, component='NEU')

Return a filtered Gts using a Vondrak filter.

Parameters:

fc (float) – Cutoff frequency in cycles per year.
in_place (bool, optional) – If True, replace the current time series; otherwise return a new Gts.
verbose (bool, optional) – Verbose mode.
component (str, optional) – Components to filter (‘NEU’ or subset).

Returns:

Filtered time series.

Return type:

wiener(in_place=False, verbose=True, my_size=15, noise=None)

Return a filtered time series using scipy.signal.wiener.

Parameters:

in_place (bool, optional) – If True, replace the current time series; otherwise return a new Gts.
verbose (bool, optional) – Verbose mode.
my_size (int, optional) – Size of the Wiener filter window.
noise (float, optional) – Noise estimate for the filter.

Returns:

Filtered time series.

Return type:

See also

scipy.signal.wiener

write_cats(idir='./cats', offsets_dates=None, add_key='')

Write a file for CATS analysis.

If offsets_dates is not None, offset date information is added at the beginning.

Parameters:

idir (str, optional) – Directory to save the file (default ‘./cats’).
offsets_dates (list, optional) – List of offset dates to write in the header.
add_key (str, optional) – Additional key in the file name (e.g. cats_<code>_<add_key>.dat).

Return type:

None

write_mb_file(idir, add_key='')

write_pos(idir='./pos', add_key='', force=None, verbose=False)

Write time series in GAMIT/GLOBK PBO pos format.

Parameters:

idir (str, optional) – Output directory. Default is ‘./pos’.
add_key (str, optional) – If non-blank, output file is CODE_add_key.pos; otherwise CODE.pos.
force (str, optional) – ‘data’ or ‘data_xyz’ to force source; None for default behavior.
verbose (bool, optional) – If True, print progress. Default is False.

Returns:

self.

Return type:

Notes

Default (force=None): if both data and data_xyz exist, both are written; if only data, uses X0,Y0,Z0 to write data_xyz; if only data_xyz, recreates data and writes.

wrms()

Return the weighted RMS of the time series (NEU).

Returns:: np.array([wrms_n, wrms_e, wrms_up]).
Return type:: ndarray

xyz2neu(corr=False, ref_xyz=None, verbose=False)

Populate .data (NEU) from .data_xyz (XYZ); lon, lat and h are also set.

Parameters:

corr (bool, optional) – If True, standard deviations and correlations are propagated to NEU.
ref_xyz (list or ndarray, optional) – [X, Y, Z] for the origin of the local NEU frame. If not provided, first position is used.
verbose (bool, optional) – Verbose mode.

Returns:

self (modified in place).

Return type:

Notes

This method always modifies in place.

pyacs.gts.Gts.get_index_from_dates(dates, data, tol=0.25)[source]

Return indices in data whose first column matches given dates within tolerance.

Parameters:

dates (list of float) – Target dates in decimal year.
data (numpy.ndarray) – 2D array with decimal dates in first column.
tol (float, optional) – Date tolerance in days to consider a match. Default is 0.25.

Returns:

Indices into data for each date (or [] if no match).

Return type:

list