pyacs.gts.lib package

Subpackages

Submodules

pyacs.gts.lib.errors module

Exception class for Gts

exception pyacs.gts.lib.errors.GtsError[source]

Bases: Exception

exception pyacs.gts.lib.errors.GtsInputDataNone(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputData_xyzNone(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputData_allNone(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputDataTypeError(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputDataBadDim(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputDataBadNcolumns(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputDataBadNrows(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputData_xyzTypeError(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputData_xyzBadDim(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputData_xyzBadNcolumns(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputData_xyzBadNrows(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputDataDiffShape(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsInputDataDiffDate(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsCDataError(method_name, lib, gts)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsReadFileError(method_name, lib, file_name)[source]

Bases: pyacs.gts.lib.errors.GtsError

exception pyacs.gts.lib.errors.GtsMethodError(method_name, lib, file_name)[source]

Bases: pyacs.gts.lib.errors.GtsError

pyacs.gts.lib.gts_dates module

This module gathers a few useful date functions, which can be used to read different time series formats. pyacs.pygts uses decimal year as time defaults. Seconds for High-Rate GPS solution can also be used.

pyacs.gts.lib.gts_dates.np_yyyy_mm_dd_hh_mm_ss_2_decyear(data)[source]

converts a numpy array including year month mday hour minute sec to decimal year returns a 1D array

pyacs.gts.lib.gts_dates.np_yyyy_mm_dd_hh_mm_ss_2_datetime(data)[source]

converts a numpy array including year month mday hour minute sec to an array of python datetime.datetime object returns a hash

pyacs.gts.lib.gts_dates.np_datetime_2_eq_time(data, leap_sec=0.0, eq_time=0.0)[source]

takes a hash of python datetime.datetime object and return a numpy array of seconds with respect to eq_time if the input array is in GPS time, providing leap_sec correct for the GPS_time - UTC delta

Parameters
  • leap_sec – number of seconds between GPS_time - UTC delta (leap_sec=17 that is GPS is ahead of UTC by 17 seconds on 13/02/2016)

  • eq_time – time of earthquake as a python datetime.datetime object (in UTC)

pyacs.gts.lib.gts_dates.decyear2days(self, ref_date='', in_place=False)[source]

Converts the dates of a time series from decimal years to days after a reference date ref_date is read by guess_date

pyacs.gts.lib.gts_estimators module

pyacs.gts.lib.gts_estimators.least_square(A, L, P=None)[source]

Least squares estimation for system equation AX + L = 0, P input: A: design matrix;L observation vector L; P: weight matrix for L (P defaut is the identity matrix) output: unknown matrix : X residuals matrix: V standard deviation sigma_0: std unknown parameters variance: s_X residuals variance: s_V model S:S = A*X

pyacs.gts.lib.message module

Print messages and log

pyacs.gts.lib.message.message(message, log_file=None)[source]
Parameters
  • message – message to be printed

  • log_file – file object where message will be written

Returns

pyacs.gts.lib.metadata module

pyacs.gts.lib.metadata.read_lon_lat(self, gmt_file, verbose=False)[source]

Reads a gmt psvelo file and populates Gts.lon & Gts.lat

Parameters
  • gmt_file – gmt psvelo file

  • verbose – verbose mode (boolean)

Returns

the current Gts instance

pyacs.gts.lib.metadata.save_velocity(self, gmt_file, verbose=True, comment=None, up=False)[source]

Appends velocity estimates (with uncertainties) to a gmt psvelo file

Parameters
  • gmt_file – output gmt psvelo file (will append if gmt_file already exists)

  • verbose – verbose mode (boolean)

  • comment – comment as a string. ‘# ‘ is pre-prended to comment if not provided

  • up – boolean. If True, then Ve, SVe and SVen are set to 0 and Vu and Vu are written as 4-th and 6-th fields

Returns

the current Gts instance

pyacs.gts.lib.metadata.save_offsets(self, ofile, verbose=True, comment='', up=False, info=False)[source]

Appends offsets values to a given text file (gmt psvelo format)

Parameters
  • ofile – output offset file

  • verbose – verbose mode (boolean)

  • comment – comment as a string. ‘# ‘ is pre-prended to comment if not provided

  • up – boolean. If True, then Ve, SVe and SVen are set to 0 and Vu and Vu are written as 4-th and 6-th fields

Returns

the current Gts instance

pyacs.gts.lib.metadata.read_eq_rename(self, eq_rename, in_place=False, verbose=False)[source]

Reads the information for the current site (code) from an eq_rename globk file.

Populates loutliers and offsets_dates Found excluded periods in the eq_rename file are added to loutliers

Parameters
  • eq_rename – eq_rename (globk format) file to be read

  • in_place – boolean. If True then the Gts instance is modified, if False the Gts instance is preserved and a new Gts instance is return

  • verbose – verbose mode (boolean)

Returns

Gts instance

pyacs.gts.lib.metadata.save_eq_rename(self, eq_rename, verbose=False, excluded_periods=None)[source]

save results of a Gts analysis in globk format eq_rename

Parameters
  • eq_rename – output eq_rename file (Golbk format)

  • verbose – verbose mode (boolean)

  • exluded_periods – periods to be excluded

Returns

Gts instance

pyacs.gts.lib.metadata.make_dynamic_apr(self, apr, time_step=30.0, pos_tol=0.03, dates=[], gap=20.0, verbose=False)[source]

Creates an apr file for GAMIT The created apr file has no velocity, but a series of coordinates at different time

Parameters
  • apr – apr file (Globk format)

  • time_step – time step for writing dates (default 30 days)

  • pos_tol – position tolerance. If exceeded, a new date will be written. (default 0.03 m)

  • dates – a list of dates in decimal years where writing will be forced

  • gap – gap in days. If there is no data during a duration greater than gap, then observation is forced to be included and tested against pos_tol

  • verbose – verbose mode (boolean)

Returns

Gts instance

pyacs.gts.lib.metadata.save_apr(self, apr, epoch=None, verbose=False, excluded_periods=None)[source]

save results of a Gts analysis in globk format apr file

Parameters
  • apr – apr file (Globk format)

  • epoch – epoch in decimal year for coordinates in apr

  • verbose – verbose mode (boolean)

  • exluded_periods – periods to be excluded

Returns

Gts instance

Note

following Globk’s convention, site will be named XXXX_1PS, XXXX_2PS etc between offset dates

pyacs.gts.lib.metadata.read_offset_dates(self, offset_file)[source]

Reads an offset file and populates offsets_dates (pyacs format) attribute of the current Gts instance. format is simply a code dates. dates can be any format read by pyacs.guess_date

Parameters

offset_file – offset_file to be read

Returns

the current Gts instance

pyacs.gts.lib.metadata.info(self, info=2)[source]

Print various informations about the time series

pyacs.gts.lib.noise module

Pyacs noise module This is a wrapper to Globk tsfit and the CATS Software Reference Herring, T. A., King, R. W., Floyd, M. A. & McClusky, S. C. (2015). GAMIT/GLOBK Reference manual, 10.6. Williams, S. D. (2008). CATS: GPS coordinate time series analysis software. GPS solutions, 12(2), 147-153.

pyacs.gts.lib.noise.wrms(self)[source]

Return the wrms :return wrms: return(np.array([wrms_n,wrms_e,wrms_up]))

pyacs.gts.lib.noise.realistic_sigma(self, option='tsfit', in_place=False, verbose=False)[source]

Calculates realistic sigmas on velocity components :param option:

  • tsfit: globk T. Herring realistic sigma

  • cats_pl: CATS estimates with noise type estimated (i.e. –model=pl:)

  • cats_seasonal_pl: CATS estimates with seasonal terms and noise type estimated (i.e. –model=pl: –sinusoid=1y1)

  • cats_flicker: CATS estimates assuming flicker noise (i.e. –model=pl:k-1)

  • cats_seasonal_flicker: CATS estimates with seasonal terms and assuming flicker noise (i.e. –model=pl:k-1 –sinusoid=1y1)

pyacs.gts.lib.noise.sigma_cats(self, in_place=False, verbose=False, k='k-1', seasonal='')[source]

runs CATS for getting realistic sigma

pyacs.gts.lib.noise.sigma_vel_tsfit(self, in_place=False, verbose=False)[source]

runs tsfit for getting realistic sigma

pyacs.gts.lib.noise.get_spectral_index(cats_file)[source]
pyacs.gts.lib.noise.get_vel_and_sigma(cats_file)[source]

pyacs.gts.lib.offset module

pyacs.gts.lib.offset.print_no_return(str)[source]
pyacs.gts.lib.offset.suspect_offsets_mf(self, threshold=3, verbose=True, lcomponent='NE', n_max_offsets=5, in_place=False, debug=False)[source]

Tries to find offsets in a time series using a median filter

pyacs.gts.lib.offset.test_offset_significance(self, date, conf_level=95, lcomponent='NE', verbose=True, debug=False, mode='local')[source]

test the significance of an offset

Param

date : date of the offset to be tested

Param

conf_level : confidence level in percent (default=95)

Param

lcomponent : component to be tested (default=’NE’)

Param

mode : choose among ‘local’,’detrend’,’detrend_seasonal’ to test significance

Param

verbose : verbose mode

Returns

: True if significant, else False

pyacs.gts.lib.offset.local_offset_robust(self, date, n, verbose=False, debug=False)[source]

estimate a local offset (no velocity) with a robust method :param date: date in decimal year :param n: number of dates before and after the dates used in the estimation :return : a 1D numpy array with [date, north, east, up, s_north, s_east, s_up] :note: the offsets are estimated using the difference between the median position before and after the earthquake using i days for all i <=n. Then the median of the estimates is returned.

pyacs.gts.lib.offset.apply_offsets(self, np_offset, opposite=False, in_place=False, verbose=False)[source]

Applies given offsets to a times series np_offset is a 1D np.array with lines [dates,north,east,up] if in_place = True then replace the current time series

Parameters
  • np_offset – 1D or 3D numpy array or list or list of list with column offset_dates, north, east, up, s_north, s_east, s_up

  • opposite – boolean, if True apply the oppsite of provided offsets

  • in_place – if True, will make change in place, if False, return s a new time series

  • verbose – boolean, verbose mode

pyacs.gts.lib.offset.find_offsets_t_scan(self, threshold=0.8, window=250, in_place=False, lcomponent='NE', verbose=True, debug=True)[source]
pyacs.gts.lib.offset.get_suspected_dates(diff_data, threshold, lcomponent='NEU', verbose=False)[source]

get the list of the largest values; these are suspected offsets

pyacs.gts.lib.offset.check_suspected_offsets(lindex, verbose=False)[source]

Check that the list of suspected index does not contain two successive values. In this case, it is certainly an isolated outlier and the suspected offsets are removed from the list.

pyacs.gts.lib.offset.find_offsets(self, threshold=3, n_max_offsets=9, conf_level=95, lcomponent='NE', verbose=True, in_place=False)[source]

A simple empirical procedure to find offsets.

Parameters
  • threshold – threshold value for offset preliminary detection

  • n_max_offset – maximum number of offsets to be detected simultaneously

  • conf_level – confidence level for a suspected offset to be accepted

  • lcomponent – components used for offset detection

Returns

a new Gts instance with offsets_dates and outliers now populated

pyacs.gts.lib.offset.suspect_offsets(self, threshold=3, verbose=True, lcomponent='NE', n_max_offsets=10, in_place=False)[source]

Tries to find offsets in a time series

pyacs.gts.lib.offset.find_time_offsets(self, option=None, ndays=7, th_detection_rms=3, th_detection_offset=3)[source]

Find the time of suspected offsets by rms time series calculated over ndays Then check the time of offsets: if one offset is too small/None then it is removed input:

  • ndays: number of positions in the windows. rms time series are calculated over ndays.

  • th_detection_rms: the threshold for detecting the anomalous windows rms(t) > th_detection_rms*median(rms(ts)).

  • th_detection_offset: the threshold for detecting the offsets,

    for each anomalous time windows, differentiate positions then test whether it is a suspected offset (differentiated(t) > threshold * median(differentiated))

output:

add the time of offsets in to self.offsets

pyacs.gts.lib.offset.delete_small_offsets(self, offsets, del_by_pricise=False)[source]

The aim for test_offset modul. Estimate the offsets with clean data. Then delete the offsets which their values are so small input: list of time offsets output: list of time offsets tested

pyacs.gts.lib.offset.test_offsets(self, verbose=False, debug=True, window_length=None)[source]
Test the offset:
  • delete the small offset (1mm for East/North, 2mm for Up)

  • then make a F ratio test

  • re-check to delete the small offsets

pyacs.gts.lib.offset.estimate_local_offset(self, window_length=4, in_place=False)[source]

Estimate the local offset, just used window_length positions before & window_length positions behind of offset output: amplitude of local offsets

pyacs.gts.lib.outliers_old module

pyacs.gts.lib.outliers_old.remove_outliers(self, periods=None, in_place=False)[source]

removes outliers provided in self.outliers return a new Gts without the outliers if in_place = True then self has the outliers removed as well (in _place)

pyacs.gts.lib.outliers_old.find_outlier_around_date(self, date, conf_level=95, n=3, lcomponent='NE', verbose=True)[source]

Find an outlier around a given date returns the index of the outlier, returns [] if no outlier found :param date : given date :param conf_level : confidence level for F_ratio test of outlier significance (default 95%%) :param n : number of dates either sides of date (default n=3) :param lcomponent : components ‘N’,’E’,’U’,’NE’,’NEU’ (default ‘NE’)

pyacs.gts.lib.outliers_old.find_outliers_percentage(self, percentage=0.03, in_place=False, verbose=False, component='NEU', periods=None, excluded_periods=None)[source]

detrend a time series and ranks the residuals by increasing absolute value populate the outliers with the x % largest ones on each component

pyacs.gts.lib.outliers_old.find_outliers_simple(self, threshold=100, window_length=10, in_place=False, verbose=False, component='NEU', periods=None, excluded_periods=None)[source]
pyacs.gts.lib.outliers_old.find_outliers_and_offsets_through_differentiation(self, th=100)[source]

find outliers and offsets using differenciation

pyacs.gts.lib.outliers_old.find_outliers_by_RMS_ts(self, ndays=7, th_detection=5, th_rejection=2)[source]
Find index of outliers in a time series, populate self.outliers.
  • rms time series are first calculated over ndays

  • time windows are kept for further inspection if rms(t) > th_detection * median(rms(ts))

  • for each anomalous time windows, differentiate positions, find the max

  • test whether it is a true outlier (differentiated(t) > th_rejection * median(differentiated))

output:

None

pyacs.gts.lib.outliers_old.find_outliers_by_residuals(self, threshold=5, model='detrend_seasonal', component='NE', in_place=False)[source]

Find index of outliers by trendline/trendline_annual/trendline_seasonal (the complete model) Then the outliers are detected if their residuals are greater than th_rejection*standard_deviation

output:

Add the list of outlier index into self.outliers

pyacs.gts.lib.outliers_old.find_outliers_vondrak(self, threshold=10, fc=2.0, in_place=False, verbose=True, periods=[[]], excluded_periods=[[]], component='NE')[source]

Find outliers using a Vondrak filter

pyacs.gts.lib.outliers_old.find_outliers_sliding_window(self, threshold=3, in_place=False, verbose=True, periods=[[]], excluded_periods=[[]], component='NE', window_len=15, automatic=True)[source]

Find outliers using sliding windows

pyacs.gts.lib.step_detect module

Thomas Kahn thomas.b.kahn@gmail.com

pyacs.gts.lib.step_detect.t_scan(L, window=1000.0, num_workers=- 1)[source]

Computes t statistic for i to i+window points versus i-window to i points for each point i in input array. Uses multiple processes to do this calculation asynchronously. Array is decomposed into window number of frames, each consisting of points spaced at window intervals. This optimizes the calculation, as the drone function need only compute the mean and variance for each set once. Parameters ———- L : numpy array

1 dimensional array that represents time series of datapoints

windowint / float

Number of points that comprise the windows of data that are compared

num_workersint

Number of worker processes for multithreaded t_stat computation Defult value uses num_cpu - 1 workers

t_statnumpy array

Array which holds t statistic values for each point. The first and last (window) points are replaced with zero, since the t statistic calculation cannot be performed in that case.

pyacs.gts.lib.step_detect.mz_fwt(x, n=2)[source]

Computes the multiscale product of the Mallat-Zhong discrete forward wavelet transform up to and including scale n for the input data x. If n is even, the spikes in the signal will be positive. If n is odd the spikes will match the polarity of the step (positive for steps up, negative for steps down). This function is essentially a direct translation of the MATLAB code provided by Sadler and Swami in section A.4 of the following: http://www.dtic.mil/dtic/tr/fulltext/u2/a351960.pdf Parameters ———- x : numpy array

1 dimensional array that represents time series of data points

nint

Highest scale to multiply to

prodnumpy array

The multiscale product for x

pyacs.gts.lib.step_detect.find_steps(array, threshold)[source]

Finds local maxima by segmenting array based on positions at which the threshold value is crossed. Note that this thresholding is applied after the absolute value of the array is taken. Thus, the distinction between upward and downward steps is lost. However, get_step_sizes can be used to determine directionality after the fact. Parameters ———- array : numpy array

1 dimensional array that represents time series of data points

thresholdint / float

Threshold value that defines a step

stepslist

List of indices of the detected steps

pyacs.gts.lib.step_detect.get_step_sizes(array, indices, window=1000)[source]

Calculates step size for each index within the supplied list. Step size is determined by averaging over a range of points (specified by the window parameter) before and after the index of step occurrence. The directionality of the step is reflected by the sign of the step size (i.e. a positive value indicates an upward step, and a negative value indicates a downward step). The combined standard deviation of both measurements (as a measure of uncertainty in step calculation) is also provided. Parameters ———- array : numpy array

1 dimensional array that represents time series of data points

indiceslist

List of indices of the detected steps (as provided by find_steps, for example)

windowint, optional

Number of points to average over to determine baseline levels before and after step.

step_sizeslist

List of the calculated sizes of each step

step_error : list

pyacs.gts.lib.step_detect_edge_filter module

pyacs.gts.lib.step_detect_edge_filter.find_offsets_edge_filter(self, 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)[source]

Module contents