Compute the H/V spectral ratio from the imaginary part of the CCFs

Import & getting the data from the computed autocorrelations (ZZ, EE, NN)

import os
if "SPHINX_DOC_BUILD" in os.environ:
    if "MSNOISE_DOC" in os.environ:
        os.chdir(os.environ["MSNOISE_DOC"])

import matplotlib
matplotlib.use("agg")

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pandas.plotting import register_matplotlib_converters
import datetime
register_matplotlib_converters()

plt.style.use("ggplot")

from msnoise.api import connect, get_results, build_movstack_datelist, get_params, get_t_axis, xr_get_ccf


# connect to the database
db = connect()

# Obtain a list of dates between ``start_date`` and ``enddate``
start, end, datelist = build_movstack_datelist(db)

# Get the list of parameters from the DB:
params = get_params(db)

# Get the time axis for plotting the CCF:
taxis = get_t_axis(db)

# Get the first mov_stack configured:
mov_stack = params.mov_stack[0]

# Get the results for two station, filter id=1, ZZ component, mov_stack=1 and the results as a 2D array:
ZZ = xr_get_ccf("PF.FJS.00", "PF.FJS.00", "ZZ", 2, mov_stack, taxis, format="dataframe")
EE = xr_get_ccf("PF.FJS.00", "PF.FJS.00", "EE", 2, mov_stack, taxis, format="dataframe")
NN = xr_get_ccf("PF.FJS.00", "PF.FJS.00", "NN", 2, mov_stack, taxis, format="dataframe")

Checking the data has the same time base (index) and joining

r = "1D"
rZZ = ZZ.resample(r).mean()
rEE = EE.resample(r).mean()
rNN = NN.resample(r).mean()

merged = pd.concat({'ZZ': rZZ, 'EE': rEE, 'NN': rNN}, axis=0)
merged.index.names = ['channel', 'date']

# Swap the levels of the MultiIndex
result = merged.swaplevel('channel', 'date').sort_index(level='date') #.iloc[:500]

Helper functions

