{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Compute the H/V spectral ratio from the ratios of PSDs\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\n\nif \"SPHINX_DOC_BUILD\" in os.environ:\n    if \"MSNOISE_DOC\" in os.environ:\n        os.chdir(os.environ[\"MSNOISE_DOC\"])\n\nimport matplotlib\n# matplotlib.use(\"agg\")\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nfrom pandas.plotting import register_matplotlib_converters\nimport datetime\n\nregister_matplotlib_converters()\n\nplt.style.use(\"ggplot\")\n\nfrom msnoise.api import *\n\n\n# connect to the database\ndb = connect()\n\n# Obtain a list of dates between ``start_date`` and ``enddate``\nstart, end, datelist = build_movstack_datelist(db)\n\n# Get the list of parameters from the DB:\nparams = get_params(db)\n\n# Get the time axis for plotting the CCF:\ntaxis = get_t_axis(db)\n\n# Get the first mov_stack configured:\nmov_stack = params.mov_stack[0]\n\n\n\ndef convert_to_velocity(df):\n    df = df.resample(\"30Min\").mean()\n    df.columns = 1. / df.columns\n    df = df.sort_index(axis=1)\n    df = df.dropna(axis=1, how='all')\n\n    w2f = (2.0 * np.pi * df.columns)\n    # The acceleration amplitude spectrum and velocity spectral amplitude (not power)\n    vamp = np.sqrt(10.0 ** (df / 10.) / w2f ** 2)\n    return vamp\n\n\n\nZ = hdf_open_store(\"PF.FJS.00.HHZ\", mode=\"r\").PSD\nZ = convert_to_velocity(Z)\nE = hdf_open_store(\"PF.FJS.00.HHE\", mode=\"r\").PSD\nE = convert_to_velocity(E)\nN = hdf_open_store(\"PF.FJS.00.HHN\", mode=\"r\").PSD\nN = convert_to_velocity(N)\n\n\n\nr = \"6h\"\nrZ = Z.resample(r).mean()\nrE = E.resample(r).mean()\nrN = N.resample(r).mean()\n\nmerged = pd.concat({'Z': rZ, 'E': rE, 'N': rN}, axis=0)\nmerged.index.names = ['channel', 'date']\n\n# Swap the levels of the MultiIndex\nresult = merged.swaplevel('channel', 'date').sort_index(level='date')  # .iloc[:500]\nresult\n\n\n\nHVSRs = {}\nfor date, group in result.groupby(level='date'):\n    print(f\"Date: {date}\")\n    group = group.droplevel(0)\n    try:\n        iZ = group.loc[\"Z\"].values\n        iE = group.loc[\"E\"].values\n        iN = group.loc[\"N\"].values\n    except:\n        continue\n    hvsr = np.sqrt((iE + iN) / iZ)\n    HVSRs[date] = hvsr\n\n\n\nHVSR = pd.DataFrame(HVSRs, index=result.columns)\n\n\n\nX = HVSR.copy().fillna(0.0)\nX = X.T.resample(r).mean().T"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.subplots(figsize=(18, 18))\nplt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=4)\nplt.colorbar()\nplt.grid()\nplt.xlabel(\"Frequency (Hz)\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting & zooming around 0.5-4 Hz\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.subplots(figsize=(18, 18))\nplt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=4)\nplt.colorbar()\nplt.xlim(0.5, 4)\nplt.grid()\nplt.xlabel(\"Frequency (Hz)\")\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting the HVSR curve (truncating the low frequency here):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "X.loc[0.2:20].mean(axis=1).plot()\nplt.xlabel(\"Frequency (Hz)\")\nplt.ylabel(\"Amplitude\")\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.10"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}