example_local.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1c5e48a2-922e-4fff-a800-78b2ab11bfac",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 0. Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc7f3d3a-b46c-4252-973e-73be4e1dd687",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "from urllib.request import urlretrieve\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "from obspy import UTCDateTime\n",
    "from waveform_collection import gather_waveforms\n",
    "\n",
    "from rtm import (calculate_time_buffer, define_grid, get_peak_coordinates,\n",
    "                 grid_search, plot_record_section, plot_st, plot_time_slice,\n",
    "                 process_waveforms, produce_dem)\n",
    "\n",
    "# Ignore benign Matplotlib backend warning due to fig.show()\n",
    "warnings.filterwarnings(action='ignore', message='Matplotlib is currently using module')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49efcb76-3580-46ca-b645-96a9ead2000b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1. Define grid\n",
    "\n",
    "To obtain `DEM_WGS84.tif` from OpenTopography (~200 MB GeoTIFF), run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cf71665",
   "metadata": {},
   "outputs": [],
   "source": [
    "DEM_URL = 'https://opentopography.s3.sdsc.edu/dataspace/OTDS.072019.4326.1/raster/DEM_WGS84.tif'\n",
    "dem_file = DEM_URL.split('/')[-1]\n",
    "if not Path(dem_file).is_file():\n",
    "    print('Downloading DEM...')\n",
    "    urlretrieve(DEM_URL, dem_file)\n",
    "    print('Done')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "567ec2a4-5ac5-43d6-b504-56adf9658426",
   "metadata": {},
   "source": [
    "Alternatively, specify `EXTERNAL_FILE = None` below to automatically download and use 1 arc-second SRTM data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afee6300-e5a1-4ec1-9aa6-c60f81f1573a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "EXTERNAL_FILE = dem_file\n",
    "\n",
    "LON_0 = 169.448212  # [deg] Longitude of grid center\n",
    "LAT_0 = -19.527908  # [deg] Latitude of grid center\n",
    "\n",
    "X_RADIUS = 600  # [m] E-W grid radius (half of grid \"width\")\n",
    "Y_RADIUS = 650  # [m] N-S grid radius (half of grid \"height\")\n",
    "SPACING = 10  # [m] Grid spacing\n",
    "\n",
    "grid = define_grid(\n",
    "    lon_0=LON_0,\n",
    "    lat_0=LAT_0,\n",
    "    x_radius=X_RADIUS,\n",
    "    y_radius=Y_RADIUS,\n",
    "    spacing=SPACING,\n",
    "    projected=True,\n",
    ")\n",
    "\n",
    "dem = produce_dem(grid, external_file=EXTERNAL_FILE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03175d31-110a-4eb2-b888-3269c5adfcf2",
   "metadata": {},
   "source": [
    "## 2. Grab and process the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26af6c3b-d05a-4fef-978f-15b14d51ee4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start and end of time window containing (suspected) events\n",
    "STARTTIME = UTCDateTime('2016-07-30T05:22:45')\n",
    "ENDTIME = STARTTIME + 10\n",
    "\n",
    "# Data collection parameters\n",
    "SOURCE = 'IRIS'\n",
    "NETWORK = '3E'\n",
    "STATION = 'YIF?'\n",
    "LOCATION = '*'\n",
    "CHANNEL = '*'\n",
    "\n",
    "MAX_STATION_DIST = 0.8  # [km] Max. dist. from grid center to station (approx.)\n",
    "\n",
    "FREQ_MIN = 0.5  # [Hz] Lower bandpass corner\n",
    "FREQ_MAX = 10  # [Hz] Upper bandpass corner\n",
    "\n",
    "DECIMATION_RATE = 10  # [Hz] New sampling rate to use for decimation\n",
    "SMOOTH_WIN = 1  # [s] Smoothing window duration\n",
    "\n",
    "# Automatically determine appropriate time buffer in s\n",
    "time_buffer = calculate_time_buffer(grid, MAX_STATION_DIST)\n",
    "\n",
    "st = gather_waveforms(\n",
    "    source=SOURCE,\n",
    "    network=NETWORK,\n",
    "    station=STATION,\n",
    "    location=LOCATION,\n",
    "    channel=CHANNEL,\n",
    "    starttime=STARTTIME,\n",
    "    endtime=ENDTIME,\n",
    "    time_buffer=time_buffer,\n",
    ")\n",
    "\n",
    "st_proc = process_waveforms(\n",
    "    st,\n",
    "    freqmin=FREQ_MIN,\n",
    "    freqmax=FREQ_MAX,\n",
    "    envelope=True,\n",
    "    smooth_win=SMOOTH_WIN,\n",
    "    decimation_rate=DECIMATION_RATE,\n",
    "    normalize=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16dc4bc8-cfb6-4b50-9c36-1d5abfca419b",
   "metadata": {},
   "source": [
    "## 3. Perform grid search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae85d524-8a0f-4850-a6d2-c87a98a5116f",
   "metadata": {},
   "outputs": [],
   "source": [
    "STACK_METHOD = 'sum'  # Choose either 'sum', 'product', or 'semblance'\n",
    "TIME_METHOD = 'celerity'  # Choose either 'celerity' or 'fdtd'\n",
    "TIME_KWARGS = {'celerity': 343, 'dem': dem}\n",
    "\n",
    "S = grid_search(\n",
    "    processed_st=st_proc,\n",
    "    grid=grid,\n",
    "    time_method=TIME_METHOD,\n",
    "    starttime=STARTTIME,\n",
    "    endtime=ENDTIME,\n",
    "    stack_method=STACK_METHOD,\n",
    "    **TIME_KWARGS\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7a5681c-fd16-48b9-8f02-265e2d128757",
   "metadata": {},
   "source": [
    "## 4. Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3702cc84-3b4b-4c5d-8f6f-c15811c8598d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig_st = plot_st(\n",
    "    st,\n",
    "    filt=[FREQ_MIN, FREQ_MAX],\n",
    "    equal_scale=False,\n",
    "    remove_response=True,\n",
    "    label_waveforms=True,\n",
    ")\n",
    "\n",
    "fig_slice = plot_time_slice(\n",
    "    S, st_proc, label_stations=True, dem=dem, plot_peak=True, xy_grid=600\n",
    ")\n",
    "\n",
    "time_max, y_max, x_max, peaks, props = get_peak_coordinates(\n",
    "    S, global_max=False, height=3, min_time=2, unproject=S.UTM\n",
    ")\n",
    "\n",
    "peak_ind = np.argmax(props['peak_heights'])\n",
    "\n",
    "fig = plot_record_section(\n",
    "    st_proc,\n",
    "    origin_time=time_max[peak_ind],\n",
    "    source_location=(y_max[peak_ind], x_max[peak_ind]),\n",
    "    plot_celerity=S.celerity,\n",
    "    label_waveforms=True,\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}