def get_imag(d, fs):
    NFFT = 1024*32
    iX = np.imag(np.fft.fft(d, n=NFFT)[:NFFT//2])
    freqs = np.fft.fftfreq(NFFT, d=1/fs)[:NFFT//2]

    return freqs, iX

Doing the job !

HVSRs = {}
for date, group in result.groupby(level='date'):
    print(f"Date: {date}")
    group= group.droplevel(0)
    try:
        f, iZZ = get_imag(group.loc["ZZ"].values, 20)
        f, iEE = get_imag(group.loc["EE"].values, 20)
        f, iNN = get_imag(group.loc["NN"].values, 20)
    except:
        continue
    hvsr = np.sqrt((iEE+iNN)/iZZ)
    HVSRs[date] = hvsr

HVSR = pd.DataFrame(HVSRs, index=f)
X = HVSR.copy().fillna(0.0)
X = X.T.resample(r).mean().T

Date: 2014-03-21 00:00:00
C:\Users\thlec\PycharmProjects\MSNoise\examples\plot_compute_hvsr.py:91: RuntimeWarning: invalid value encountered in divide
  hvsr = np.sqrt((iEE+iNN)/iZZ)
C:\Users\thlec\PycharmProjects\MSNoise\examples\plot_compute_hvsr.py:91: RuntimeWarning: invalid value encountered in sqrt
  hvsr = np.sqrt((iEE+iNN)/iZZ)
Date: 2014-03-22 00:00:00
Date: 2014-03-23 00:00:00
Date: 2014-03-24 00:00:00
Date: 2014-03-25 00:00:00
C:\Users\thlec\PycharmProjects\MSNoise\examples\plot_compute_hvsr.py:91: RuntimeWarning: divide by zero encountered in divide
  hvsr = np.sqrt((iEE+iNN)/iZZ)
Date: 2014-03-26 00:00:00
Date: 2014-03-27 00:00:00
Date: 2014-03-28 00:00:00
Date: 2014-03-29 00:00:00
Date: 2014-03-30 00:00:00
Date: 2014-03-31 00:00:00
Date: 2014-04-01 00:00:00
Date: 2014-04-02 00:00:00
Date: 2014-04-03 00:00:00
Date: 2014-04-04 00:00:00
Date: 2014-04-05 00:00:00
Date: 2014-04-06 00:00:00
Date: 2014-04-07 00:00:00
Date: 2014-04-08 00:00:00
Date: 2014-04-09 00:00:00
Date: 2014-04-10 00:00:00
Date: 2014-04-11 00:00:00
Date: 2014-04-12 00:00:00
Date: 2014-04-13 00:00:00
Date: 2014-04-14 00:00:00
Date: 2014-04-15 00:00:00
Date: 2014-04-16 00:00:00
Date: 2014-04-17 00:00:00
Date: 2014-04-18 00:00:00
Date: 2014-04-19 00:00:00
Date: 2014-04-20 00:00:00
Date: 2014-04-21 00:00:00
Date: 2014-04-22 00:00:00
Date: 2014-04-23 00:00:00
Date: 2014-04-24 00:00:00
Date: 2014-04-25 00:00:00
Date: 2014-04-26 00:00:00
Date: 2014-04-27 00:00:00
Date: 2014-04-28 00:00:00
Date: 2014-04-29 00:00:00
Date: 2014-04-30 00:00:00
Date: 2014-05-01 00:00:00
Date: 2014-05-02 00:00:00
Date: 2014-05-03 00:00:00
Date: 2014-05-04 00:00:00
Date: 2014-05-05 00:00:00
Date: 2014-05-06 00:00:00
Date: 2014-05-07 00:00:00
Date: 2014-05-08 00:00:00
Date: 2014-05-09 00:00:00
Date: 2014-05-10 00:00:00
Date: 2014-05-11 00:00:00
Date: 2014-05-12 00:00:00
Date: 2014-05-13 00:00:00
Date: 2014-05-14 00:00:00
Date: 2014-05-15 00:00:00
Date: 2014-05-16 00:00:00
Date: 2014-05-17 00:00:00
Date: 2014-05-18 00:00:00
Date: 2014-05-19 00:00:00
Date: 2014-05-20 00:00:00
Date: 2014-05-21 00:00:00
Date: 2014-05-22 00:00:00
Date: 2014-05-23 00:00:00
Date: 2014-05-24 00:00:00
Date: 2014-05-25 00:00:00
Date: 2014-05-26 00:00:00
Date: 2014-05-27 00:00:00
Date: 2014-05-28 00:00:00
Date: 2014-05-29 00:00:00
Date: 2014-05-30 00:00:00
Date: 2014-05-31 00:00:00
Date: 2014-06-01 00:00:00
Date: 2014-06-02 00:00:00
Date: 2014-06-03 00:00:00
Date: 2014-06-04 00:00:00
Date: 2014-06-05 00:00:00
Date: 2014-06-06 00:00:00
Date: 2014-06-07 00:00:00
Date: 2014-06-08 00:00:00
Date: 2014-06-09 00:00:00
Date: 2014-06-10 00:00:00
Date: 2014-06-11 00:00:00
Date: 2014-06-12 00:00:00
Date: 2014-06-13 00:00:00
Date: 2014-06-14 00:00:00
Date: 2014-06-15 00:00:00
Date: 2014-06-16 00:00:00
Date: 2014-06-17 00:00:00
Date: 2014-06-18 00:00:00
Date: 2014-06-19 00:00:00
Date: 2014-06-20 00:00:00
Date: 2014-06-21 00:00:00
Date: 2014-06-22 00:00:00
Date: 2014-06-23 00:00:00
Date: 2014-06-24 00:00:00
Date: 2014-06-25 00:00:00
Date: 2014-06-26 00:00:00
Date: 2014-06-27 00:00:00
Date: 2014-06-28 00:00:00
Date: 2014-06-29 00:00:00
Date: 2014-06-30 00:00:00
Date: 2014-07-01 00:00:00
Date: 2014-07-02 00:00:00
Date: 2014-07-03 00:00:00
Date: 2014-07-04 00:00:00
Date: 2014-07-05 00:00:00
Date: 2014-07-06 00:00:00
Date: 2014-07-07 00:00:00
Date: 2014-07-08 00:00:00
Date: 2014-07-09 00:00:00
Date: 2014-07-10 00:00:00
Date: 2014-07-11 00:00:00
Date: 2014-07-12 00:00:00
Date: 2014-07-13 00:00:00
Date: 2014-07-14 00:00:00
Date: 2014-07-15 00:00:00
Date: 2014-07-16 00:00:00
Date: 2014-07-17 00:00:00
Date: 2014-07-18 00:00:00
Date: 2014-07-19 00:00:00
Date: 2014-07-20 00:00:00

Plotting

plt.subplots(figsize=(8,10))
plt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=5)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Date")
plt.grid()
plt.xlim(1,10)
plt.show()

Plotting & zooming around 1-4 Hz

plt.subplots(figsize=(8,10))
plt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=3)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Date")
plt.grid()
plt.xlim(1,4)
plt.show()

Plotting the HVSR curve

plt.subplots(figsize=(8,10))
hvsr = X.mean(axis=1)
plt.plot(hvsr.index, hvsr)
plt.show()


#EOF