pyacs.gts.lib.l1trend package

Submodules

• pyacs.gts.lib.l1trend.breakpoints module
  • l1trend_to_breakpoints()
• pyacs.gts.lib.l1trend.check_trend module
  • check_l1_trend()
  • compute_measure_bias()
  • make_node_ts()
  • measure_bias()
• pyacs.gts.lib.l1trend.cleaning module
  • clean_l1trend()
• pyacs.gts.lib.l1trend.l1trend2d module
  • find_optimal_lambda()
  • l1_trendfilter_2d()
  • plot_lambda_selection()
  • second_diff_matrix()
• pyacs.gts.lib.l1trend.l1trendi module
  • l1trendi()
• pyacs.gts.lib.l1trend.optimal_l1trend_workflow module
  • l1trend_optimal_workflow()
• pyacs.gts.lib.l1trend.optimal_refinement module
  • optimal_pwlf_refinement()
  • optimal_pwlf_refinement_fast()
  • preprocess_timeseries_for_optimization()
• pyacs.gts.lib.l1trend.optimization module
  • best_l1trend_custom()
  • best_l1trend_golden()
• pyacs.gts.lib.l1trend.outlier_flagging module
  • flag_outliers_using_l1trend()
• pyacs.gts.lib.l1trend.preprocessing module
  • pre_process_test()
  • pre_process_ts()
• pyacs.gts.lib.l1trend.refinement module
  • refine_l1trend()
• pyacs.gts.lib.l1trend.simplification module
  • simplify_l1trend()
  • simplify_l1trend_with_fisher_test()
• pyacs.gts.lib.l1trend.statistics module
  • get_stats_l1_model()
• pyacs.gts.lib.l1trend.sum_l1_trend module
  • sum_l1_trend()

Module contents

L1-trend analysis and refinement tools for pyacs.

This package contains functions for L1-trend filtering, refinement, and analysis of GPS time series data.

pyacs.gts.lib.l1trend.best_l1trend_custom(x, y, criterion_idx, logger=None, component_mask=None)

Find the optimal hyperparameter alpha in l1trend using a custom search algorithm.

Parameters:

• x (numpy.ndarray) – Input time array

• y (numpy.ndarray) – Input data array

• criterion_idx (int) – Index of the criterion to use (-1 for BIC, -2 for AICc, -3 for Cp)

• logger (logging.Logger, optional) – Logger instance for logging messages

Returns:

(optimal filtered data, history dictionary, optimal alpha)

Return type:

tuple

pyacs.gts.lib.l1trend.best_l1trend_golden(x, y, criterion_idx, bounds=[-2, 1], tol=0.01, logger=None, component_mask=None)

Find the optimal hyperparameter alpha in l1trend using golden section search algorithm.

Parameters:

• x (numpy.ndarray) – Input time array

• y (numpy.ndarray) – Input data array

• criterion_idx (int) – Index of the criterion to use (-1 for BIC, -2 for AICc, -3 for Cp)

• bounds (list) – Bounds for the search [lower, upper]

• tol (float) – Tolerance for convergence

• logger (logging.Logger, optional) – Logger instance for logging messages

Returns:

Optimally filtered data

Return type:

numpy.ndarray

