Compute the H/V spectral ratio from the imaginary part of the CCFs
Import & getting the data from the computed autocorrelations (ZZ, EE, NN)
[1]:
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
plt.style.use("ggplot")
from msnoise.core.db import connect
from msnoise.results import MSNoiseResult
# connect to the database
db = connect()
# Build a result object at the stack step (filter 2 for HVSR)
result = MSNoiseResult.from_ids(db, preprocess=1, cc=1, filter=2, stack=1)
params = result.params
# Get the first mov_stack configured:
mov_stack = params.stack.mov_stack[0]
# Get the autocorrelation CCFs (ZZ, EE, NN) for station PF.FJS.00
ZZ = result.get_ccf(pair="PF.FJS.00:PF.FJS.00", components="ZZ",
mov_stack=mov_stack)
EE = result.get_ccf(pair="PF.FJS.00:PF.FJS.00", components="EE",
mov_stack=mov_stack)
NN = result.get_ccf(pair="PF.FJS.00:PF.FJS.00", components="NN",
mov_stack=mov_stack)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File D:\PythonForSource\MSNoise_Stack\MSNoise\msnoise\core\db.py:153, in read_db_inifile(inifile)
152 try:
--> 153 f = open(inifile, 'rb')
154 #except FileNotFoundError: # This is better but only for python3
FileNotFoundError: [Errno 2] No such file or directory: 'D:\\PythonForSource\\MSNoise_Stack\\MSNoise\\doc\\auto_examples\\db.ini'
During handling of the above exception, another exception occurred:
DBConfigNotFoundError Traceback (most recent call last)
Cell In[1], line 20
16 from msnoise.results import MSNoiseResult
19 # connect to the database
---> 20 db = connect()
22 # Build a result object at the stack step (filter 2 for HVSR)
23 result = MSNoiseResult.from_ids(db, preprocess=1, cc=1, filter=2, stack=1)
File D:\PythonForSource\MSNoise_Stack\MSNoise\msnoise\core\db.py:104, in connect(inifile)
101 if not inifile:
102 inifile = os.path.join(os.getcwd(), 'db.ini')
--> 104 engine = get_engine(inifile)
105 session = sessionmaker(bind=engine)
106 return session()
File D:\PythonForSource\MSNoise_Stack\MSNoise\msnoise\core\db.py:67, in get_engine(inifile)
65 from sqlalchemy import create_engine
66 from sqlalchemy.pool import NullPool
---> 67 dbini = read_db_inifile(inifile)
68 if dbini.tech == 1:
69 engine = create_engine('sqlite:///%s' % dbini.hostname, echo=False,
70 connect_args={'check_same_thread': False})
File D:\PythonForSource\MSNoise_Stack\MSNoise\msnoise\core\db.py:156, in read_db_inifile(inifile)
154 #except FileNotFoundError: # This is better but only for python3
155 except IOError:
--> 156 raise DBConfigNotFoundError(
157 "No db.ini file in this directory, please run "
158 "'msnoise db init' in this folder to initialize it as "
159 "an MSNoise project folder.")
160 try:
161 # New ini file with prefix support
162 tech, hostname, database, username, password, prefix = pickle.load(f)
DBConfigNotFoundError: No db.ini file in this directory, please run 'msnoise db init' in this folder to initialize it as an MSNoise project folder.
Checking the data has the same time base (index) and joining
[2]:
r = "1D"
rZZ = MSNoiseResult.to_dataframe(ZZ).resample(r).mean()
rEE = MSNoiseResult.to_dataframe(EE).resample(r).mean()
rNN = MSNoiseResult.to_dataframe(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_df = merged.swaplevel("channel", "date").sort_index(level="date")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 2
1 r = "1D"
----> 2 rZZ = MSNoiseResult.to_dataframe(ZZ).resample(r).mean()
3 rEE = MSNoiseResult.to_dataframe(EE).resample(r).mean()
4 rNN = MSNoiseResult.to_dataframe(NN).resample(r).mean()
NameError: name 'ZZ' is not defined
Helper functions
[3]:
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!
[4]:
HVSRs = {}
for date, group in result_df.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 Exception:
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
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 2
1 HVSRs = {}
----> 2 for date, group in result_df.groupby(level="date"):
3 print(f"Date: {date}")
4 group = group.droplevel(0)
NameError: name 'result_df' is not defined
Plotting
[5]:
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()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 2
1 plt.subplots(figsize=(8, 10))
----> 2 plt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=5)
3 plt.xlabel("Frequency (Hz)")
4 plt.ylabel("Date")
NameError: name 'X' is not defined
Plotting & zooming around 1-4 Hz
[6]:
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()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 2
1 plt.subplots(figsize=(8, 10))
----> 2 plt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=3)
3 plt.xlabel("Frequency (Hz)")
4 plt.ylabel("Date")
NameError: name 'X' is not defined
Plotting the HVSR curve
[7]:
plt.subplots(figsize=(8, 10))
hvsr = X.mean(axis=1)
plt.plot(hvsr.index, hvsr)
plt.show()
#EOF
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 2
1 plt.subplots(figsize=(8, 10))
----> 2 hvsr = X.mean(axis=1)
3 plt.plot(hvsr.index, hvsr)
4 plt.show()
NameError: name 'X' is not defined