.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples\plot_compute_hvsr.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.py>`
        to download the full example code

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

.. _sphx_glr_auto_examples_plot_compute_hvsr.py:


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

.. GENERATED FROM PYTHON SOURCE LINES 9-10

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

.. GENERATED FROM PYTHON SOURCE LINES 10-51

.. 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, 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")



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

.. code-block:: pytb

    Traceback (most recent call last):
      File "D:\PythonForSource\MSNoise_Stack\MSNoise\examples\plot_compute_hvsr.py", line 32, 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 52-53

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

.. GENERATED FROM PYTHON SOURCE LINES 53-66

.. code-block:: default



    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]


.. GENERATED FROM PYTHON SOURCE LINES 67-68

Helper functions

.. GENERATED FROM PYTHON SOURCE LINES 68-77

.. code-block:: default



    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


.. GENERATED FROM PYTHON SOURCE LINES 78-79

Doing the job !

.. GENERATED FROM PYTHON SOURCE LINES 79-99

.. code-block:: default



    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



.. GENERATED FROM PYTHON SOURCE LINES 100-101

Plotting

.. GENERATED FROM PYTHON SOURCE LINES 101-111

.. code-block:: default


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



.. GENERATED FROM PYTHON SOURCE LINES 112-113

Plotting & zooming around 1-4 Hz

.. GENERATED FROM PYTHON SOURCE LINES 113-122

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 123-124

Plotting the HVSR curve

.. GENERATED FROM PYTHON SOURCE LINES 124-131

.. code-block:: default


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


    #EOF

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

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


.. _sphx_glr_download_auto_examples_plot_compute_hvsr.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.py <plot_compute_hvsr.py>`

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

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


.. only:: html

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

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