diff --git a/04-PPSD/1-probabilistic_power_spectral_densities.ipynb b/04-PPSD/1-probabilistic_power_spectral_densities.ipynb new file mode 100644 index 0000000..94635a5 --- /dev/null +++ b/04-PPSD/1-probabilistic_power_spectral_densities.ipynb @@ -0,0 +1,223 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "
\n", + "
Noise
\n", + "
Lab: Probabilistic Power Spectral Densities
\n", + "
\n", + "
\n", + "
\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 +} diff --git a/04-PPSD/data/GR.FUR..BHN.D.2015.361 b/04-PPSD/data/GR.FUR..BHN.D.2015.361 new file mode 100644 index 0000000..124664d Binary files /dev/null and b/04-PPSD/data/GR.FUR..BHN.D.2015.361 differ diff --git a/04-PPSD/data/PPSD_FUR_HHN.npz b/04-PPSD/data/PPSD_FUR_HHN.npz new file mode 100644 index 0000000..646926a Binary files /dev/null and b/04-PPSD/data/PPSD_FUR_HHN.npz differ diff --git a/04-PPSD/data/station_FUR.stationxml b/04-PPSD/data/station_FUR.stationxml new file mode 100644 index 0000000..e353863 --- /dev/null +++ b/04-PPSD/data/station_FUR.stationxml @@ -0,0 +1,232 @@ + + + Jane + Jane + JANE WEB SERVICE: fdsnws-station | Jane version: 0.0.0-gd6e1 + http://jane/fdsnws/station/1/query?network=GR&station=FUR&starttime=2015-12-01&endtime=2016-02-01&level=response + 2016-05-03T10:14:19+00:00 + + GRSN + 1 + 4 + + 48.162899 + 11.2752 + 565.0 + + Fuerstenfeldbruck, Bavaria, GR-Net + + 2006-12-16T00:00:00.000 + 12 + 12 + + 48.162899 + 11.2752 + 565.0 + 0.0 + 0.0 + 0.0 + TRIGGERED + GEOPHYSICAL + 100.0 + 0.02 + + A + Amperes + + + Streckeisen STS-2/N seismometer + + + + 9.4368E8 + 0.02 + + M/S + Velocity in Meters per Second + + + COUNTS + Digital Counts + + + + + + M/S + Velocity in Meters per Second + + + V + Volts + + LAPLACE (RADIANS/SECOND) + 6.0077E7 + 1.0 + + 0.0 + 0.0 + + + 0.0 + 0.0 + + + -0.037004 + 0.037016 + + + -0.037004 + -0.037016 + + + -251.33 + 0.0 + + + -131.04 + -467.29 + + + -131.04 + 467.29 + + + + 1500.0 + 0.02 + + + + + + V + Volts + + + COUNTS + Digital Counts + + DIGITAL + + + 200.0 + 1 + 0 + 0.0 + 0.0 + + + 629121.0 + 0.0 + + + + + + 48.162899 + 11.2752 + 565.0 + 0.0 + 0.0 + 0.0 + TRIGGERED + GEOPHYSICAL + 20.0 + 0.02 + + A + Amperes + + + Streckeisen STS-2/N seismometer + + + + 9.4368E8 + 0.02 + + M/S + Velocity in Meters per Second + + + COUNTS + Digital Counts + + + + + + M/S + Velocity in Meters per Second + + + V + Volts + + LAPLACE (RADIANS/SECOND) + 6.0077E7 + 1.0 + + 0.0 + 0.0 + + + 0.0 + 0.0 + + + -0.037004 + 0.037016 + + + -0.037004 + -0.037016 + + + -251.33 + 0.0 + + + -131.04 + -467.29 + + + -131.04 + 467.29 + + + + 1500.0 + 0.02 + + + + + + V + Volts + + + COUNTS + Digital Counts + + DIGITAL + + + 200.0 + 1 + 0 + 0.0 + 0.0 + + + 629121.0 + 0.0 + + + + + + +