.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples\plot_compute_hvsr_psd.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_plot_compute_hvsr_psd.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_plot_compute_hvsr_psd.py:


Compute the H/V spectral ratio from the ratios of PSDs
======================================================

.. GENERATED FROM PYTHON SOURCE LINES 13-111

.. code-block:: default


    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 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



.. rst-class:: sphx-glr-script-out

.. code-block:: pytb

    Traceback (most recent call last):
      File "D:\PythonForSource\MSNoise_Stack\MSNoise\examples\plot_compute_hvsr_psd.py", line 37, in <module>
        db = connect()
      File "d:\pythonforsource\msnoise_stack\msnoise\msnoise\api.py", line 109, in connect
        engine = get_engine(inifile)
      File "d:\pythonforsource\msnoise_stack\msnoise\msnoise\api.py", line 73, in get_engine
        dbini = read_db_inifile(inifile)
      File "d:\pythonforsource\msnoise_stack\msnoise\msnoise\api.py", line 159, in read_db_inifile
        raise DBConfigNotFoundError(
    msnoise.DBConfigNotFoundError: No db.ini file in this directory, please run 'msnoise db init' in this folder to initialize it as an MSNoise project folder.




.. GENERATED FROM PYTHON SOURCE LINES 112-113

Plotting

.. GENERATED FROM PYTHON SOURCE LINES 113-120

.. code-block:: default

    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()


.. GENERATED FROM PYTHON SOURCE LINES 121-122

Plotting & zooming around 0.5-4 Hz

.. GENERATED FROM PYTHON SOURCE LINES 122-130

.. code-block:: default

    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()


.. GENERATED FROM PYTHON SOURCE LINES 131-132

Plotting the HVSR curve (truncating the low frequency here):

.. GENERATED FROM PYTHON SOURCE LINES 132-139

.. code-block:: default

    X.loc[0.2:20].mean(axis=1).plot()
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Amplitude")
    plt.show()





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  0.005 seconds)


.. _sphx_glr_download_auto_examples_plot_compute_hvsr_psd.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_compute_hvsr_psd.py <plot_compute_hvsr_psd.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_compute_hvsr_psd.ipynb <plot_compute_hvsr_psd.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_