adding examples

doc/conf.py

@@ -145,7 +145,7 @@ def process_command_or_group(out, commands, chain=[]):
      'examples_dirs': '../examples',   # path to your example scripts
      'gallery_dirs': 'auto_examples',  # path where to save gallery generated examples
      'show_memory': False,
-     'thumbnail_size': (250, 250),
+     'thumbnail_size': (400, 400),
 }
 intersphinx_cache_limit = 5
 

examples/README.rst

@@ -4,7 +4,7 @@ Interaction Examples & Gallery
 The following examples are meant to show you how to interact with MSNoise using
 its API, thus avoiding having to dive into the folder structure.
 
-Users should try examples while checking the :doc:`../api`. (application
+Users should try examples while checking the :doc:`../api` (application
 programming interface) for understanding the calls to different functions.
 
 In a nutshell, all examples start with the following Python code:

examples/plot_compute_hvsr.py

@@ -0,0 +1,130 @@
+# -*- coding: utf-8 -*-
+"""
+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 or 1:
+    os.chdir(r"C:\tmp\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
+
+
+##############################################################
+# 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

examples/plot_compute_hvsr_psd.py

@@ -0,0 +1,137 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[4]:
+
+
+# -*- coding: utf-8 -*-
+"""
+Compute the H/V spectral ratio from the ratios of PSDs
+======================================================
+
+"""
+
+import os
+
+if "SPHINX_DOC_BUILD" in os.environ or 1:
+    os.chdir(r"C:\tmp\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 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]
+
+
+
+def convert_to_velocity(df):
+    df = df.resample("30Min").mean()
+    df.columns = 1. / df.columns
+    df = df.sort_index(axis=1)
+    df = df.dropna(axis=1, how='all')
+
+    w2f = (2.0 * np.pi * df.columns)
+    # The acceleration amplitude spectrum and velocity spectral amplitude (not power)
+    vamp = np.sqrt(10.0 ** (df / 10.) / w2f ** 2)
+    return vamp
+
+
+
+Z = hdf_open_store("PF.FJS.00.HHZ", mode="r").PSD
+Z = convert_to_velocity(Z)
+E = hdf_open_store("PF.FJS.00.HHE", mode="r").PSD
+E = convert_to_velocity(E)
+N = hdf_open_store("PF.FJS.00.HHN", mode="r").PSD
+N = convert_to_velocity(N)
+
+
+
+r = "6h"
+rZ = Z.resample(r).mean()
+rE = E.resample(r).mean()
+rN = N.resample(r).mean()
+
+merged = pd.concat({'Z': rZ, 'E': rE, 'N': rN}, axis=0)
+merged.index.names = ['channel', 'date']
+
+# Swap the levels of the MultiIndex
+result = merged.swaplevel('channel', 'date').sort_index(level='date')  # .iloc[:500]
+result
+
+
+
+HVSRs = {}
+for date, group in result.groupby(level='date'):
+    print(f"Date: {date}")
+    group = group.droplevel(0)
+    try:
+        iZ = group.loc["Z"].values
+        iE = group.loc["E"].values
+        iN = group.loc["N"].values
+    except:
+        continue
+    hvsr = np.sqrt((iE + iN) / iZ)
+    HVSRs[date] = hvsr
+
+
+
+HVSR = pd.DataFrame(HVSRs, index=result.columns)
+
+
+
+X = HVSR.copy().fillna(0.0)
+X = X.T.resample(r).mean().T
+
+##############################################################
+# Plotting
+plt.subplots(figsize=(18, 18))
+plt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=4)
+plt.colorbar()
+plt.grid()
+plt.xlabel("Frequency (Hz)")
+plt.show()
+
+##############################################################
+# Plotting & zooming around 0.5-4 Hz
+plt.subplots(figsize=(18, 18))
+plt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=4)
+plt.colorbar()
+plt.xlim(0.5, 4)
+plt.grid()
+plt.xlabel("Frequency (Hz)")
+plt.show()
+
+##############################################################
+# Plotting the HVSR curve (truncating the low frequency here):
+X.loc[0.2:20].mean(axis=1).plot()
+plt.xlabel("Frequency (Hz)")
+plt.ylabel("Amplitude")
+plt.show()
+
+
+

msnoise/s04_stack2.py

@@ -49,7 +49,7 @@
 
 * |ref_begin|
 * |ref_end|
-* |mov_stack|
+* |mov_stack| | * changed in 2.0 !!
 * |stack_method| | *new in 1.4*
 * |pws_timegate| | *new in 1.4*
 * |pws_power| | *new in 1.4*