{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Template Matching\n",
"\n",
"Use the events detected with the backprojection method and previously formatted in `6_build_template_database.ipynb` as template events for a matched-filter search."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import BPMF\n",
"import h5py as h5\n",
"import os\n",
"import sys\n",
"\n",
"import beampower as bp\n",
"import glob\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"from BPMF.data_reader_examples import data_reader_mseed\n",
"from obspy import UTCDateTime as udt\n",
"from time import time as give_time\n",
"\n",
"# choose the number of threads you want to limit the computation to\n",
"n_CPUs = 24\n",
"os.environ[\"OMP_NUM_THREADS\"] = str(n_CPUs)\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%config InlineBackend.figure_formats = ['svg']"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# PROGRAM PARAMETERS\n",
"TEMPLATE_DB = \"template_db\"\n",
"MATCHED_FILTER_DB = \"matched_filter_db\"\n",
"PATH_MF = os.path.join(BPMF.cfg.OUTPUT_PATH, MATCHED_FILTER_DB)\n",
"DATA_FOLDER = \"preprocessed_2_12\"\n",
"NETWORK_FILENAME = \"network.csv\"\n",
"\n",
"# TEMPLATE MATCHING PARAMETERS\n",
"# MIN_INTEREVENT_TIME_SEC: minimum inter-event time, in seconds, allowed between 2 detections,\n",
"# keep the highest CC otherwise\n",
"MIN_INTEREVENT_TIME_SEC = 5.0\n",
"# THRESHOLD_WINDOW_DUR_SEC: duration of the sliding window, in seconds, over which the time-dependent\n",
"# detection threshold is computed\n",
"THRESHOLD_WINDOW_DUR_SEC = 1800.0\n",
"# OVERLAP: overlap, between 0 and 1, between sliding windows. This parameter is mainly here for legacy\n",
"# reason but in the later versions of BPMF it does not seem necessary to use overlaping windows.\n",
"OVERLAP = 0.0\n",
"# THRESHOLD_TYPE: 'rms' or 'mad'. Use a time-dependent threshold based on `BPMF.cfg.N_DEV_MF_THRESHOLD` times\n",
"# the rms or the mad. Mad stands for Median absolute deviation and is equal to 1.48xRMS for\n",
"# the normal distribution.\n",
"THRESHOLD_TYPE = \"rms\"\n",
"# WEIGHT_TYPE: 'simple'. BPMF does not support any other weight types for now.\n",
"WEIGHT_TYPE = \"simple\"\n",
"# DEVICE: 'gpu' or 'cpu'. Use 'gpu' if you have Nvidia GPUs and have compiled fast_matched_filter and beampower\n",
"# for you GPUs.\n",
"DEVICE = \"cpu\"\n",
"# MAX_MEMORY_GB: The maximum memory, in gigabytes, that BPMF can use at once. When running a matched-filter search,\n",
"# computing a correlation coefficient time series for each template event can be very memory consuming.\n",
"# Thus, a memory efficient trick is to loop over small batches of the template database. BPMF automatically\n",
"# determines the appropriate number of templates per batch to satisfy this memory constraint.\n",
"MAX_MEMORY_GB = 32.0 # Gb\n",
"# SANITY_CHECK: If true, test the gaussianity of the correlation coefficient time series. A problem-free CC time series\n",
"# is expected to be mostly gaussian distributed. If large deviations are observed, we simply discard the\n",
"# CC time series for the given template.\n",
"SANITY_CHECK = True\n",
"# MAX_KURTOSIS: If SANITY_CHECK is True, the CC time series with kurtosis higher than MAX_KURTOSIS are discarded.\n",
"MAX_KURTOSIS = 100\n",
"# N_SELECTED_STATIONS: For large seismic networks, limit the number of stations over which the matched-filter search\n",
"# is run to the N_SELECTED_STATIONS closest to a given template.\n",
"N_SELECTED_STATIONS = 10\n",
"# MIN_CHANNELS: If the data availability for a given template and a given day drops below MIN_CHANNELS channels, the\n",
"# template is not used for the matched-filter search.\n",
"MIN_CHANNELS = 6\n",
"# MIN_STATIONS: If the data availability for a given template and a given day drops below MIN_STATIONS stations, the\n",
"# template is not used for the matched-filter search.\n",
"MIN_STATIONS = 3\n",
"# NORMALIZE: If True, the continuous waveforms are normalized by their rms before starting the matched-filter search.\n",
"# This avoids issues with numerical imprecisions when dealing with small numbers (note that all the algebra is\n",
"# done in float32 on GPUs and, generally, in BPMF).\n",
"NORMALIZE = True\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Read the template database"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# read the network metadata\n",
"net = BPMF.dataset.Network(NETWORK_FILENAME)\n",
"net.read()\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"template_filenames = glob.glob(\n",
" os.path.join(BPMF.cfg.OUTPUT_PATH, TEMPLATE_DB, \"template*\")\n",
")\n",
"template_filenames.sort()\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.47sec to read the 8 templates.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"ERROR 1: PROJ: proj_create_from_database: Open of /home/ebeauce/miniconda3/envs/hy7_py310/share/proj failed\n"
]
}
],
"source": [
"# initialize the template group\n",
"t_start = give_time()\n",
"template_group = BPMF.dataset.TemplateGroup.read_from_files(template_filenames, net)\n",
"# read and normalize waveforms\n",
"template_group.read_waveforms()\n",
"template_group.normalize(method=\"rms\")\n",
"t_end = give_time()\n",
"print(f\"{t_end-t_start:.2f}sec to read the {template_group.n_templates}\" \" templates.\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Initialize the `BPMF.similarity_search.MatchedFilter` instance."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# initialize the MatchedFilter object\n",
"MF = BPMF.similarity_search.MatchedFilter(\n",
" template_group,\n",
" normalize=NORMALIZE,\n",
" threshold_type=THRESHOLD_TYPE,\n",
" max_memory=MAX_MEMORY_GB,\n",
" max_workers=n_CPUs,\n",
" max_kurto=MAX_KURTOSIS,\n",
" min_channels=MIN_CHANNELS,\n",
" min_stations=MIN_STATIONS,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load the data"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Days to process: DatetimeIndex(['2012-07-07 00:00:00+00:00'], dtype='datetime64[ns, UTC]', freq='D')\n"
]
}
],
"source": [
"date_list = net.datelist()\n",
"print(\"Days to process: \", date_list)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Processing 2012-07-07\n",
"24 Trace(s) in Stream:\n",
"YH.SAUV..HHZ | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DE07..BHE | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DC07..BHE | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DD06..BHZ | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DC07..BHN | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DE07..BHN | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.SPNC..BHZ | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DE08..BHZ | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DC08..BHZ | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DC06..BHZ | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DC08..BHE | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.SPNC..BHE | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DE08..BHE | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.SAUV..HHN | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DC06..BHE | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DD06..BHN | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DC07..BHZ | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DE07..BHZ | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DD06..BHE | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.SPNC..BHN | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DE08..BHN | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DC08..BHN | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.SAUV..HHE | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n",
"YH.DC06..BHN | 2012-07-07T00:00:00.000000Z - 2012-07-08T00:00:00.000000Z | 25.0 Hz, 2160001 samples\n"
]
}
],
"source": [
"# this could be written as a loop:\n",
"#for date in date_list:\n",
"date = date_list[0]\n",
"date = pd.Timestamp(date)\n",
"date_str = date.strftime(\"%Y-%m-%d\")\n",
"data_root_folder = os.path.join(BPMF.cfg.INPUT_PATH, str(date.year), date.strftime(\"%Y%m%d\"))\n",
"t_start_day = give_time()\n",
"print(f\"Processing {date_str}\")\n",
"# -------------------------------------------\n",
"# Load the data\n",
"# -------------------------------------------\n",
"data = BPMF.dataset.Data(\n",
" date,\n",
" data_root_folder,\n",
" data_reader_mseed,\n",
" duration=24.0 * 3600.0,\n",
" sampling_rate=BPMF.cfg.SAMPLING_RATE_HZ,\n",
" db_path=BPMF.cfg.INPUT_PATH,\n",
")\n",
"data.read_waveforms(data_folder=DATA_FOLDER)\n",
"print(data.traces.__str__(extended=True))\n",
"# get the data in a numpy.ndarray\n",
"data_arr = data.get_np_array(net.stations)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# compute data availability\n",
"data.set_availability(net.stations)\n",
"# re-arrange selected stations on each template\n",
"template_group.n_closest_stations(\n",
" N_SELECTED_STATIONS,\n",
" available_stations=data.availability.index[data.availability]\n",
" )\n",
"# attach the `Data` instance to the `MatchedFilter` instance\n",
"MF.set_data(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Run the matched-filter search"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"27 correlation computations were skipped on the 0-th template. Can be caused by zeros in data, or too low amplitudes (try to increase the gain).\n",
"Template 0 detected 30 events.\n",
"Template 1 detected 29 events.\n",
"Template 2 detected 33 events.\n",
"Template 3 detected 17 events.\n",
"Template 4 detected 32 events.\n",
"Template 5 detected 1 events.\n",
"Template 6 detected 15 events.\n",
"Template 7 detected 1 events.\n",
"Total time spent on computing CCs: 1.18sec\n",
"Total time spent on finding detections: 0.40sec\n"
]
}
],
"source": [
"# we use verbose = 1 for this example but you might want to use verbose = 0 for large-scale applications\n",
"detections = MF.run_matched_filter_search(\n",
" MIN_INTEREVENT_TIME_SEC,\n",
" weight_type=WEIGHT_TYPE,\n",
" device=DEVICE,\n",
" threshold_window_dur=THRESHOLD_WINDOW_DUR_SEC,\n",
" overlap=OVERLAP,\n",
" sanity_check=SANITY_CHECK,\n",
" verbose=1,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot the detected events"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# choose the template you want to plot\n",
"TEMPLATE_ID = 2"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": "\n\n\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig = MF.plot_cc(TEMPLATE_ID, figsize=(15, 5))"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.48sec to read the waveforms for all 8 detections.\n"
]
}
],
"source": [
"# read the detected events' waveforms for plotting purposes\n",
"# note: on a large-scale application, we generally just want to save the detected events' metadata to the output database\n",
"\n",
"# event extraction parameters\n",
"\n",
"# PHASE_ON_COMP: dictionary defining which moveout we use to extract the waveform\n",
"PHASE_ON_COMP = {\"N\": \"S\", \"1\": \"S\", \"E\": \"S\", \"2\": \"S\", \"Z\": \"P\"}\n",
"# OFFSET_PHASE: dictionary defining the time offset taken before a given phase\n",
"# for example OFFSET_PHASE[\"P\"] = 1.0 means that we extract the window\n",
"# 1 second before the predicted P arrival time\n",
"OFFSET_PHASE = {\"P\": 1.0, \"S\": 4.0}\n",
"# TIME_SHIFTED: boolean, if True, use moveouts to extract relevant windows\n",
"TIME_SHIFTED = True\n",
"\n",
"t1 = give_time()\n",
"for i in range(len(detections)):\n",
" detections[TEMPLATE_ID][i].read_waveforms(\n",
" BPMF.cfg.TEMPLATE_LEN_SEC,\n",
" phase_on_comp=PHASE_ON_COMP,\n",
" offset_phase=OFFSET_PHASE,\n",
" time_shifted=TIME_SHIFTED,\n",
" data_folder=DATA_FOLDER,\n",
" )\n",
"t2 = give_time()\n",
"print(\n",
" f\"{t2-t1:.2f}sec to read the waveforms for all {len(detections):d} detections.\"\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using `BPMF.similarity.MatchedFilter.plot_detection` method, we plot the waveforms of a detected event (black) on top of which we plot the template waveforms (red). The location in time of the detected event is also shown by an arrow on the CC time series (top plot)."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": "\n\n\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# choose the event idx you want to plot\n",
"EVENT_IDX = 0\n",
"fig = MF.plot_detection(\n",
" detections[TEMPLATE_ID][EVENT_IDX],\n",
" figsize=(15, 20),\n",
")\n"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": "\n\n\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# choose the event idx you want to plot\n",
"EVENT_IDX = 4\n",
"fig = MF.plot_detection(\n",
" detections[TEMPLATE_ID][EVENT_IDX],\n",
" figsize=(15, 20),\n",
")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Save the detected events to database\n",
"\n",
"Here, we save the detected events' metadata in a hdf5 file that is unique to each template. Each detected event is stored in a hdf5 `Group` which is identified by a string with its origin time."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing took 1.81sec.\n"
]
}
],
"source": [
"# save detections in their parent-template file\n",
"if not os.path.isdir(PATH_MF):\n",
" os.mkdir(PATH_MF)\n",
"n_total = 0\n",
"t_writing_start = give_time()\n",
"# serial version\n",
"for tid in detections.keys():\n",
" parent_tp_filename = f\"detections_template{tid}.h5\"\n",
" fullpath = os.path.join(PATH_MF, parent_tp_filename)\n",
" # output_file = h5.File(fullpath, mode='a')\n",
" for i in range(len(detections[tid])):\n",
" tid = detections[tid][i].aux_data[\"tid\"]\n",
" gid = detections[tid][i].origin_time.strftime(\"%Y-%m-%dT%H:%M:%S\")\n",
" detections[tid][i].write(\n",
" parent_tp_filename, db_path=PATH_MF, gid=gid\n",
" ) # , hdf5_file=output_file)\n",
" n_total += 1\n",
"t_writing_end = give_time()\n",
"print(f\"Writing took {t_writing_end-t_writing_start:.2f}sec.\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are almost done! Note that neighboring templates are very likely to detect some events in common. Thus, to build a meaningful earthquake catalog with no (or, more realistically, few) multiple detections of the same events, we still need to clean up the database (see next notebook)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.10.4 ('hy7_py310')",
"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.4"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "221f0e5b1b98151b07a79bf3b6d0c1d306576197d2c4531763770570a29e708e"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}