{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n# Compute the H/V spectral ratio from the imaginary part of the CCFs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "source": [
    "import os\n\nif \"SPHINX_DOC_BUILD\" in os.environ or 1:\n    os.chdir(r\"C:\\tmp\\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\n\n\nplt.subplots(figsize=(18, 18))\nplt.pcolormesh(X.index, X.columns, X.T, rasterized=True, vmax=4)\nplt.colorbar()\nplt.grid()\nplt.xlabel(\"Frequency (Hz)\")\n\n\nplt.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)\")\n\n\nX.loc[0.2:20].mean(axis=1).plot()\nplt.xlabel(\"Frequency (Hz)\")\nplt.ylabel(\"Amplitude\")"
   ],
   "outputs": []
  }
 ],
 "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}