Tidal Calibration
This notebook is intended to demonstrate the tidal calibration process for Gladwin Tensor borehole Strainmeters (GTSM) following the methodology outlined in Hodgkinson et al. (2013) and applied in Hanagan et al. (In Prep). Note that tidal calibration requires a careful eye on the data for processing; so, when available, we suggest using the tidal calibrations provided in the metadata unless the user has a strong reason to belive the calibration requires an update.
The workflow in this notebook can be summarized into three main parts: data aquisition, offset correction, and tidal time series analysis and modelling to acquire the final tidal calibration matrix. An additional benefit of the workflow implemented here is that we will also estimate the gauge response to barometric pressure. We correct the data for spurious offsets or pulses prior to tidal analysis because they are known to hinder accurate characterization of the tides and barometric pressure repsonse.
The first portion of this notebook includes steps from the level 2 processing notebook. Then, new steps specific to the calibraiton process are completed. Sections in this notebook include:
Aquire station metadata and decimated timeseries
Offset correction
Observed tidal analysis
Predicted tides
Orientation matrix
Aquire station metadata and decimated timeseries
[1]:
# Import relevant modules from the earscopestraintools package
from earthscopestraintools.mseed_tools import ts_from_mseed
from earthscopestraintools.gtsm_metadata import GtsmMetadata
from earthscopestraintools.timeseries import plot_timeseries_comparison
import numpy as np
import pandas as pd
# Allow logged output to be printed in the notebook as code cells are run
import logging
logger = logging.getLogger()
logging.basicConfig(
format="%(message)s", level=logging.INFO
)
[2]:
# Linearize 1 hz strain from the IRIS/EarthScope DMC, filter/decimate to 5 min
# and acquire 30 min pressure data
# Metadata
network = 'IV'
station = 'TSM6'
meta = GtsmMetadata(network,station)
# Provide the start and end times for calibration
start = "2022-10-21T00:00:00.000000"
end = "2023-03-09T23:59:59.000000"
# load data
strain_raw = ts_from_mseed(network=network, station=station, location='T0', channel='LS*', start=start, end=end)
# Print stats and plot the data
strain_raw.stats()
strain_raw.plot(type='line')
IV TSM6 Loading T0 LS* from 2022-10-21T00:00:00.000000 to 2023-03-09T23:59:59.000000 from Earthscope DMC miniseed
Trace 1. 2022-10-21T00:00:00.000000Z:2023-03-09T23:59:59.000000Z mapping LS1 to CH0
Trace 2. 2022-10-21T00:00:00.000000Z:2023-03-09T23:59:59.000000Z mapping LS2 to CH1
Trace 3. 2022-10-21T00:00:00.000000Z:2023-03-09T23:59:59.000000Z mapping LS3 to CH2
Trace 4. 2022-10-21T00:00:00.000000Z:2023-03-09T23:59:59.000000Z mapping LS4 to CH3
Found 0 epochs with nans, 29.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Converting missing data from 999999 to nan
Converting 999999 values to nan
Found 29 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
IV.TSM6.T0.LS*
| Channels: ['CH0', 'CH1', 'CH2', 'CH3']
| TimeRange: 2022-10-21 00:00:00 - 2023-03-09 23:59:59 | Period: 1s
| Series: raw| Units: counts| Level: 0| Gaps: 0.0%
| Epochs: 12096000| Good: 12095971.0| Missing: 29.0| Interpolated: 0.0
| Samples: 48384000| Good: 48383884| Missing: 116| Interpolated: 0
[3]:
# Filter to 5 min, then to hourly in a 2 step filter and convert to nstrain
# Save the hourly data in a separate Timeseries with an index in datetime seconds
# filter and decimate to 5 min
decimated_counts = strain_raw.decimate_1s_to_300s(method='linear',limit=3600)
filt_cutoff_s = 2*60*60 # 2 hr cutoff period
# Applies a lowpass butterworth filter via the scipy butterworth filter function
filtered_hr_counts = decimated_counts.butterworth_filter(filter_type='lowpass',filter_order=5,
filter_cutoff_s=filt_cutoff_s)
# Decimates the data to hourly by taking the first minute/second of each hour
decimated_hr_counts = filtered_hr_counts.decimate_to_hourly()
# Linearize to microstrain
decimated_hr_strain = decimated_hr_counts.linearize(reference_strains=meta.reference_strains,gap=meta.gap)
Decimating to 300s
Interpolating data using method=linear and limit=3600
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Applying Butterworth Filter
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Decimating to hourly
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Converting raw counts to microstrain
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Offset correction
We will correct for offsets via a simple first differencing method. To enable better estimates of the offsets, we first correct the data for a “first estimate” of tides, trend, and pressure. For this time series, we use Baytap08 (more information on the program is described in the next section) to retreive an initial estimate of the tides and pressure response for the purpose of preliminary correction.
[4]:
# Acquire atmospheric pressure data and filter to hourly
# Pressure data from miniseed
atmp_raw = ts_from_mseed(network=network, station=station, location='*', channel='RDO',
start=start, end=end, period=60*30, scale_factor=0.001, units='hpa')
# Interpolate, filter and decimate to hourly
atmp_interp = atmp_raw.interpolate()
atmp_filt = atmp_interp.butterworth_filter(filter_type='lowpass',filter_order=5,
filter_cutoff_s=filt_cutoff_s)
atmp = atmp_filt.decimate_to_hourly()
# Run baytap preliminarily to retrieve the pressure response and tidal amplitudes and phases for the full metadata suite
prelim_baytap_results = decimated_hr_strain.baytap_analysis(atmp,latitude=meta.latitude,longitude=meta.longitude,elevation=meta.elevation,dmin=0.001)
IV TSM6 Loading * RDO from 2022-10-21T00:00:00.000000 to 2023-03-09T23:59:59.000000 from Earthscope DMC miniseed
Trace 1. 2022-10-21T00:00:00.000000Z:2023-03-09T23:30:00.000000Z mapping RDO to atmp
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Converting missing data from 999999 to nan
Converting 999999 values to nan
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Interpolating data using method=linear and limit=2
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Applying Butterworth Filter
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Decimating to hourly
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
bda488faa0b854188cbc29c5b15092d7d3f2a114022340fb6414b9875dab3f21
Docker container started.
Atmospheric pressure responses in microstrain/hPa) and tidal parameters in degrees/nanostrain
baytap
Docker processes finished. Container removed.
[5]:
# Corrections prior to offset correcteion
atmp_c = atmp.calculate_pressure_correction(response_coefficients=prelim_baytap_results['atmp_response'])
trend_c = decimated_hr_strain.linear_trend_correction(method='median')
tide_c = decimated_hr_strain.calculate_tide_correction(tidal_parameters=prelim_baytap_results['tidal_params'],longitude=meta.longitude)
pre_offset_corrected = decimated_hr_strain.apply_corrections([trend_c,tide_c,atmp_c])
offset_c = pre_offset_corrected.calculate_offsets(limit_multiplier=10)
offset_corrected = decimated_hr_strain.apply_corrections([offset_c])
offset_c.plot()
plot_timeseries_comparison([decimated_hr_strain,offset_corrected],names=['Original','Offset Corrected'],zero=True,detrend='linear')
Calculating pressure correction
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Calculating linear trend correction
Trend Start: 2022-10-21 00:00:00
Trend End: 2023-03-09 23:00:00
Median trend calculated with points 24.0 hr apart.
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Calculating tide correction
WARNING: The requested image's platform (linux/amd64) does not match the detected host platform (linux/arm64/v8) and no specific platform was requested
WARNING: The requested image's platform (linux/amd64) does not match the detected host platform (linux/arm64/v8) and no specific platform was requested
WARNING: The requested image's platform (linux/amd64) does not match the detected host platform (linux/arm64/v8) and no specific platform was requested
WARNING: The requested image's platform (linux/amd64) does not match the detected host platform (linux/arm64/v8) and no specific platform was requested
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Applying corrections
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Calculating offsets using cutoff percentile of 0.75 and limit multiplier of 10.
Using offset limits of [0.008823, 0.013784, 0.008222, 0.008635]
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Applying corrections
Found 0 epochs with nans, 0.0 epochs with 999999s, and 0 missing epochs.
Total missing data is 0.0%
Observed tidal analysis
The tidal constituents (amplitudes and phases of daily and bi-daily tides) will be determined from the observed strain timeseries in a program called BAYTAP08 (Tamura 1991; Tamura and Agnew 2008). The original BAYTAP-G program was rewritten by Duncan Agnew as BAYTAP08: http://igppweb.ucsd.edu/~agnew/Baytap/baytap.html. The primary function of BAYTAP08 is to isolate distinct tidal amplitudes and phases from data streams contaminated by other signals through assuming a smoothly varying drift, which is empirically inferred by minimizing Akaike’s Bayesian Information Criterion (ABIC).
To summarize, the program assumes a strain signal (\(y_i\)) composed of the following components:
\[y_i = t_i + d_i + c_i + s_i\]
\(t_i\): tidal signal - \(d_i\): long term trend (instrument drift) - \(c_i\): response to other effects (e.g., barometric pressure) - \(s_i\): data offsets
The tidal signal is composed of the tidal amplitudes and phases we seek (\(A_m\) and \(B_m\)), and the known theoretical tidal group values (\(C_{mi}\) and \(S_{mi}\)) that are the sum of the constituents with similar frequencies: >
The barometric response is considered linear with the barometric pressure change. Where, \(b_i\) is the barometric response coefficient: >
The residual time series is left once these signals are removed: >
BAYTAP solves for the model parameters by minimizing S: >
D (drift) is an input smoothness parameter, where a very large D provides a linear drift with time. - W is an input parameter that controls how much the tidal admittance can vary over frequency bands. We allow W to assume the default value of 1.
Typically, the strainmeters have already been analyzed with BAYTAP, and those outputs are what we used to perform the preliminary corrections prior to detecting offsets and pulses in the previous step. Here, the goal is to reanalyze the tidal amplitudes and phases to perform tidal calibration. More specifically, we are interested in the M2 and O1 tidal constituent amplitudes and phases per gauge.
[6]:
# Run baytap to retrieve the pressure response and tidal amplitudes and phases for the full metadata suite
baytap_results= offset_corrected.baytap_analysis(atmp,latitude=meta.latitude,longitude=meta.longitude,elevation=meta.elevation,dmin=0.001)
55222a7e1f6c1757836c78145af65a76f681e91a927ac93e052cc8c4bdc130eb
Docker container started.
Atmospheric pressure responses in microstrain/hPa) and tidal parameters in degrees/nanostrain
baytap
Docker processes finished. Container removed.
[7]:
# Print the results for the M2 and O1 tidal constituents
for ch in offset_corrected.columns:
print(f'{ch}:')
print('M2 amp (nstrain) and phase (degree):',
baytap_results['tidal_params'][(ch,'M2','amp')],baytap_results['tidal_params'][(ch,'M2','phz')])
print('O1 amp (nstrain) and phase (degree):',
baytap_results['tidal_params'][(ch,'O1','amp')],baytap_results['tidal_params'][(ch,'O1','phz')])
CH0:
M2 amp (nstrain) and phase (degree): 5.941 -46.037
O1 amp (nstrain) and phase (degree): 3.983 132.863
CH1:
M2 amp (nstrain) and phase (degree): 8.598 63.713
O1 amp (nstrain) and phase (degree): 4.622 -116.351
CH2:
M2 amp (nstrain) and phase (degree): 13.152 -165.709
O1 amp (nstrain) and phase (degree): 2.564 10.488
CH3:
M2 amp (nstrain) and phase (degree): 8.066 -120.395
O1 amp (nstrain) and phase (degree): 3.702 60.080
Predicted tides
In this section, we will gather the predicted M2 and O1 tidal amplitudes and phases for the areal, differential shear, and engineering shear strains for the station using SPOTL (Some Programs for Ocean Tide Loading; Agnew, 2013). SPOTL, like BAYTAP, will be run from a Docker container.
[8]:
from earthscopestraintools.processing import spotl_predict_tides
# Parameters for the tidal model
latitude= meta.latitude
longitude = meta.longitude
elevation = meta.elevation
# Regional ocean load model, global ocean load model, and Green's functions to decribe the elastic Earth
# available from the SPOTL package
reg_oc = 'osu.mediterranean.2011'
glob_oc= 'osu.tpxo72.2010' # change to osu.tpxo72atlas.2011 ?
greenf = 'green.contap.std'
pred_tides = spotl_predict_tides(latitude,longitude,elevation,glob_oc,reg_oc,greenf)
pred_tides
Docker started
WARNING: The requested image's platform (linux/amd64) does not match the detected host platform (linux/arm64/v8) and no specific platform was requested
spotl
Docker container stopped and removed.
[8]:
{'M2': {'areal': (15.260654206772385+1.559118536665793j),
'differential': (-10.46438760598575-3.393253671866749j),
'engineering': (0.9177684867587321-5.360603600781594j)},
'O1': {'areal': (11.46017566155345+0.08470101319819293j),
'differential': (2.8783089246691294+0.03687517960914067j),
'engineering': (-0.24828220829683426+3.638939777606006j)}}
Orientation matrix
The orientation matrix is determined through simple least squares with the M2 and O1 tidal results in the code cells below, as described in Hodgkinson et al., (2013). Essentially, the orientation matrix that transforms the observed tidal to match the predicted tides is determined.
The forward problem, with gauge strains \(e_i\), regional strains \(E_j\), and the coupling matrix \(A_{ij}\) has the form:
\[e_{i} = A_{ij} E_{j}\]\[\begin{split}\begin{pmatrix} e_{0} \\ e_{1} \\ e_{2} \\ e_{3} \end{pmatrix} = \begin{pmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \end{pmatrix} \begin{pmatrix} e_{EE} + e_{NN} \\ e_{EE} - e_{NN} \\ 2e_{EN} \end{pmatrix} = \begin{pmatrix} E_A \\ E_D \\ E_S \end{pmatrix}\end{split}\]
The orientation matrix is the Moore-Penrose pseudoinverse of the coupling matrix, and is applied to the gauge strains as follows to attain the regional areal, differential, and engineering shear strains:
\[A^{-1}_{ji} e_{i} = E_{j}\]\[\begin{split}\begin{pmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \end{pmatrix}^{-1} \begin{pmatrix} e_{0} \\ e_{1} \\ e_{2} \\ e_{3} \end{pmatrix} = \begin{pmatrix} e_{EE} + e_{NN} \\ e_{EE} - e_{NN} \\ 2e_{EN} \end{pmatrix} = \begin{pmatrix} E_A \\ E_D \\ E_S \end{pmatrix}\end{split}\]
[9]:
# Compute the observed real and imaginary components of the tidal amplitudes and phases,
# and place in a 4x4 matrix
def make_complex(ampstr, phzstr):
amp, phz = float(ampstr), float(phzstr)
return complex(amp*np.cos(phz*np.pi/180),amp*np.sin(phz*np.pi/180))
obs_tide_matrix = np.zeros([4,4])
for i,ch in enumerate(offset_corrected.columns):
m2 = make_complex(baytap_results['tidal_params'][(ch,'M2','amp')], baytap_results['tidal_params'][(ch,'M2','phz')])
o1 = make_complex(baytap_results['tidal_params'][(ch,'O1','amp')], baytap_results['tidal_params'][(ch,'O1','phz')])
obs_tide_matrix[i,:] = np.real(m2), np.imag(m2), np.real(o1), np.imag(o1)
print('Observed tidal matrix:')
obs_tide_matrix
Observed tidal matrix:
[9]:
array([[ 4.12420475, -4.27626194, -2.70942648, 2.91946865],
[ 3.80777711, 7.70885448, -2.05156249, -4.14173579],
[-12.74500502, -3.24652907, 2.52116339, 0.46672387],
[ -4.08106121, -6.95739142, 1.84652174, 3.2086074 ]])
[10]:
# For the predicted tide matrix, 3x4
pred_tide_matrix = np.array([[np.real(pred_tides['M2']['areal']),np.imag(pred_tides['M2']['areal']),
np.real(pred_tides['O1']['areal']),np.imag(pred_tides['O1']['areal'])],
[np.real(pred_tides['M2']['differential']),np.imag(pred_tides['M2']['differential']),
np.real(pred_tides['O1']['differential']),np.imag(pred_tides['O1']['differential'])],
[np.real(pred_tides['M2']['engineering']),np.imag(pred_tides['M2']['engineering']),
np.real(pred_tides['O1']['engineering']),np.imag(pred_tides['O1']['engineering'])]])
print('Predicted tidal matrix:')
pred_tide_matrix
Predicted tidal matrix:
[10]:
array([[ 15.26065421, 1.55911854, 11.46017566, 0.08470101],
[-10.46438761, -3.39325367, 2.87830892, 0.03687518],
[ 0.91776849, -5.3606036 , -0.24828221, 3.63893978]])
[11]:
# Calculate the orentation matrix from the coupling matrix
coupling_mat = np.linalg.lstsq(pred_tide_matrix.T,obs_tide_matrix.T,rcond=-1)[0].T
res2 = np.linalg.lstsq(pred_tide_matrix.T,obs_tide_matrix.T,rcond=-1)[1]
orient_mat = np.linalg.pinv(coupling_mat)
print('Coupling Matrix:\n',coupling_mat)
print('****************************')
print('Tidal Orientation Matrix:\n',orient_mat)
print('****************************')
print('Residuals (obs-modelled)')
print(obs_tide_matrix - np.dot(coupling_mat,pred_tide_matrix))
# Sum of squared residuals,
# divided by number of data (4 per gauge)
print('RMSE, per gauge:\n',np.sqrt(res2/4))
print('RMSE:\n',np.sqrt(sum(res2)/16))
Coupling Matrix:
[[-0.09759912 -0.44106484 0.9756906 ]
[-0.06350349 -0.55316905 -1.11388831]
[-0.06603285 1.10820783 -0.04401844]
[ 0.04723211 0.54496425 0.93881452]]
****************************
Tidal Orientation Matrix:
[[-4.8884553 -2.59321853 -4.13009602 1.81001004]
[-0.27530833 -0.13289096 0.6483382 0.1588481 ]
[ 0.41141518 -0.31490917 -0.12047325 0.25831448]]
****************************
Residuals (obs-modelled)
[[ 0.10269965 -0.39044769 -0.07915587 -0.60647961]
[ 0.0105982 -0.04029265 -0.00816857 -0.06258628]
[-0.10018539 0.38088889 0.07721801 0.59163199]
[ 0.03924873 -0.14921743 -0.03025101 -0.23177836]]
RMSE, per gauge:
[0.36642864 0.03781398 0.35745786 0.14003806]
RMSE:
0.26602942711926336
References
Agnew, D. C. (2012). SPOTL: Some Programs for Ocean-Tide Loading, SIO Technical Report, Scripps Institution of Oceanography. From https://igppweb.ucsd.edu/~agnew/Spotl/spotlmain.html
Hodgkinson, K., J. Langbein, B. Henderson, D. Mencin, and A. Borsa (2013), Tidal calibration of plate boundary observatory borehole strainmeters, J. Geophys. Res. Solid Earth, 118, 447–458, doi:10.1029/2012JB009651.
Tamura, Y., & Agnew, D. C. (2008). Baytap08 User’s Manual. UC San Diego: Library – Scripps Digital Collection. From https://escholarship.org/uc/item/4c27740c