09_Anak_Krakatau_eruption_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/seismo-live/seismo_live/refs/heads/master/notebooks/Workshops/2017_Baku_STCU_IRIS_ObsPy_course/images/obspy_logo_full_524x179px.png\" width=90%> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebooks we will look at the Anak Krakatau eruption from 2018 in more detail. We will use the data that you downloaded in the FDSN exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "  import obspy\n",
    "  import cartopy\n",
    "except ModuleNotFoundError:\n",
    "  !pip -qq install obspy cartopy\n",
    "  import cartopy\n",
    "  import obspy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://github.com/yannikbehr/BMKG_ObsPy_workshop/raw/refs/heads/main/data/indonesia_waveforms.tar.gz\n",
    "!tar -xzf indonesia_waveforms.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start again by first loading the catalogue around the eruption time at 2018-12-22T13:55:00. We will use Geofon's catalogue as IRIS' catalogue doesn't contain this event."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy import read, read_events, UTCDateTime\n",
    "from obspy.clients.fdsn import Client\n",
    "\n",
    "t = UTCDateTime(\"2018-12-22T13:55:00\")\n",
    "client = Client(\"Geofon\")\n",
    "catalog = client.get_events(starttime=t-180, endtime=t + 1800,\n",
    "                            minmagnitude=2)\n",
    "print(catalog)\n",
    "catalog.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are interested in the M5.1 event. Select the event from the catalog and print its latitude, longitude and depth using the `preferred_origin()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event = catalog[-1]\n",
    "event_lat = event.preferred_origin().latitude\n",
    "event_lon = event.preferred_origin().longitude\n",
    "event_depth = event.preferred_origin().depth\n",
    "print(event_lat, event_lon, event_depth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now read in the waveform data `data/anak_krakatau_2018.mseed` into a `stream`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st = read('data/anak_krakatau_2018.mseed')\n",
    "print(st)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we have to download the instrument response for each station which we will need for the magnitude computation later on. To do so use the `level=response` keyword argument to [get_stations](https://docs.obspy.org/packages/autogen/obspy.clients.fdsn.client.Client.get_stations.html#obspy.clients.fdsn.client.Client.get_stations)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy import Inventory\n",
    "inv = Inventory()\n",
    "for tr in st:\n",
    "    inv_ = client.get_stations(network=tr.stats.network, station=tr.stats.station,\n",
    "                               channel=tr.stats.channel, location=tr.stats.location,\n",
    "                               starttime=tr.stats.starttime, endtime=tr.stats.endtime,\n",
    "                               level=\"response\")\n",
    "    inv += inv_\n",
    "print(inv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to plot a record section of the data as shown [here](https://docs.obspy.org/tutorial/code_snippets/waveform_plotting_tutorial.html#plotting-a-record-section). To do so we first need to compute the distance of each station to the event in meters. You can use the function [gps2dist_azimuth](https://docs.obspy.org/packages/autogen/obspy.geodetics.base.gps2dist_azimuth.html#obspy.geodetics.base.gps2dist_azimuth) from the [obspy.geodetics](https://docs.obspy.org/packages/obspy.geodetics.html#module-obspy.geodetics) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy.geodetics import gps2dist_azimuth\n",
    "\n",
    "for tr in st:\n",
    "    station_lat = inv.select(station=tr.stats.station).networks[0].stations[0].latitude\n",
    "    station_lon = inv.select(station=tr.stats.station).networks[0].stations[0].longitude\n",
    "    dist, az, baz = gps2dist_azimuth(event_lat, event_lon,\n",
    "                                     station_lat, station_lon)\n",
    "    tr.stats.distance = dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let's calculate some rough theoretical travel times for the P-wave using the ak135 earth model. You can do that using the [obspy.taup](https://docs.obspy.org/packages/obspy.taup.html#module-obspy.taup) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy.taup import TauPyModel\n",
    "model = TauPyModel(model=\"ak135\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy.geodetics import kilometer2degrees\n",
    "print(event.preferred_origin().time)\n",
    "for tr in st:\n",
    "    distance = kilometer2degrees(tr.stats.distance/1000)\n",
    "    arrivals = model.get_ray_paths(source_depth_in_km=0,\n",
    "                                    distance_in_degree=distance,\n",
    "                                    phase_list=[\"P\"])\n",
    "    tr.stats.p_arrival = event.preferred_origin().time + arrivals[0].time\n",
    "    print(tr.id, event.preferred_origin().time + arrivals[0].time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may have noticed that you get five different travel times for the P-wave. To find out why plot the ray paths using the [get_ray_paths](https://docs.obspy.org/packages/obspy.taup.html#module-obspy.taup) and [plot_rays](https://docs.obspy.org/packages/obspy.taup.html#module-obspy.taup) methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ray_paths = model.get_ray_paths(0, distance, phase_list=[\"P\"])\n",
    "ray_paths.plot_rays(plot_type=\"cartesian\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you have everything to plot a record section for the data. For an extra challenge try to add station names and theoretical first arrival times. You can get some ideas how to do that by looking at the source code of the [example](https://docs.obspy.org/tutorial/code_snippets/waveform_plotting_tutorial.html#plotting-a-record-section)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.transforms import blended_transform_factory\n",
    "\n",
    "fig = plt.figure()\n",
    "st.detrend(\"demean\")\n",
    "st.plot(type='section', time_down=True, fig=fig, norm_method='stream')\n",
    "ax = fig.axes[0]\n",
    "transform = blended_transform_factory(ax.transData, ax.transAxes)\n",
    "for tr in st:\n",
    "    ax.text(tr.stats.distance / 1e3, 1.0, tr.stats.station, rotation=270,\n",
    "            va=\"bottom\", ha=\"center\", transform=transform, zorder=10)\n",
    "    ax.plot([tr.stats.distance/1e3 - 10, tr.stats.distance/1e3 + 10],\n",
    "            [tr.stats.p_arrival - tr.stats.starttime, tr.stats.p_arrival - tr.stats.starttime],\n",
    "            color='red', lw=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate the magnitude first plot the [spectrograms](https://docs.obspy.org/packages/obspy.imaging.html#module-obspy.imaging) so you can get an idea of the frequency content of the surface waves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st.spectrogram(wlen=100);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you have all the information to compute the magnitude of the event. We will use [estimate_magnitude](https://docs.obspy.org/packages/autogen/obspy.signal.invsim.estimate_magnitude.html#obspy.signal.invsim.estimate_magnitude) to do so. Can you think of reasons why the magnitude was underestimated?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from obspy.signal.invsim import estimate_magnitude\n",
    "\n",
    "for station in ['BBJI', 'PMBI']:\n",
    "    response = inv.get_response(seed_id=f\"GE.{station}..LHZ\", datetime=event.preferred_origin().time)\n",
    "    tr = st.select(station=station)[0]\n",
    "    max_amplitude = tr.data.max()\n",
    "    magnitude = estimate_magnitude(response, 2*abs(max_amplitude), 1.5, tr.stats.distance/1e3)\n",
    "    print(magnitude)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pyvolprob",
   "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}