diff --git a/Fig4_PyGMT_pandas.ipynb b/Fig4_PyGMT_pandas.ipynb
new file mode 100644
index 00000000..0dde8c37
--- /dev/null
+++ b/Fig4_PyGMT_pandas.ipynb
@@ -0,0 +1,94 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "533bc728-db2c-4f3a-ab3f-1f4f2c9ecd69",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import io\n",
+    "\n",
+    "import pandas as pd\n",
+    "import pygmt\n",
+    "import requests\n",
+    "from pygmt.params import Box\n",
+    "\n",
+    "params = {\n",
+    "    \"format\": \"csv\",\n",
+    "    \"starttime\": \"2000-01-01\",\n",
+    "    \"endtime\": \"2025-10-30\",\n",
+    "    \"mindepth\": 70,\n",
+    "    \"minmagnitude\": 5,\n",
+    "}\n",
+    "r = requests.get(\"https://earthquake.usgs.gov/fdsnws/event/1/query\", params=params)\n",
+    "df_eqs = pd.read_csv(io.StringIO(r.text))\n",
+    "df_eqs = df_eqs[\n",
+    "    (df_eqs[\"longitude\"] >= 91)\n",
+    "    & (df_eqs[\"longitude\"] <= 134)\n",
+    "    & (df_eqs[\"latitude\"] >= -22)\n",
+    "    & (df_eqs[\"latitude\"] <= 18)\n",
+    "].sort_values(by=\"mag\", ascending=False)\n",
+    "\n",
+    "fig = pygmt.Figure()\n",
+    "fig.basemap(region=[91, 134, -22, 18], projection=\"M15c\", frame=True)\n",
+    "fig.coast(land=\"gray95\", shorelines=\"gray50\")\n",
+    "\n",
+    "# Plot epicenters with color (hypocentral depth) or size (magnitude)\n",
+    "pygmt.makecpt(cmap=\"SCM/navia\", series=[0, 700], reverse=True, transparency=30)\n",
+    "fig.plot(\n",
+    "    x=df_eqs.longitude,\n",
+    "    y=df_eqs.latitude,\n",
+    "    style=\"c\",\n",
+    "    size=0.005 * 2**df_eqs.mag,\n",
+    "    fill=df_eqs.depth,\n",
+    "    cmap=True,\n",
+    "    pen=\"gray10\",\n",
+    ")\n",
+    "fig.colorbar(frame=[\"xaf+lHypocentral depth\", \"y+lkm\"])\n",
+    "# Add legend for size-coding\n",
+    "legend = io.StringIO(\n",
+    "    \"\\n\".join(f\"S 0.4 c {0.005 * 2**mag:.2f} - 1p 1 M {mag}\" for mag in [5, 6, 7])\n",
+    ")\n",
+    "fig.legend(spec=legend, position=\"jBR+o0.2c+l2\", box=Box(fill=\"white\", pen=\"black\"))\n",
+    "\n",
+    "# Add histogram for magnitude distribution\n",
+    "with fig.inset(position=\"jBL+w7c/4c+o0.2c\", margin=(1.3, 0.2, 1.1, 0.2), box=True):\n",
+    "    fig.histogram(\n",
+    "        region=[4.8, 9.2, 0, 0],\n",
+    "        projection=\"X?/?\",\n",
+    "        frame=[\"WSrt\", \"xa1f0.2+lMagnitude\", \"yaf+lCounts\"],\n",
+    "        data=df_eqs.mag,\n",
+    "        series=0.2,\n",
+    "        fill=\"darkgray\",\n",
+    "        pen=\"lightgray\",\n",
+    "        histtype=0,\n",
+    "    )\n",
+    "\n",
+    "fig.show()\n",
+    "# fig.savefig(\"Fig4_PyGMT_pandas.png\")"
+   ]
+  }
+ ],
+ "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.12.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}