{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Miscellaneous Utilities\n", "The follow demonstrates some miscellaneous utilities included in ObsPlus.\n", "\n", "## Geodetics\n", "It is often necessary to calculate geometric parameters (distance, azimuth, etc.) for pairs of entities in two different groups. For example, distance from each event in a catalog to each receiver in an inventory. ObsPlus provides a simple class for completing this task." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:36.768046Z", "iopub.status.busy": "2025-01-09T18:05:36.767868Z", "iopub.status.idle": "2025-01-09T18:05:38.682551Z", "shell.execute_reply": "2025-01-09T18:05:38.681866Z" } }, "outputs": [], "source": [ "# Load the catalog and inventory from the crandall dataset\n", "import obsplus\n", "from obsplus.utils.geodetics import SpatialCalculator\n", "\n", "crandall = obsplus.load_dataset(\"crandall_test\")\n", "cat = crandall.event_client.get_events()\n", "inv = crandall.station_client.get_stations()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.684754Z", "iopub.status.busy": "2025-01-09T18:05:38.684564Z", "iopub.status.idle": "2025-01-09T18:05:38.687180Z", "shell.execute_reply": "2025-01-09T18:05:38.686725Z" } }, "outputs": [], "source": [ "# init a SpatialCalculator instance (defaults to Earth's params)\n", "spatial_calc = SpatialCalculator()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.688924Z", "iopub.status.busy": "2025-01-09T18:05:38.688597Z", "iopub.status.idle": "2025-01-09T18:05:38.829577Z", "shell.execute_reply": "2025-01-09T18:05:38.828970Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/runner/micromamba/envs/test/lib/python3.10/site-packages/obspy/taup/taup_geo.py:105: UserWarning: Assuming spherical planet when calculating epicentral distance. Install the Python module 'geographiclib' to solve this.\n", " warnings.warn(msg)\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
distance_mazimuthback_azimuthdistance_degreesvertical_distance_m
id1id2
smi:local/248828UU.SRU..HHE72882.074942302.534256122.0821660.6547118374.0
smi:local/248843UU.SRU..HHE72425.955953302.739971122.2917440.6506133854.0
smi:local/248839UU.SRU..HHE72254.546800303.034657122.5889540.6490742214.0
smi:local/248925UU.SRU..HHE71648.952510303.260251122.8194280.6436336424.0
smi:local/248891UU.SRU..HHE71821.459138302.940209122.4967370.6451835044.0
\n", "
" ], "text/plain": [ " distance_m azimuth back_azimuth \\\n", "id1 id2 \n", "smi:local/248828 UU.SRU..HHE 72882.074942 302.534256 122.082166 \n", "smi:local/248843 UU.SRU..HHE 72425.955953 302.739971 122.291744 \n", "smi:local/248839 UU.SRU..HHE 72254.546800 303.034657 122.588954 \n", "smi:local/248925 UU.SRU..HHE 71648.952510 303.260251 122.819428 \n", "smi:local/248891 UU.SRU..HHE 71821.459138 302.940209 122.496737 \n", "\n", " distance_degrees vertical_distance_m \n", "id1 id2 \n", "smi:local/248828 UU.SRU..HHE 0.654711 8374.0 \n", "smi:local/248843 UU.SRU..HHE 0.650613 3854.0 \n", "smi:local/248839 UU.SRU..HHE 0.649074 2214.0 \n", "smi:local/248925 UU.SRU..HHE 0.643633 6424.0 \n", "smi:local/248891 UU.SRU..HHE 0.645183 5044.0 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# create distance dataframe\n", "df = spatial_calc(entity_1=cat, entity_2=inv)\n", "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since a multi-index is used it provides a fairly intuitive way to look up particular event-channel pairs using a tuple of (event_id, seed_id) in conjunction with the `.loc` DataFrame attribute like so:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.857031Z", "iopub.status.busy": "2025-01-09T18:05:38.856690Z", "iopub.status.idle": "2025-01-09T18:05:38.859671Z", "shell.execute_reply": "2025-01-09T18:05:38.859150Z" } }, "outputs": [], "source": [ "event_id = str(cat[0].resource_id)\n", "seed_id = \"UU.MPU..HHZ\"\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.861528Z", "iopub.status.busy": "2025-01-09T18:05:38.861162Z", "iopub.status.idle": "2025-01-09T18:05:38.863578Z", "shell.execute_reply": "2025-01-09T18:05:38.863099Z" } }, "outputs": [], "source": [ "# or just get a particular parameter" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The distances can be converted to km and the distrtibutions can be described." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.865377Z", "iopub.status.busy": "2025-01-09T18:05:38.865025Z", "iopub.status.idle": "2025-01-09T18:05:38.880766Z", "shell.execute_reply": "2025-01-09T18:05:38.880299Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
azimuthback_azimuthdistance_degreesdistance_kmvertical_distance_km
count408.00408.00408.00408.00408.00
mean201.06211.690.8291.245.22
std104.3798.060.3336.551.85
min9.1648.810.1618.261.70
25%126.47122.450.6370.463.69
50%163.87189.790.8290.855.50
75%293.44307.271.12124.366.25
max358.12344.771.29144.009.31
\n", "
" ], "text/plain": [ " azimuth back_azimuth distance_degrees distance_km \\\n", "count 408.00 408.00 408.00 408.00 \n", "mean 201.06 211.69 0.82 91.24 \n", "std 104.37 98.06 0.33 36.55 \n", "min 9.16 48.81 0.16 18.26 \n", "25% 126.47 122.45 0.63 70.46 \n", "50% 163.87 189.79 0.82 90.85 \n", "75% 293.44 307.27 1.12 124.36 \n", "max 358.12 344.77 1.29 144.00 \n", "\n", " vertical_distance_km \n", "count 408.00 \n", "mean 5.22 \n", "std 1.85 \n", "min 1.70 \n", "25% 3.69 \n", "50% 5.50 \n", "75% 6.25 \n", "max 9.31 " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Convert add km columns and delete m columns\n", "m_columns = [x for x in df.columns if x.endswith(\"_m\")]\n", "km_columns = [x.replace(\"_m\", \"_km\") for x in m_columns]\n", "\n", "df_km = df.assign(**{x: df[y] / 1000.0 for x, y in zip(km_columns, m_columns)}).drop(\n", " columns=m_columns\n", ")\n", "\n", "# Calculate stats for source reseiver distances\n", "df_km.describe().round(decimals=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Time\n", "\n", "Working with dates/times can be frustrating, especially since (as of 2020) ObsPy, numpy, and pandas all use slightly different methods for working with time. ObsPlus provides some utilities to make things a little easier." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.882664Z", "iopub.status.busy": "2025-01-09T18:05:38.882305Z", "iopub.status.idle": "2025-01-09T18:05:38.885199Z", "shell.execute_reply": "2025-01-09T18:05:38.884637Z" } }, "outputs": [], "source": [ "\n", "import numpy as np\n", "import obspy\n", "\n", "import obsplus\n", "from obsplus.utils import to_datetime64, to_timedelta64, to_utc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note:** All ObsPlus datafames use numpy/pandas datatypes. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.886970Z", "iopub.status.busy": "2025-01-09T18:05:38.886640Z", "iopub.status.idle": "2025-01-09T18:05:38.910611Z", "shell.execute_reply": "2025-01-09T18:05:38.909992Z" } }, "outputs": [], "source": [ "df = obsplus.events_to_df(obspy.read_events())" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.912587Z", "iopub.status.busy": "2025-01-09T18:05:38.912135Z", "iopub.status.idle": "2025-01-09T18:05:38.916373Z", "shell.execute_reply": "2025-01-09T18:05:38.915811Z" } }, "outputs": [ { "data": { "text/plain": [ "0 2012-04-04 14:21:42.300\n", "1 2012-04-04 14:18:37.000\n", "2 2012-04-04 14:08:46.000\n", "Name: time, dtype: datetime64[ns]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df[\"time\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numpy and ObsPy time differences" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One difference between numpy's `datatime64` and ObsPy's `UTCDateTime` is how offsets are applied. For ObsPy, numbers are simply taken as seconds, but numpy requires explicitly using `timedeta64` instances." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.918327Z", "iopub.status.busy": "2025-01-09T18:05:38.917982Z", "iopub.status.idle": "2025-01-09T18:05:38.920922Z", "shell.execute_reply": "2025-01-09T18:05:38.920411Z" } }, "outputs": [], "source": [ "time_str = \"2020-01-03T11:00:00\"\n", "utc = obspy.UTCDateTime(time_str)\n", "dt64 = np.datetime64(time_str, \"ns\")" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.922683Z", "iopub.status.busy": "2025-01-09T18:05:38.922327Z", "iopub.status.idle": "2025-01-09T18:05:38.924768Z", "shell.execute_reply": "2025-01-09T18:05:38.924280Z" } }, "outputs": [], "source": [ "# add 1/2 second to UTCDateTime\n", "utc2 = utc + 0.5" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.926628Z", "iopub.status.busy": "2025-01-09T18:05:38.926279Z", "iopub.status.idle": "2025-01-09T18:05:38.929006Z", "shell.execute_reply": "2025-01-09T18:05:38.928535Z" } }, "outputs": [], "source": [ "# however doing the same thing with datetime64 raises a TypeError\n", "try:\n", " dt64 + 0.5\n", "except TypeError:\n", " pass" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.930583Z", "iopub.status.busy": "2025-01-09T18:05:38.930401Z", "iopub.status.idle": "2025-01-09T18:05:38.933045Z", "shell.execute_reply": "2025-01-09T18:05:38.932567Z" } }, "outputs": [], "source": [ "# so you need to use a timedelta64\n", "dt64_2 = dt64 + np.timedelta64(500, \"ms\")" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.934963Z", "iopub.status.busy": "2025-01-09T18:05:38.934629Z", "iopub.status.idle": "2025-01-09T18:05:38.938909Z", "shell.execute_reply": "2025-01-09T18:05:38.938322Z" } }, "outputs": [ { "data": { "text/plain": [ "0 2012-04-04 14:21:42.800\n", "1 2012-04-04 14:18:37.500\n", "2 2012-04-04 14:08:46.500\n", "Name: time, dtype: datetime64[ns]" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# This is, of course, also the case with datetime64 columns\n", "df[\"time\"] + np.timedelta64(500, \"ms\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ObsPlus provides some common converters between ObsPy and Pandas/Numpy time objects." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Converting between ObsPy and Numpy time datatypes\n", "Conversion between ObsPy `UTCDateTime` and numpy `datetime64` objects when needed." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.940871Z", "iopub.status.busy": "2025-01-09T18:05:38.940438Z", "iopub.status.idle": "2025-01-09T18:05:38.943402Z", "shell.execute_reply": "2025-01-09T18:05:38.942934Z" } }, "outputs": [], "source": [ "# Convert a time column to an array of ObsPy UTCDateTime objects\n", "utc_array = to_utc(df[\"time\"])" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.945068Z", "iopub.status.busy": "2025-01-09T18:05:38.944890Z", "iopub.status.idle": "2025-01-09T18:05:38.947932Z", "shell.execute_reply": "2025-01-09T18:05:38.947450Z" } }, "outputs": [], "source": [ "# Convert back to datetime64\n", "dt64_array = to_datetime64(utc_array)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**Note**: `datetime64` arrays are much more efficient in terms of memory usage and computational efficiency than arrays of `UTCDateTime` *objects*.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you prefer not to manually define `timedelta64` to perform offsets, `to_timedelta64` simply converts a real number to an offset in seconds." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "execution": { "iopub.execute_input": "2025-01-09T18:05:38.949757Z", "iopub.status.busy": "2025-01-09T18:05:38.949580Z", "iopub.status.idle": "2025-01-09T18:05:38.953355Z", "shell.execute_reply": "2025-01-09T18:05:38.952874Z" } }, "outputs": [ { "data": { "text/plain": [ "numpy.timedelta64(3255000000,'ns')" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "to_timedelta64(3.255)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.16" } }, "nbformat": 4, "nbformat_minor": 4 }