{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Compute the H/V spectral ratio from the imaginary part of the CCFs\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import & getting the data from the computed autocorrelations (ZZ, EE, NN)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\nif \"SPHINX_DOC_BUILD\" in os.environ:\n    if \"MSNOISE_DOC\" in os.environ:\n        os.chdir(os.environ[\"MSNOISE_DOC\"])\n\nimport matplotlib\nmatplotlib.use(\"agg\")\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\n\nplt.style.use(\"ggplot\")\n\nfrom msnoise.core.db import connect\nfrom msnoise.results import MSNoiseResult\n\n\n# connect to the database\ndb = connect()\n\n# Build a result object at the stack step (filter 2 for HVSR)\nresult = MSNoiseResult.from_ids(db, preprocess=1, cc=1, filter=2, stack=1)\nparams = result.params\n\n# Get the first mov_stack configured:\nmov_stack = params.stack.mov_stack[0]\n\n# Get the autocorrelation CCFs (ZZ, EE, NN) for station PF.FJS.00\nZZ = result.get_ccf(pair=\"PF.FJS.00:PF.FJS.00\", components=\"ZZ\",\n                    mov_stack=mov_stack)\nEE = result.get_ccf(pair=\"PF.FJS.00:PF.FJS.00\", components=\"EE\",\n                    mov_stack=mov_stack)\nNN = result.get_ccf(pair=\"PF.FJS.00:PF.FJS.00\", components=\"NN\",\n                    mov_stack=mov_stack)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Checking the data has the same time base (index) and joining\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "r = \"1D\"\nrZZ = MSNoiseResult.to_dataframe(ZZ).resample(r).mean()\nrEE = MSNoiseResult.to_dataframe(EE).resample(r).mean()\nrNN = MSNoiseResult.to_dataframe(NN).resample(r).mean()\n\nmerged = pd.concat({\"ZZ\": rZZ, \"EE\": rEE, \"NN\": rNN}, axis=0)\nmerged.index.names = [\"channel\", \"date\"]\n\n# Swap the levels of the MultiIndex\nresult_df = merged.swaplevel(\"channel\", \"date\").sort_index(level=\"date\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Helper functions\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def get_imag(d, fs):\n    NFFT = 1024 * 32\n    iX = np.imag(np.fft.fft(d, n=NFFT)[:NFFT // 2])\n    freqs = np.fft.fftfreq(NFFT, d=1 / fs)[:NFFT // 2]\n    return freqs, iX"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Doing the job!\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "HVSRs = {}\nfor date, group in result_df.groupby(level=\"date\"):\n    print(f\"Date: {date}\")\n    group = group.droplevel(0)\n    try:\n        f, iZZ = get_imag(group.loc[\"ZZ\"].values, 20)\n        f, iEE = get_imag(group.loc[\"EE\"].values, 20)\n        f, iNN = get_imag(group.loc[\"NN\"].values, 20)\n    except Exception:\n        continue\n    hvsr = np.sqrt((iEE + iNN) / iZZ)\n    HVSRs[date] = hvsr\n\nHVSR = pd.DataFrame(HVSRs, index=f)\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=(8, 10))\nplt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=5)\nplt.xlabel(\"Frequency (Hz)\")\nplt.ylabel(\"Date\")\nplt.grid()\nplt.xlim(1, 10)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting & zooming around 1-4 Hz\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.subplots(figsize=(8, 10))\nplt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=3)\nplt.xlabel(\"Frequency (Hz)\")\nplt.ylabel(\"Date\")\nplt.grid()\nplt.xlim(1, 4)\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting the HVSR curve\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.subplots(figsize=(8, 10))\nhvsr = X.mean(axis=1)\nplt.plot(hvsr.index, hvsr)\nplt.show()\n\n\n#EOF"
      ]
    }
  ],
  "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.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}