{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"id": "bba6933e",
"metadata": {},
"source": [
"# 7. Seismogram atlas\n",
"\n",
"We have tested PCA and ICA for reducing the dimensions of the scattering coefficients and clustering tasks. We will now investigate Uniform Manifold Projection and Approximation (UMAP) for similar purposes. PCA and ICA preserve quite well the pair-wise distances of data points at the costs of losing information about the local structure. Because they preserve the distances, the techniques are often used for clustering tasks later on. UMAP basically learns the high-dimensional manifold of the data and maps it into a lower dimensional space. By doing so, it preserves much better the local structure but at the costs of global structures, meaning that distances between distant data points might be distorted. Therefore, more care has to be taken if we perform clustering in the UMAP space.\n",
"\n",
"In the literature, these UMAP spaces are called atlases, since they resemble a map of the data such as the activation atlas of a convolutional neural network. Therefore, we decided to call the UMAP spaces of a continuous seismogram a \"seismogram atlas\". Seismogram atlases were introduced in Steinmann et al. 2023."
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "713f93cf",
"metadata": {},
"source": [
"> **Note**: This notebook requires the installation of the matplotlib, sklearn and umap-learn libraries, which are not included in the default installation of the scatseisnet library. To install them, please visit the [matplotlib website](https://matplotlib.org/stable/users/installing.html), the [scikit-learn website](https://scikit-learn.org/stable/install.html) and the [umap-learn website](https://umap-learn.readthedocs.io/en/latest/index.html)."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "cc24ead9",
"metadata": {},
"outputs": [],
"source": [
"import pickle\n",
"\n",
"from matplotlib import pyplot as plt\n",
"import numpy as np\n",
"from sklearn.cluster import DBSCAN, KMeans\n",
"import umap\n",
"import obspy\n",
"\n",
"plt.rcParams[\"date.converter\"] = \"concise\"\n",
"%config InlineBackend.figure_format = \"svg\""
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "db10ba58",
"metadata": {},
"source": [
"## Load scattering coefficients\n",
"\n",
"First, we load the scattering coefficients and reshape them for any dimensionality reduction model (here UMAP) of the `scikit-learn` package. The shape of the scattering coefficients are given by the following tuples:\n",
"\n",
"- Order 1: `(n_times, n_channel, octaves[0] * resolution[0])`\n",
"- Order 2: `(n_times, n_channel, octaves[0] * octaves[1] * resolution[0] * resolution[1])`\n",
"- Order n: `(n_times, n_channel, np.prod(octaves) * np.prod(resolution))`\n",
"\n",
"We then need to collect the all-order scattering coefficients into a two-dimensional matrix for the UMAP application."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "6c652505",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Collected 4319 samples of 504 dimensions each.\n"
]
}
],
"source": [
"# Load data from file\n",
"with np.load(\n",
" \"../example/scattering_coefficients.npz\", allow_pickle=True\n",
") as data:\n",
" order_1 = data[\"order_1\"]\n",
" order_2 = data[\"order_2\"]\n",
" times = data[\"times\"]\n",
"\n",
"# Reshape and stack scattering coefficients of all orders\n",
"order_1 = order_1.reshape(order_1.shape[0], -1)\n",
"order_2 = order_2.reshape(order_2.shape[0], -1)\n",
"scattering_coefficients = np.hstack((order_1, order_2))\n",
"\n",
"# transform into log\n",
"scattering_coefficients = np.log10(scattering_coefficients)\n",
"\n",
"# print info about shape\n",
"n_times, n_coeff = scattering_coefficients.shape\n",
"print(\"Collected {} samples of {} dimensions each.\".format(n_times, n_coeff))"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "15c67bbf",
"metadata": {},
"outputs": [],
"source": [
"# Load the corresponding time vector\n",
"with np.load(\n",
" \"../example/independent_components.npz\", allow_pickle=True\n",
") as data:\n",
" times = data[\"times\"]"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "23596a99",
"metadata": {},
"source": [
"## Apply UMAP\n",
"\n",
"After loading and stacking the scattering coefficients into a matrix, we can now apply UMAP to the scattering coefficient matrix. The `n_components` keyword argument informs the model about how many components we want to extract. Since we want to visualize an atlas, we extract 2 components. Other important hyperparameters are `n_neighbors` and `min_dist`, which control the local versus global structure preservation. We will work with the default parameters."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/seydoux/GitHub/scatseisnet/.venv/lib/python3.13/site-packages/umap/umap_.py:1952: UserWarning: n_jobs value 1 overridden to 1 by setting random_state. Use no seed for parallelism.\n",
" warn(\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" 2025-10-31T17:32:30.163488 \n",
" image/svg+xml \n",
" \n",
" \n",
" Matplotlib v3.10.7, https://matplotlib.org/ \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"umap_model = umap.UMAP(n_components=2, random_state=42)\n",
"embedding = umap_model.fit_transform(scattering_coefficients)\n",
"\n",
"# Visualize the UMAP results\n",
"fig, ax = plt.subplots()\n",
"ax.scatter(embedding[:, 0], embedding[:, 1], s=1)\n",
"ax.set_aspect(\"equal\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Interpret the atlas\n",
"\n",
"The seismogram atlas depicts two distinct clusters, which is very different from the two-dimensional ICA and PCA visualizations in notebook 4. Many cluster algorithms will be able to identify these two distinct structures, so we will use the basic k-means clustering."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" 2025-10-31T17:32:30.459956 \n",
" image/svg+xml \n",
" \n",
" \n",
" Matplotlib v3.10.7, https://matplotlib.org/ \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# clustering task\n",
"N_CLUSTERS = 2\n",
"cluster_model = KMeans(n_clusters=N_CLUSTERS)\n",
"predictions = cluster_model.fit_predict(embedding)\n",
"\n",
"# show clustering results in the atlas\n",
"fig, ax = plt.subplots()\n",
"for label in np.unique(predictions):\n",
" where = label == predictions\n",
" ax.scatter(embedding[where, 0], embedding[where, 1], s=1)\n",
"ax.set_aspect(\"equal\")"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
" 2025-10-31T17:32:30.691442 \n",
" image/svg+xml \n",
" \n",
" \n",
" Matplotlib v3.10.7, https://matplotlib.org/ \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"SMOOTH_KERNEL = 20\n",
"\n",
"# Convert predictions to one-hot encoding\n",
"one_hot = np.zeros((len(times), N_CLUSTERS + 1))\n",
"one_hot[np.arange(len(times)), predictions] = 1\n",
"\n",
"# Plot the results\n",
"fig, ax = plt.subplots(figsize=(6, 6))\n",
"\n",
"# Plot each cluster as a separate line\n",
"for i in range(N_CLUSTERS):\n",
"\n",
" # Obtain the detection rate by convolving with a boxcar kernel\n",
" detection_rate = (\n",
" np.convolve(one_hot[:, i], np.ones(SMOOTH_KERNEL), mode=\"same\")\n",
" / SMOOTH_KERNEL\n",
" )\n",
"\n",
" # Plot the detection rate\n",
" ax.plot(times, one_hot[:, i] + i, alpha=0.5)\n",
" ax.plot(times, detection_rate + i, color=\"black\")\n",
"\n",
"# Labels\n",
"ax.set_xlabel(\"Time\")\n",
"ax.set_ylabel(\"Cluster index\")\n",
"\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "61e9cf9f",
"metadata": {},
"source": [
"## Evaluate clustering results\n",
"\n",
"Our k-means clustering in the seismogram atlas separates very sharply the tremor from the non-tremor period. The results look much cleaner than the k-means clustering in the 4-dimensional feature space of ICA from notebook 6."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
}
],
"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
}