initial import of PPSD exercise from seismo live web site
This commit is contained in:
parent
d829006103
commit
dbcb4a3d8a
223
04-PPSD/1-probabilistic_power_spectral_densities.ipynb
Normal file
223
04-PPSD/1-probabilistic_power_spectral_densities.ipynb
Normal file
@ -0,0 +1,223 @@
|
||||
{
|
||||
"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
|
||||
}
|
||||
BIN
04-PPSD/data/GR.FUR..BHN.D.2015.361
Normal file
BIN
04-PPSD/data/GR.FUR..BHN.D.2015.361
Normal file
Binary file not shown.
BIN
04-PPSD/data/PPSD_FUR_HHN.npz
Normal file
BIN
04-PPSD/data/PPSD_FUR_HHN.npz
Normal file
Binary file not shown.
232
04-PPSD/data/station_FUR.stationxml
Normal file
232
04-PPSD/data/station_FUR.stationxml
Normal file
@ -0,0 +1,232 @@
|
||||
<?xml version='1.0' encoding='UTF-8'?>
|
||||
<FDSNStationXML xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.0" xsi:schemaLocation="http://www.fdsn.org/xml/station/1 http://www.fdsn.org/xml/station/fdsn-station-1.0.xsd">
|
||||
<Source>Jane</Source>
|
||||
<Sender>Jane</Sender>
|
||||
<Module>JANE WEB SERVICE: fdsnws-station | Jane version: 0.0.0-gd6e1</Module>
|
||||
<ModuleURI>http://jane/fdsnws/station/1/query?network=GR&station=FUR&starttime=2015-12-01&endtime=2016-02-01&level=response</ModuleURI>
|
||||
<Created>2016-05-03T10:14:19+00:00</Created>
|
||||
<Network startDate="2006-12-16T00:00:00.000000Z" code="GR">
|
||||
<Description>GRSN</Description>
|
||||
<SelectedNumberStations>1</SelectedNumberStations>
|
||||
<TotalNumberStations>4</TotalNumberStations>
|
||||
<Station startDate="2006-12-16T00:00:00.000000Z" code="FUR">
|
||||
<Latitude plusError="0.0" minusError="0.0">48.162899</Latitude>
|
||||
<Longitude plusError="0.0" minusError="0.0">11.2752</Longitude>
|
||||
<Elevation plusError="0.0" minusError="0.0">565.0</Elevation>
|
||||
<Site>
|
||||
<Name>Fuerstenfeldbruck, Bavaria, GR-Net</Name>
|
||||
</Site>
|
||||
<CreationDate>2006-12-16T00:00:00.000</CreationDate>
|
||||
<SelectedNumberChannels>12</SelectedNumberChannels>
|
||||
<TotalNumberChannels>12</TotalNumberChannels>
|
||||
<Channel locationCode="" code="HHN" startDate="2006-12-16T00:00:00.000">
|
||||
<Latitude plusError="0.0" minusError="0.0">48.162899</Latitude>
|
||||
<Longitude plusError="0.0" minusError="0.0">11.2752</Longitude>
|
||||
<Elevation plusError="0.0" minusError="0.0">565.0</Elevation>
|
||||
<Depth plusError="0.0" minusError="0.0">0.0</Depth>
|
||||
<Azimuth plusError="0.0" minusError="0.0">0.0</Azimuth>
|
||||
<Dip plusError="0.0" minusError="0.0">0.0</Dip>
|
||||
<Type>TRIGGERED</Type>
|
||||
<Type>GEOPHYSICAL</Type>
|
||||
<SampleRate plusError="0.0" minusError="0.0">100.0</SampleRate>
|
||||
<ClockDrift plusError="0.0" minusError="0.0">0.02</ClockDrift>
|
||||
<CalibrationUnits>
|
||||
<Name>A</Name>
|
||||
<Description>Amperes</Description>
|
||||
</CalibrationUnits>
|
||||
<Sensor>
|
||||
<Type>Streckeisen STS-2/N seismometer</Type>
|
||||
</Sensor>
|
||||
<Response>
|
||||
<InstrumentSensitivity>
|
||||
<Value>9.4368E8</Value>
|
||||
<Frequency>0.02</Frequency>
|
||||
<InputUnits>
|
||||
<Name>M/S</Name>
|
||||
<Description>Velocity in Meters per Second</Description>
|
||||
</InputUnits>
|
||||
<OutputUnits>
|
||||
<Name>COUNTS</Name>
|
||||
<Description>Digital Counts</Description>
|
||||
</OutputUnits>
|
||||
</InstrumentSensitivity>
|
||||
<Stage number="1">
|
||||
<PolesZeros>
|
||||
<InputUnits>
|
||||
<Name>M/S</Name>
|
||||
<Description>Velocity in Meters per Second</Description>
|
||||
</InputUnits>
|
||||
<OutputUnits>
|
||||
<Name>V</Name>
|
||||
<Description>Volts</Description>
|
||||
</OutputUnits>
|
||||
<PzTransferFunctionType>LAPLACE (RADIANS/SECOND)</PzTransferFunctionType>
|
||||
<NormalizationFactor>6.0077E7</NormalizationFactor>
|
||||
<NormalizationFrequency plusError="0.0" minusError="0.0">1.0</NormalizationFrequency>
|
||||
<Zero number="0">
|
||||
<Real plusError="0.0" minusError="0.0">0.0</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">0.0</Imaginary>
|
||||
</Zero>
|
||||
<Zero number="1">
|
||||
<Real plusError="0.0" minusError="0.0">0.0</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">0.0</Imaginary>
|
||||
</Zero>
|
||||
<Pole number="0">
|
||||
<Real plusError="0.0" minusError="0.0">-0.037004</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">0.037016</Imaginary>
|
||||
</Pole>
|
||||
<Pole number="1">
|
||||
<Real plusError="0.0" minusError="0.0">-0.037004</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">-0.037016</Imaginary>
|
||||
</Pole>
|
||||
<Pole number="2">
|
||||
<Real plusError="0.0" minusError="0.0">-251.33</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">0.0</Imaginary>
|
||||
</Pole>
|
||||
<Pole number="3">
|
||||
<Real plusError="0.0" minusError="0.0">-131.04</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">-467.29</Imaginary>
|
||||
</Pole>
|
||||
<Pole number="4">
|
||||
<Real plusError="0.0" minusError="0.0">-131.04</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">467.29</Imaginary>
|
||||
</Pole>
|
||||
</PolesZeros>
|
||||
<StageGain>
|
||||
<Value>1500.0</Value>
|
||||
<Frequency>0.02</Frequency>
|
||||
</StageGain>
|
||||
</Stage>
|
||||
<Stage number="2">
|
||||
<Coefficients>
|
||||
<InputUnits>
|
||||
<Name>V</Name>
|
||||
<Description>Volts</Description>
|
||||
</InputUnits>
|
||||
<OutputUnits>
|
||||
<Name>COUNTS</Name>
|
||||
<Description>Digital Counts</Description>
|
||||
</OutputUnits>
|
||||
<CfTransferFunctionType>DIGITAL</CfTransferFunctionType>
|
||||
</Coefficients>
|
||||
<Decimation>
|
||||
<InputSampleRate plusError="0.0" minusError="0.0">200.0</InputSampleRate>
|
||||
<Factor>1</Factor>
|
||||
<Offset>0</Offset>
|
||||
<Delay plusError="0.0" minusError="0.0">0.0</Delay>
|
||||
<Correction plusError="0.0" minusError="0.0">0.0</Correction>
|
||||
</Decimation>
|
||||
<StageGain>
|
||||
<Value>629121.0</Value>
|
||||
<Frequency>0.0</Frequency>
|
||||
</StageGain>
|
||||
</Stage>
|
||||
</Response>
|
||||
</Channel>
|
||||
<Channel locationCode="" code="BHN" startDate="2006-12-16T00:00:00.000">
|
||||
<Latitude plusError="0.0" minusError="0.0">48.162899</Latitude>
|
||||
<Longitude plusError="0.0" minusError="0.0">11.2752</Longitude>
|
||||
<Elevation plusError="0.0" minusError="0.0">565.0</Elevation>
|
||||
<Depth plusError="0.0" minusError="0.0">0.0</Depth>
|
||||
<Azimuth plusError="0.0" minusError="0.0">0.0</Azimuth>
|
||||
<Dip plusError="0.0" minusError="0.0">0.0</Dip>
|
||||
<Type>TRIGGERED</Type>
|
||||
<Type>GEOPHYSICAL</Type>
|
||||
<SampleRate plusError="0.0" minusError="0.0">20.0</SampleRate>
|
||||
<ClockDrift plusError="0.0" minusError="0.0">0.02</ClockDrift>
|
||||
<CalibrationUnits>
|
||||
<Name>A</Name>
|
||||
<Description>Amperes</Description>
|
||||
</CalibrationUnits>
|
||||
<Sensor>
|
||||
<Type>Streckeisen STS-2/N seismometer</Type>
|
||||
</Sensor>
|
||||
<Response>
|
||||
<InstrumentSensitivity>
|
||||
<Value>9.4368E8</Value>
|
||||
<Frequency>0.02</Frequency>
|
||||
<InputUnits>
|
||||
<Name>M/S</Name>
|
||||
<Description>Velocity in Meters per Second</Description>
|
||||
</InputUnits>
|
||||
<OutputUnits>
|
||||
<Name>COUNTS</Name>
|
||||
<Description>Digital Counts</Description>
|
||||
</OutputUnits>
|
||||
</InstrumentSensitivity>
|
||||
<Stage number="1">
|
||||
<PolesZeros>
|
||||
<InputUnits>
|
||||
<Name>M/S</Name>
|
||||
<Description>Velocity in Meters per Second</Description>
|
||||
</InputUnits>
|
||||
<OutputUnits>
|
||||
<Name>V</Name>
|
||||
<Description>Volts</Description>
|
||||
</OutputUnits>
|
||||
<PzTransferFunctionType>LAPLACE (RADIANS/SECOND)</PzTransferFunctionType>
|
||||
<NormalizationFactor>6.0077E7</NormalizationFactor>
|
||||
<NormalizationFrequency plusError="0.0" minusError="0.0">1.0</NormalizationFrequency>
|
||||
<Zero number="0">
|
||||
<Real plusError="0.0" minusError="0.0">0.0</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">0.0</Imaginary>
|
||||
</Zero>
|
||||
<Zero number="1">
|
||||
<Real plusError="0.0" minusError="0.0">0.0</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">0.0</Imaginary>
|
||||
</Zero>
|
||||
<Pole number="0">
|
||||
<Real plusError="0.0" minusError="0.0">-0.037004</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">0.037016</Imaginary>
|
||||
</Pole>
|
||||
<Pole number="1">
|
||||
<Real plusError="0.0" minusError="0.0">-0.037004</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">-0.037016</Imaginary>
|
||||
</Pole>
|
||||
<Pole number="2">
|
||||
<Real plusError="0.0" minusError="0.0">-251.33</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">0.0</Imaginary>
|
||||
</Pole>
|
||||
<Pole number="3">
|
||||
<Real plusError="0.0" minusError="0.0">-131.04</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">-467.29</Imaginary>
|
||||
</Pole>
|
||||
<Pole number="4">
|
||||
<Real plusError="0.0" minusError="0.0">-131.04</Real>
|
||||
<Imaginary plusError="0.0" minusError="0.0">467.29</Imaginary>
|
||||
</Pole>
|
||||
</PolesZeros>
|
||||
<StageGain>
|
||||
<Value>1500.0</Value>
|
||||
<Frequency>0.02</Frequency>
|
||||
</StageGain>
|
||||
</Stage>
|
||||
<Stage number="2">
|
||||
<Coefficients>
|
||||
<InputUnits>
|
||||
<Name>V</Name>
|
||||
<Description>Volts</Description>
|
||||
</InputUnits>
|
||||
<OutputUnits>
|
||||
<Name>COUNTS</Name>
|
||||
<Description>Digital Counts</Description>
|
||||
</OutputUnits>
|
||||
<CfTransferFunctionType>DIGITAL</CfTransferFunctionType>
|
||||
</Coefficients>
|
||||
<Decimation>
|
||||
<InputSampleRate plusError="0.0" minusError="0.0">200.0</InputSampleRate>
|
||||
<Factor>1</Factor>
|
||||
<Offset>0</Offset>
|
||||
<Delay plusError="0.0" minusError="0.0">0.0</Delay>
|
||||
<Correction plusError="0.0" minusError="0.0">0.0</Correction>
|
||||
</Decimation>
|
||||
<StageGain>
|
||||
<Value>629121.0</Value>
|
||||
<Frequency>0.0</Frequency>
|
||||
</StageGain>
|
||||
</Stage>
|
||||
</Response>
|
||||
</Channel>
|
||||
</Station>
|
||||
</Network>
|
||||
</FDSNStationXML>
|
||||
Loading…
Reference in New Issue
Block a user