{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"id": "e0929444",
"metadata": {},
"source": [
"# 3. Scattering transformation\n",
"\n",
"If you made up to here, it likely means that you have successfully designed a scattering network for seismic data using ScatSeisNet. In this notebook, we will apply the scattering transformation to the seismograms and visualize the resulting scattering coefficients. This process involves convolving the input signal with the wavelets defined in the scattering network and computing the modulus and pooling operations to obtain the scattering coefficients."
]
},
{
"cell_type": "code",
"execution_count": 101,
"id": "8db7a2e3",
"metadata": {},
"outputs": [],
"source": [
"import pickle\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from obspy import read"
]
},
{
"cell_type": "code",
"execution_count": 102,
"id": "d9c12109",
"metadata": {},
"outputs": [],
"source": [
"%config InlineBackend.figure_format = \"svg\""
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "7d06b446",
"metadata": {},
"source": [
"## Load scattering network and seismograms\n",
"\n",
"First, we need to load the scattering network that we designed in the previous tutorial. We will also load the seismogram data that contains the tremor signals that we introduced at the beginning of this tutorial series."
]
},
{
"cell_type": "code",
"execution_count": 103,
"id": "39c37bc4",
"metadata": {},
"outputs": [],
"source": [
"network = pickle.load(open(\"../example/scattering_network.pickle\", \"rb\"))\n",
"stream = read(\"../example/stream.mseed\")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "c38d374e",
"metadata": {},
"source": [
"## Create seismogram segments\n",
"\n",
"Before applying the scattering network, we need to divide our continuous seismogram into segments of equal length. The segment duration was defined when we created the scattering network in the previous notebook, so we can extract that information directly from the network object. We use overlapping windows to avoid losing information at segment boundaries, with 50% overlap between consecutive segments."
]
},
{
"cell_type": "code",
"execution_count": 104,
"id": "3f4acf7b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Extracted 4319 segments from 2015-12-04 00:00:00 to 2015-12-04 23:59:20\n"
]
}
],
"source": [
"segment_duration_seconds = network.bins / network.sampling_rate\n",
"segment_overlap_seconds = segment_duration_seconds / 2\n",
"\n",
"# Collect segments and timestamps (window start time)\n",
"segments = list()\n",
"timestamps = list()\n",
"for segment in stream.slide(segment_duration_seconds, segment_overlap_seconds):\n",
" timestamps.append(segment[0].stats.starttime.datetime)\n",
" segments.append(np.array([tr.data[: network.bins] for tr in segment]))\n",
"\n",
"# Log\n",
"print(\n",
" f\"Extracted {len(segments)} segments \"\n",
" f\"from {timestamps[0]} to {timestamps[-1]}\"\n",
")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "bb394086",
"metadata": {},
"source": [
"## Apply the scattering transformation\n",
"\n",
"Now we can calculate the scattering coefficients from our data segments. The network processes each segment through its filter banks and applies a pooling operation to summarize the wavelet coefficients. \n",
"\n",
"We use median pooling here, but you could also choose maximum pooling or average pooling by changing the `reduce_type` argument. This computation might take a minute or two depending on your data size."
]
},
{
"cell_type": "code",
"execution_count": 105,
"id": "a3db56cb",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Scattering transformation complete.\n",
"First-order coefficients shape:\n",
" Segments: 4319\n",
" Channels: 3\n",
" Coefficients: 24\n",
"Second-order coefficients shape:\n",
" Segments: 4319\n",
" Channels: 3\n",
" Coefficients: 144\n"
]
}
],
"source": [
"coefficients = network.transform(segments, reduce_type=np.median)\n",
"\n",
"print(f\"Scattering transformation complete.\")\n",
"print(f\"First-order coefficients shape:\")\n",
"print(f\" Segments: {coefficients[0].shape[0]}\")\n",
"print(f\" Channels: {coefficients[0].shape[1]}\")\n",
"print(f\" Coefficients: {coefficients[0].shape[2]}\")\n",
"print(\"Second-order coefficients shape:\")\n",
"print(f\" Segments: {coefficients[1].shape[0]}\")\n",
"print(f\" Channels: {coefficients[1].shape[1]}\")\n",
"print(f\" Coefficients: {np.prod(coefficients[1].shape[2:])}\")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "9ed574be",
"metadata": {},
"source": [
"## Visualize coefficients\n",
"\n",
"Let's examine the first-order scattering coefficients, which resemble a classical spectrogram showing amplitude as a function of time and frequency. In our data, we can clearly see the signature of the tremor signals starting around 10:00 AM, with increasing amplitude throughout the day."
]
},
{
"cell_type": "code",
"execution_count": 106,
"id": "8dd2be40",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Extract first-order coefficients for first channel\n",
"channel = 0\n",
"\n",
"# Get trace\n",
"trace = stream[channel]\n",
"\n",
"# Get scattering coefficients\n",
"center_frequencies_1 = network.banks[0].centers\n",
"center_frequencies_2 = network.banks[1].centers\n",
"order_2_f1 = 1\n",
"order_2_f1_index = np.argmin(center_frequencies_1 - order_2_f1)\n",
"order_1 = np.log10(coefficients[0][:, channel, :].squeeze())\n",
"order_2 = np.log10(coefficients[1][:, channel, order_2_f1_index].squeeze())\n",
"\n",
"\n",
"# Create figure and axes\n",
"fig, ax = plt.subplots(3, sharex=True, dpi=300)\n",
"\n",
"# Plot the waveform\n",
"ax[0].plot(trace.times(\"matplotlib\"), trace.data, rasterized=True, lw=0.5)\n",
"ax[0].grid()\n",
"\n",
"# First-order scattering coefficients\n",
"ax[1].pcolormesh(timestamps, center_frequencies_1, order_1.T, rasterized=True)\n",
"ax[1].grid()\n",
"\n",
"# Second-order scattering coefficients\n",
"ax[2].pcolormesh(timestamps, center_frequencies_2, order_2.T, rasterized=True)\n",
"ax[2].grid()\n",
"\n",
"# Axes labels\n",
"ax[1].set_yscale(\"log\")\n",
"ax[2].set_yscale(\"log\")\n",
"ax[0].set_ylabel(\"Amplitude (counts)\")\n",
"ax[1].set_ylabel(\"$f_1$ (Hz)\")\n",
"ax[2].set_ylabel(\"$f_2$ (Hz)\")\n",
"fig.autofmt_xdate()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "ebaac46b",
"metadata": {},
"source": [
"## Save the scattering coefficients\n",
"\n",
"Now that we've computed the scattering coefficients, we'll save them for use in the next tutorials where we'll explore dimensionality reduction and clustering. We save both the first and second-order coefficients along with the timestamps in a compressed NumPy file."
]
},
{
"cell_type": "code",
"execution_count": 107,
"id": "3d1242ba",
"metadata": {},
"outputs": [],
"source": [
"np.savez(\n",
" \"../example/scattering_coefficients.npz\",\n",
" order_1=coefficients[0],\n",
" order_2=coefficients[1],\n",
" times=timestamps,\n",
")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"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.13.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}