pyacs.gts.lib.l1trend.check_l1_trend(ts, 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

pyacs.gts.lib.l1trend.clean_l1trend(self, 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

pyacs.gts.lib.l1trend.find_optimal_lambda(Y, lambda_range=None, method='cross_validation', cv_folds=5, solver='CLARABEL', verbose=False)

Find the optimal lambda value for L1-trend filtering using various criteria.

Parameters:

• Y (array_like) – Input data of shape (T, 2)

• lambda_range (array_like, optional) – Range of lambda values to test. If None, uses a logarithmic range.

• method (str, optional) – Method to use for optimization. Options: ‘cross_validation’, ‘aic’, ‘bic’, ‘gcv’

• cv_folds (int, optional) – Number of folds for cross-validation

• solver (str, optional) – Solver to use for optimization

• verbose (bool, optional) – Whether to print progress information

Returns:

• float – Optimal lambda value

• dict – Dictionary containing optimization results

pyacs.gts.lib.l1trend.flag_outliers_using_l1trend(self, 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:

Gts

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.

pyacs.gts.lib.l1trend.get_stats_l1_model(x, y, fy, alpha, component_mask=None)

Calculate statistics for L1-trend model evaluation.

Parameters:

• x (numpy.ndarray) – Input time as 1D array

• y (numpy.ndarray) – Input raw data as 1D array

• fy (numpy.ndarray) – Input l1trend filtered data as 1D array

• alpha (float) – Hyperparameter used for l1 filtering

• component_mask (numpy.ndarray, optional) – Boolean mask indicating which data points should be considered for statistics. If None, all points are considered.

Returns:

(alpha, cchi2, sigma, cp, cchi2_probability, Cp, AICc, BIC)

Return type:

tuple

pyacs.gts.lib.l1trend.l1_trendfilter_2d(Y, lam=1.0, solver='CLARABEL', verbose=False, scs_kwargs=None, clarabel_kwargs=None, warm_start=None)

pyacs.gts.lib.l1trend.l1trend_optimal_workflow(self, 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:

Gts

pyacs.gts.lib.l1trend.l1trend_to_breakpoints(self, 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:

record

pyacs.gts.lib.l1trend.l1trendi(self, 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:

Gts

Raises:

• ValueError – If invalid parameters are provided

• RuntimeError – If processing fails

pyacs.gts.lib.l1trend.optimal_pwlf_refinement(t, y, y0, yn, nbp, weights=None)

Find the optimal piecewise linear function with nbp breakpoints fitting y.

The function finds the optimal piecewise linear function (opwlf) fitting y with nbp breakpoints according to the weighted L2 norm. The obtained optimal piecewise linear function verifies opwlf[0]=y0 and opwlf[-1]=yn, that is end and starting point of y are kept fixed.

The problem is applied to small dimension (n<=6), so that a systematic search can be performed and does not require complex optimization approaches.

For large time series (>100 points), preprocessing is applied to decimate the data while preserving essential structure.

Parameters:

• t (numpy.ndarray) – Time array expressed in integer seconds

• y (numpy.ndarray) – Values of the function to fit

• y0 (float) – Value of the function at t[0] (fixed)

• yn (float) – Value of the function at t[-1] (fixed)

• nbp (int) – Number of breakpoints (excluding start and end points)

• weights (numpy.ndarray, optional) – Weights for each data point. If None, all weights are set to 1.0.

Returns:

(optimal_breakpoints, optimal_values, min_error) - optimal_breakpoints: array of breakpoint times - optimal_values: array of function values at breakpoints - min_error: minimum weighted L2 error achieved

Return type:

tuple

pyacs.gts.lib.l1trend.optimal_pwlf_refinement_fast(t, y, y0, yn, weights=None)

Fast version of optimal_pwlf_refinement using l1trendi with criterion AICc.

Parameters:

• t (numpy.ndarray) – Time array expressed in integer seconds

• y (numpy.ndarray) – Values of the function to fit

• y0 (float) – Value of the function at t[0] (fixed)

• yn (float) – Value of the function at t[-1] (fixed)

• weights (numpy.ndarray, optional) – Weights for each data point. If None, all weights are set to 1.0.

Returns:

(optimal_breakpoints, optimal_values, min_error)

Return type:

tuple

pyacs.gts.lib.l1trend.plot_lambda_selection(results, title=None)

Plot lambda selection results.

Parameters:

• results (dict) – Results from find_optimal_lambda

• title (str, optional) – Title for the plot

pyacs.gts.lib.l1trend.pre_process_test(gts, component, threshold=5000.0, logger=None)

Pre-process. Test whether a large offset is present in the time series. Find the largest offset exceeding threshold found in the time series. Returns the index in the time series immediately after the largest offset or a void list if no offset is found

Parameters:

• gts (Gts) – Input time series

• component (str) – Component to test (‘E’, ‘N’, or ‘U’)

• threshold (float) – Threshold for offset detection

• logger (logging.Logger, optional) – Logger instance for logging messages

Returns:

List of indices where offsets are found

Return type:

list

pyacs.gts.lib.l1trend.pre_process_ts(self, threshold=5000.0)

Pre-process a time series to avoid convergence problems in l1trend. It performs change point detection analysis, correct for unrealistically large velocity and return a cleaned time series

pyacs.gts.lib.l1trend.preprocess_timeseries_for_optimization(t, y, target_length=50)

Preprocess time series for optimization by applying median filter and decimation.

This function reduces the number of points while preserving the essential structure of the time series, making optimization faster.

Parameters:

• t (numpy.ndarray) – Time array

• y (numpy.ndarray) – Values array

• target_length (int) – Target length for the decimated series (default: 50)

Returns:

(t_decimated, y_decimated, decimation_factor) - t_decimated: decimated time array - y_decimated: decimated values array - decimation_factor: factor used for decimation

Return type:

tuple

pyacs.gts.lib.l1trend.refine_l1trend(self, 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

pyacs.gts.lib.l1trend.simplify_l1trend(self, 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

pyacs.gts.lib.l1trend.simplify_l1trend_with_fisher_test(self, 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

pyacs.gts.lib.l1trend.sum_l1_trend(self)

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).