03_Obspy_forme_onda.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](images/obspy_logo_full_524x179px.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"font-size: xx-large ; font-weight: 900 ; color: rgba(0 , 0 , 0 , 0.9) ; line-height: 100%\">Breve guida a ObsPy: gestione dei dati della forma d'onda</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Questa guida trae spunto dal Tutorial di:\n",
    "\n",
    "* Lion Krischer ([@krischer](https://github.com/krischer))\n",
    "* Tobias Megies ([@megies](https://github.com/megies))\n",
    "\n",
    "Seismo-Live: http://seismo-live.org\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import print_function\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.figsize'] = 12, 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/Stream_Trace.svg\" width=90%>\n",
    "\n",
    "* la lettura dei dati delle forme d'onda restituisce un oggetto **`Stream`**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy import read\n",
    "st = read(\"./data/waveform_PFO.mseed\", format=\"mseed\")\n",
    "print(st)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Il carattere `*` di UNIX può essere utilizzato per leggere più file contemporaneamente\n",
    "- Rilevamento automatico del formato dei file, non è necessario preoccuparsi dei formati di file\n",
    "\n",
    "   - formati attualmente supportati: **MiniSEED, SAC, SEG-Y, SEG-2, GSE1 / 2, Seisan, SH, DataMark, CSS, wav, y, Q, ASCII, Guralp Compressed Format, Kinemetrics EVT, K-NET e KiK, PDAS, Reftek 130, WIN (continua a crescere ...)**\n",
    "   - sono inclusi nuovi formati ogni volta che viene fornita una routine di lettura di base (o, per esempio, documentazione sufficiente sulla compressione dei dati, ecc.)\n",
    "   - È possibile aggiungere formati di file personalizzati specifici dell'utente (tramite plug-in) per essere rilevati automaticamente in un'installazione ObsPy locale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy import read\n",
    "st = read(\"./data/waveform_*\")\n",
    "print(st)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- per i file MiniSEED, è supportata solo la lettura di brevi parti di file senza che i file vengano letti e tenuti interamente in memoria (consente di risparmiare tempo e memoria quando si lavora su grandi raccolte di file di grandi dimensioni)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy import UTCDateTime\n",
    "\n",
    "t = UTCDateTime(\"2011-03-11T05:46:23.015400Z\")\n",
    "st = read(\"./data/waveform_*\", starttime=t + 10 * 60, endtime=t + 12 * 60)\n",
    "print(st)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### L'oggetto Stream\n",
    "\n",
    " - Un oggetto Stream è una collezione di oggetti Trace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy import read\n",
    "st = read(\"./data/waveform_PFO.mseed\")\n",
    "print(type(st))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(st)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(st.traces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(st[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Più tracce possono essere assemblate usando l'operatore **`+`** (o usando i metodi `.append()` e `.extend()`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st1 = read(\"./data/waveform_PFO.mseed\")\n",
    "st2 = read(\"./data/waveform_PFO_synthetics.mseed\")\n",
    "\n",
    "st = st1 + st2\n",
    "print(st)\n",
    "\n",
    "st3 = read(\"./data/waveform_BFO_BHE.sac\")\n",
    "\n",
    "st += st3\n",
    "print(st)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " - loop  su tracce comodo (e ben leggibile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for tr in st:\n",
    "    print(tr.id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " - Lo stream è utile per applicare la stessa elaborazione a un grande maggiore di forme d'onda diverse o per raggruppare tracce per l'elaborazione (ad esempio tre componenti di una stazione in uno stream)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### L'oggetto Trace (traccia)\n",
    "\n",
    "- un oggetto Trace è un blocco singolo contiguo di dati di forme d'onda (cioè serie temporali a intervalli regolari, senza spazi vuoti)\n",
    "- un oggetto Trace contiene una quantità limitata di metadati in un oggetto simile al dizionario (come **`Trace.stats`**) che descrive completamente le serie temporali specificando..\n",
    "   * posizione / strumento di registrazione (rete, stazione, posizione e codice canale)\n",
    "   * ora di inizio\n",
    "   * frequenza di campionamento"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st = read(\"./data/waveform_PFO.mseed\")\n",
    "tr = st[0]  # get the first Trace in the Stream\n",
    "print(tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tr.stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Per le applicazioni personalizzate a volte è necessario manipolare direttamente i metadati di una Traccia.\n",
    "- I metadati della Traccia rimarranno **coerenti**, poiché tutti i valori sono derivati dall'orario di inizio, i dati e la frequenza di campionamento e sono **aggiornati automaticamente**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tr.stats.delta, \"|\", tr.stats.endtime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tr.stats.sampling_rate = 5.0\n",
    "print(tr.stats.delta, \"|\", tr.stats.endtime)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tr.stats.npts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tr.data = tr.data[:100]\n",
    "print(tr.stats.npts, \"|\", tr.stats.endtime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Alcuni metodi rendono semplici le manipolazioni di base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tr = read(\"./data/waveform_PFO.mseed\")[0]\n",
    "tr.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tr)\n",
    "tr.resample(sampling_rate=100.0)\n",
    "print(tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tr)\n",
    "tr.trim(tr.stats.starttime + 12 * 60, tr.stats.starttime + 14 * 60)\n",
    "print(tr)\n",
    "tr.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tr.detrend(\"linear\")\n",
    "tr.taper(max_percentage=0.05, type='cosine')\n",
    "tr.filter(\"lowpass\", freq=0.1)\n",
    "tr.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# try tr.<Tab> for other methods defined for Trace\n",
    "tr.detrend?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Dati grezzi sono disponibili come [**`numpy.ndarray`**](http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html) (come **`Trace.data`**)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tr.data[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- I dati possono essere modificati direttamente per esempio...\n",
    "\n",
    "...facendo operazioni aritmetiche (veloce, gestito in C da NumPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tr.data ** 2 + 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...usando [**`numpy.ndarray`** metodo predefinito](http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html) (sempre gestito in C da NumPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tr.data.max())\n",
    "print(tr.data.mean())\n",
    "print(tr.data.ptp())\n",
    "# try tr.data.<Tab> per una lista di metodi numpy definiti su ndarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...usando le funzioni **`numpy`** (sempre gestito in C da NumPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "print(np.abs(tr.data))\n",
    "# puoi provare np.<Tab> ma qui c'è molto\n",
    "# try np.a<Tab>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...inserendo puntatori a routine esistenti C/Fortran dall'interno di Python!\n",
    "\n",
    "Questo viene fatto internamente in diversi posti, ad es. per correlazioni incrociate, beamforming o in librerie di tipi di file di terze parti come ad es. libmseed.\n",
    "\n",
    "- Gli oggetti traccia possono anche essere generati manualmente con i dati in un [**`numpy.ndarray`**] (http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html) (ad es. quando è necessario analizzare le forme d'onda da file ascii non standard)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy import Trace\n",
    "\n",
    "x = np.random.randint(-100, 100, 500)\n",
    "tr = Trace(data=x)\n",
    "tr.stats.station = \"XYZ\"\n",
    "tr.stats.starttime = UTCDateTime()\n",
    "\n",
    "tr.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Ogetti Stream possono essere assemblati da Tracce generate manualmente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy import Stream\n",
    "\n",
    "tr2 = Trace(data=np.random.randint(-300, 100, 1000))\n",
    "tr2.stats.starttime = UTCDateTime()\n",
    "tr2.stats.sampling_rate = 10.0\n",
    "st = Stream([tr, tr2])\n",
    "\n",
    "st.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Metodi pre-definiti in  **`Stream`** / **`Trace`**\n",
    "\n",
    "- La maggior parte dei metodi che funziona su un oggetto Trace funziona anche per un oggetto Stream. Sono semplicemente eseguiti su ciascuna traccia. [Vedere la ducumentazione ObsPy per una panoramica dei metodi disponibili](http://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.html) (oppure prova **`st.<Tab>`**).\n",
    " - **`st.filter()`** - Filtra tutte le tracce allegate.\n",
    " - **`st.trim()`** - taglia tutte le tracce.\n",
    " - **`st.resample()`** / **`st.decimate()`** - Cambia il passo di campionamento.\n",
    " - **`st.trigger()`** - Avvia un algoritmo di trigger.\n",
    " - **`st.plot()`** / **`st.spectrogram()`** - Visualizza i dati.\n",
    " - **`st.attach_response()`**/**`st.remove_response()`**, **`st.simulate()`** - Correzione per lo strumento\n",
    " - **`st.merge()`**, **`st.normalize()`**, **`st.detrend()`**, **`st.taper()`**, ...\n",
    "- Un oggetto **`Stream`** può essere anche esportato in molti formati, quindi ObsPy può essere per converire tra diversi formati di file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st = read(\"./data/waveform_*.sac\")\n",
    "st.write(\"output_file.mseed\", format=\"MSEED\")\n",
    "!ls -l output_file*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/Stream_Trace.svg\" width=90%>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Esercizi con le Trace\n",
    " - Crea un **`numpy.ndarray`** con tutti zero (ad esempio usa **`numpy.zeros()`**) e metti un impulso ideale da qualche parte in esso\n",
    "  - Inizializza un oggetto **`Trace`** con il tuo array di dati\n",
    "  - Inserisci alcune informazioni sulla stazione (ad es. Rete, stazione, ..)\n",
    "  - Stampa il riepilogo della traccia e fai un plot la traccia\n",
    "  - Modifica la frequenza di campionamento a 20 Hz\n",
    "  - Modifica l'ora di inizio della traccia all'ora di inizio di questa sessione\n",
    "  - Stampa il riepilogo traccia e fai di nuovo un plot la traccia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.zeros(300)\n",
    "x[100] = 1.0\n",
    "tr = Trace(data=x)\n",
    "tr.stats.station = \"ABC\"\n",
    "tr.plot()\n",
    "tr.stats.sampling_rate = 20.0\n",
    "tr.stats.starttime = UTCDateTime(2018, 4, 13, 15, 0, 0)\n",
    "print(tr)\n",
    "tr.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Usa **`tr.filter(...)`** e applica un filtro passa-basso con una frequanza d'angolo di 1 Hertz.\n",
    "- Visualizza il grafico di anteprima, ci sono alcuni secondi di zeri che possiamo tagliare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tr.filter(\"lowpass\", freq=1)\n",
    "tr.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Usa **`tr.trim(...)`** per rimuovere alcuni degli zero all'inizio ed alla fine\n",
    "- visualizza di nuovo il grafico"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tr.trim(tr.stats.starttime + 3, tr.stats.endtime - 5)\n",
    "tr.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Scala le ampiezze della traccia di un fattore 500\n",
    "- Aggiungi rumore standard gaussiano alla traccia (usa [**`np.random.randn()`**](http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.randn.html))\n",
    "- Visualizza di nuovo la traccia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tr.data = tr.data * 500\n",
    "tr.data = tr.data + np.random.randn(len(tr))\n",
    "tr.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Esercizi con gli Stream\n",
    "\n",
    "- Leggi tutti i dati di terremoto di Tohoku in un oggetto Stream (\"./data/waveform\\_\\*\")\n",
    "- Stampa il riepilogo del flusso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st = read(\"./data/waveform_*\")\n",
    "print(st)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Usa **`st.select()`** per tenere solo le tracce della stazione BFO nello stream. \n",
    "- Visualizza le tracce"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st = st.select(station=\"BFO\")\n",
    "st.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- trim the data to a 10 minute time window around the first arrival (just roughly looking at the preview plot)\n",
    "- display the preview plot and spectrograms for the stream (with logarithmic frequency scale, use `wlen=50` for the spectrogram plot)\n",
    "\n",
    "- Taglia i dati in una finestra temporale di 10 minuti intorno al primo arrivo (guardando solo approssimativamente la figura precedente)\n",
    "- visualizza le tracce e gli spettrogrammi per lo stream (con scala di frequenza logaritmica, usa `wlen=50` per il grafico dello spettrogramma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = UTCDateTime(2011, 3, 11, 5, 55)\n",
    "st.trim(t1, t1 + 10 * 60)\n",
    "st.plot()\n",
    "st.spectrogram(log=True, wlen=50);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- rimuovi l'andamento lineare dai dati, applica un tapering e un filtro passa-basso a 0,1 Hertz\n",
    "- mostra di nuovo il grafico"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st.detrend(\"linear\")\n",
    "st.filter(\"lowpass\", freq=0.1)\n",
    "st.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}