DataAnalysis2022/04-PPSD/1-probabilistic_power_spectral_densities.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style='background-image: url(\"../images/header.svg\") ; padding: 0px ; background-size: cover ; border-radius: 5px ; height: 250px'>\n",
    "    <div style=\"float: right ; margin: 50px ; padding: 20px ; background: rgba(255 , 255 , 255 , 0.7) ; width: 50% ; height: 150px\">\n",
    "        <div style=\"position: relative ; top: 50% ; transform: translatey(-50%)\">\n",
    "            <div style=\"font-size: xx-large ; font-weight: 900 ; color: rgba(0 , 0 , 0 , 0.8) ; line-height: 100%\">Noise</div>\n",
    "            <div style=\"font-size: large ; padding-top: 20px ; color: rgba(0 , 0 , 0 , 0.5)\">Lab: Probabilistic Power Spectral Densities</div>\n",
    "        </div>\n",
    "    </div>\n",
    "</div>\n",
    "\n",
    "\n",
    "Seismo-Live: http://seismo-live.org\n",
    "\n",
    "##### Authors:\n",
    "* Tobias Megies ([@megies](https://github.com/megies))\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use(\"bmh\")\n",
    "plt.rcParams['figure.figsize'] = 10, 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * read waveform data from file `data/GR.FUR..BHN.D.2015.361` (station `FUR`, [LMU geophysical observatory in Fürstenfeldbruck](https://www.geophysik.uni-muenchen.de/observatory/seismology))\n",
    " * read corresponding station metadata from file `data/station_FUR.stationxml`\n",
    " * print info on both waveforms and station metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy import read, read_inventory\n",
    "\n",
    "st = read(\"data/GR.FUR..BHN.D.2015.361\")\n",
    "inv = read_inventory(\"data/station_FUR.stationxml\")\n",
    "\n",
    "print(st)\n",
    "print(inv)\n",
    "inv.plot(projection=\"ortho\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * compute probabilistic power spectral densities using `PPSD` class from obspy.signal, see http://docs.obspy.org/tutorial/code_snippets/probabilistic_power_spectral_density.html (but use the inventory you read from StationXML as metadata)\n",
    " * plot the processed `PPSD` (`plot()` method attached to `PPSD` object)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy.signal import PPSD\n",
    "\n",
    "tr = st[0]\n",
    "ppsd = PPSD(stats=tr.stats, metadata=inv)\n",
    "\n",
    "ppsd.add(tr)\n",
    "ppsd.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since longer term stacks would need too much waveform data and take way too long to compute, we prepared one year continuous data preprocessed for a single channel of station `FUR` to play with..\n",
    "\n",
    " * load long term pre-computed PPSD from file `PPSD_FUR_HHN.npz` using `PPSD`'s `load_npz()` staticmethod (i.e. it is called directly from the class, not an instance object of the class)\n",
    " * plot the PPSD (default is full time-range, depending on how much data and spread is in the data, adjust `max_percentage` option of `plot()` option)  (might take a couple of minutes..!)\n",
    " * do a cumulative plot (which is good to judge non-exceedance percentage dB thresholds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy.signal import PPSD\n",
    "\n",
    "ppsd = PPSD.load_npz(\"data/PPSD_FUR_HHN.npz\", allow_pickle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppsd.plot(max_percentage=10)\n",
    "ppsd.plot(cumulative=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * do different stacks of the data using the [`calculate_histogram()` (see docs!)](http://docs.obspy.org/packages/autogen/obspy.signal.spectral_estimation.PPSD.calculate_histogram.html) method of `PPSD` and visualize them\n",
    " * compare differences in different frequency bands qualitatively (anthropogenic vs. \"natural\" noise)..\n",
    "   * nighttime stack, daytime stack\n",
    " * advanced exercise: Use the `callback` option and use some crazy custom callback function in `calculate_histogram()`, e.g. stack together all data from birthdays in your family.. or all German holidays + Sundays in the time span.. or from dates of some bands' concerts on a tour.. etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppsd.calculate_histogram(time_of_weekday=[(-1, 0, 2), (-1, 22, 24)])\n",
    "ppsd.plot(max_percentage=10)\n",
    "ppsd.calculate_histogram(time_of_weekday=[(-1, 8, 16)])\n",
    "ppsd.plot(max_percentage=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * do different stacks of the data using the [`calculate_histogram()` (see docs!)](http://docs.obspy.org/packages/autogen/obspy.signal.spectral_estimation.PPSD.calculate_histogram.html) method of `PPSD` and visualize them\n",
    " * compare differences in different frequency bands qualitatively (anthropogenic vs. \"natural\" noise)..\n",
    "   * weekdays stack, weekend stack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppsd.calculate_histogram(time_of_weekday=[(1, 0, 24), (2, 0, 24), (3, 0, 24), (4, 0, 24), (5, 0, 24)])\n",
    "ppsd.plot(max_percentage=10)\n",
    "ppsd.calculate_histogram(time_of_weekday=[(6, 0, 24), (7, 0, 24)])\n",
    "ppsd.plot(max_percentage=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * do different stacks of the data using the [`calculate_histogram()` (see docs!)](http://docs.obspy.org/packages/autogen/obspy.signal.spectral_estimation.PPSD.calculate_histogram.html) method of `PPSD` and visualize them\n",
    " * compare differences in different frequency bands qualitatively (anthropogenic vs. \"natural\" noise)..\n",
    "   * seasonal stacks (e.g. northern hemisphere autumn vs. spring/summer, ...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppsd.calculate_histogram(month=[10, 11, 12, 1])\n",
    "ppsd.plot(max_percentage=10)\n",
    "ppsd.calculate_histogram(month=[4, 5, 6, 7])\n",
    "ppsd.plot(max_percentage=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " * do different stacks of the data using the [`calculate_histogram()` (see docs!)](http://docs.obspy.org/packages/autogen/obspy.signal.spectral_estimation.PPSD.calculate_histogram.html) method of `PPSD` and visualize them\n",
    " * compare differences in different frequency bands qualitatively (anthropogenic vs. \"natural\" noise)..\n",
    "   * stacks by specific month\n",
    "   * maybe even combine several of above restrictions.. (e.g. only nighttime on weekends)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "encoding": "# -*- coding: utf-8 -*-"
  },
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}