diff --git a/notebooks/tel_and_site/subsys_req_ver/m1m3/SITCOM-1799_wo_vms_Adding_Earthquake_events_to_the_TN81_on_vibration_on_M1M3.ipynb b/notebooks/tel_and_site/subsys_req_ver/m1m3/SITCOM-1799_wo_vms_Adding_Earthquake_events_to_the_TN81_on_vibration_on_M1M3.ipynb new file mode 100644 index 00000000..dba62cdf --- /dev/null +++ b/notebooks/tel_and_site/subsys_req_ver/m1m3/SITCOM-1799_wo_vms_Adding_Earthquake_events_to_the_TN81_on_vibration_on_M1M3.ipynb @@ -0,0 +1,607 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "72cef637-a558-4af8-99ad-fd87380468dd", + "metadata": {}, + "source": [ + "# SITCOM-1399: Adding Earthquake events without needs of VMS data, to the TN81, on vibration on M1M3\n", + "\n", + "2025-01-08 HyeYun Park" + ] + }, + { + "cell_type": "markdown", + "id": "cccba396-fc91-4462-8622-0a8665e43ef1", + "metadata": {}, + "source": [ + "**Descripción**\n", + "\n", + "Sub-ticket related to the [SITCOM-918: Write technote on hardpoint oscillations during tma slews](https://rubinobs.atlassian.net/browse/SITCOM-918), to study earthquake events seen by M1M3. \n", + "We can read the TN81 [here](https://sitcomtn-081.lsst.io/v/SITCOM-918/index.html).\n", + "\n", + "This notebook is to update SITCOM-1399 to avoid downloading big data of VMS, and to update new event below which was big and near from the Rubin observatory which faulted the M1M3 during the event.\n", + "- 2024-11-26 UTC 12:33:09" + ] + }, + { + "cell_type": "markdown", + "id": "a9ca43c3-0c39-42fa-95b4-82100be4937a", + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-12T15:47:36.513702Z", + "iopub.status.busy": "2024-06-12T15:47:36.513386Z", + "iopub.status.idle": "2024-06-12T15:47:36.515840Z", + "shell.execute_reply": "2024-06-12T15:47:36.515467Z", + "shell.execute_reply.started": "2024-06-12T15:47:36.513688Z" + } + }, + "source": [ + "## Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb9bbe64-0eac-4125-9469-2d5c6be88201", + "metadata": {}, + "outputs": [], + "source": [ + "import sys, time, os, asyncio\n", + "import scipy.stats as stats\n", + "from scipy.signal import find_peaks\n", + "from scipy import signal\n", + "import h5py\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "from astropy.time import Time\n", + "from lsst.summit.utils.tmaUtils import TMAEventMaker, TMAState\n", + "from lsst.summit.utils.efdUtils import getEfdData, makeEfdClient, clipDataToEvent, calcNextDay\n", + "import matplotlib.dates as mdates\n", + "%matplotlib inline\n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7bf0c802-0745-4a76-9250-9740c3118c82", + "metadata": {}, + "outputs": [], + "source": [ + "# Functions to get data\n", + "# Modified functions of ticket SITCOM-761\n", + "\n", + "key_m1m3_dict={'1 X': 'm1m3_x_1', \n", + " '1 Y': 'm1m3_y_1', \n", + " '1 Z': 'm1m3_z_1', \n", + " '2 X': 'm1m3_x_2', \n", + " '2 Y': 'm1m3_z_2', # note these two have been \n", + " '2 Z': 'm1m3_y_2', # switched pending SUMMIT-7911\n", + " '3 X': 'm1m3_x_3', \n", + " '3 Y': 'm1m3_y_3', \n", + " '3 Z': 'm1m3_z_3'\n", + " }\n", + "\n", + "def get_efd_data(begin, end, client):\n", + "\n", + " \"\"\"Extract all the MTMount data from the EFD and add to dict.\n", + "\n", + " Args:\n", + " begin (str): The start time of the query.\n", + " end (str): The end time of the query.\n", + " client (object): influx client\n", + "\n", + " Returns:\n", + " dict: A dictionary containing the MTMount data.\n", + " \"\"\"\n", + "\n", + " query_dict = {}\n", + "\n", + " query_dict[\"el\"] = getEfdData(\n", + " client,\n", + " \"lsst.sal.MTMount.elevation\",\n", + " columns=[\"private_sndStamp\", \"private_efdStamp\", \"actualPosition\", \"actualVelocity\", \"actualTorque\"],\n", + " begin=begin,\n", + " end=end,\n", + " prePadding=0,\n", + " postPadding=0,\n", + " warn=False,\n", + " )\n", + " query_dict[\"az\"] = getEfdData(\n", + " client,\n", + " \"lsst.sal.MTMount.azimuth\",\n", + " columns=[\"private_sndStamp\", \"private_efdStamp\", \"actualPosition\", \"actualVelocity\", \"actualTorque\"],\n", + " begin=begin,\n", + " end=end,\n", + " prePadding=0,\n", + " postPadding=0,\n", + " warn=False,\n", + " )\n", + " return query_dict\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c5f9c19e-ec3e-4314-8c47-865bc24acb58", + "metadata": {}, + "outputs": [], + "source": [ + "def get_freq_psd(vals, timestep):\n", + "\n", + " \"\"\"\n", + " Calculates the frequency power spectrum of a signal.\n", + "\n", + " Args:\n", + " vals (np.array): The signal values.\n", + " timestep (float): The time step between samples.\n", + "\n", + " Returns:\n", + " tuple: The frequencies and power spectral density.\n", + " \"\"\"\n", + "\n", + " # Remove the mean from the signal.\n", + " meanval = np.mean(vals)\n", + " signal = vals - meanval\n", + "\n", + " # Calculate the length of the signal.\n", + " N = len(signal)\n", + "\n", + " # Calculate the power spectral density.\n", + " psd = np.abs(np.fft.rfft(np.array(signal) * 1)) ** 2\n", + "\n", + " # Calculate the frequencies.\n", + " frequencies = np.fft.rfftfreq(N, timestep)\n", + "\n", + " return (frequencies, psd)\n", + "\n", + "def get_peak_points(freq, psd, height=0.01):\n", + " \"\"\"\n", + " Get the peak points of the power spectral density (PSD).\n", + "\n", + " Args:\n", + " freq (numpy.ndarray): The frequency vector.\n", + " psd (numpy.ndarray): The power spectral density.\n", + " height (float): The minimum peak height.\n", + "\n", + " Returns:\n", + " numpy.ndarray: The peak points.\n", + " \"\"\"\n", + "\n", + " # Find the peak indices and heights.\n", + " peak_ind, peak_dict = find_peaks(psd, height=height)\n", + " peaks = freq[peak_ind]\n", + "\n", + " # If there are no peaks, return None.\n", + " if len(peaks) < 1:\n", + " return None\n", + "\n", + " # Find the sub-peaks within each group of peaks that are close in frequency.\n", + " points = []\n", + " for i, peak in enumerate(peaks):\n", + " sel = (abs(peaks - peak) < 1)\n", + " sub_peaks = peaks[sel]\n", + " sub_heights = peak_dict['peak_heights'][sel]\n", + " points.append(sub_peaks[np.argmax(sub_heights)])\n", + "\n", + " # Return the unique peak points.\n", + " return np.unique(np.array(points))" + ] + }, + { + "cell_type": "markdown", + "id": "5fce44ff-7543-4217-be8e-3d9124807629", + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T07:40:08.325036Z", + "iopub.status.busy": "2024-06-18T07:40:08.324549Z", + "iopub.status.idle": "2024-06-18T07:40:08.326931Z", + "shell.execute_reply": "2024-06-18T07:40:08.326607Z", + "shell.execute_reply.started": "2024-06-18T07:40:08.325021Z" + } + }, + "source": [ + "## Analysis for the evening's event: 2024-11-26 UTC 12:33:09" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f809ec10-b1b6-4b92-86a0-79b8371038c9", + "metadata": {}, + "outputs": [], + "source": [ + "client=makeEfdClient()\n", + "begin_time=Time('2024-11-26 12:30:00', format=\"iso\", scale=\"utc\")\n", + "end_time=Time('2024-11-26 12:40:00', format=\"iso\", scale=\"utc\")\n", + "efd_dict=get_efd_data(begin_time, end_time, client)" + ] + }, + { + "cell_type": "markdown", + "id": "9b94d778-3b5f-43f8-8142-1c4ed0ece84c", + "metadata": {}, + "source": [ + "### HP Forces" + ] + }, + { + "cell_type": "markdown", + "id": "7dd81833-2f57-4a75-beb3-c4a3316375d0", + "metadata": {}, + "source": [ + "The limit required to keep the mirror safe during a mag 6 earthquake must be less than 3000 N. Let's study the and HP forces along the same axes and compare them with the accelerations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea011ad6-5a29-4c6e-add3-76e4ab0dad6b", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Select data from a given date\n", + "dayObs = 20241126\n", + "eventMaker = TMAEventMaker()\n", + "event = eventMaker.getEvents(dayObs)" + ] + }, + { + "cell_type": "markdown", + "id": "1ddd10cf-2685-4a7e-88a6-9fd3c6c817e2", + "metadata": { + "execution": { + "iopub.execute_input": "2024-08-05T07:08:38.001985Z", + "iopub.status.busy": "2024-08-05T07:08:38.001698Z", + "iopub.status.idle": "2024-08-05T07:08:38.048269Z", + "shell.execute_reply": "2024-08-05T07:08:38.047860Z", + "shell.execute_reply.started": "2024-08-05T07:08:38.001971Z" + } + }, + "source": [ + "#### Earthquake events using the EFD client" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f5a4064-167e-4668-8792-997062ee9c3c", + "metadata": {}, + "outputs": [], + "source": [ + "from lsst_efd_client import EfdClient\n", + "client = EfdClient('usdf_efd', db_name=\"lsst.backpack\")\n", + "\n", + "topics = await client.get_topics()\n", + "topics" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4bd701c7-f7fb-4319-830a-c8b80359fa93", + "metadata": {}, + "outputs": [], + "source": [ + "query = f'''SELECT * FROM \"lsst.backpack.usgs_earthquake_data\"'''\n", + "await client.influx_client.query(query)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ee720434-41c5-4057-ab8f-ef016e0ecfd5", + "metadata": {}, + "outputs": [], + "source": [ + "data = getEfdData(client, \"lsst.backpack.usgs_earthquake_data\", dayObs=20241126)\n", + "data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "62ff90a4-e840-485c-8e83-240db265e49b", + "metadata": {}, + "outputs": [], + "source": [ + "# Pick which earthquake to have a look. This case first ([0]) one on the list is picked.\n", + "data.index[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2dd62c55-434b-4f84-8448-c88042ff61d8", + "metadata": {}, + "outputs": [], + "source": [ + "client = EfdClient(\"usdf_efd\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0619035a-cba8-4730-9805-b3a42dd10c92", + "metadata": {}, + "outputs": [], + "source": [ + "df_hp = getEfdData(\n", + " client, \"lsst.sal.MTM1M3.hardpointActuatorData\", begin=begin_time, end=end_time\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fbaf4239-899b-4876-a6bb-252db87eb035", + "metadata": {}, + "outputs": [], + "source": [ + "def compare_mount_hardpoints(\n", + " df_mtmount_ele,\n", + " df_mtmount_azi,\n", + " df_hp,\n", + " begin,\n", + " end,\n", + "):\n", + " fig, axs = plt.subplots(1, 1, dpi=125, figsize=(15, 8))\n", + " ax = axs # [0]\n", + " df_plot = df_hp[\"measuredForce0\"][begin:end]\n", + " ax.plot(df_plot, color=\"red\", lw=\"0.5\", label=\"HP 0\")\n", + " df_plot = df_hp[\"measuredForce1\"][begin:end]\n", + " ax.plot(df_plot, color=\"blue\", lw=\"0.5\", label=\"HP 1\")\n", + " df_plot = df_hp[\"measuredForce2\"][begin:end]\n", + " ax.plot(df_plot, color=\"black\", lw=\"0.5\", label=\"HP 2\")\n", + " df_plot = df_hp[\"measuredForce3\"][begin:end]\n", + " ax.plot(df_plot, color=\"green\", lw=\"0.5\", label=\"HP 3\")\n", + " df_plot = df_hp[\"measuredForce4\"][begin:end]\n", + " ax.plot(df_plot, color=\"orange\", lw=\"0.5\", label=\"HP 4\")\n", + " df_plot = df_hp[\"measuredForce5\"][begin:end]\n", + " ax.plot(df_plot, color=\"yellow\", lw=\"0.5\", label=\"HP 5\")\n", + " ax.set_ylabel(\"HP Force \\n[N]\")\n", + " # ax.legend()\n", + " ax2 = ax.twinx()\n", + " # ax = axs[1]\n", + " df_plot = df_mtmount_ele[\"actualPosition\"][begin:end]\n", + " ax2.plot(df_plot, color=\"green\", lw=\"0.5\", label=\"el\")\n", + " # ax.axvline(begin, lw=\"0.5\", c=\"k\", label=\"Slew start\")\n", + " # ax.axvline(end, lw=\"0.5\", c=\"b\", label=\"Slew stop\")\n", + "\n", + " # ax = axs[2]\n", + " df_plot = df_mtmount_azi[\"actualPosition\"][begin:end]\n", + " ax2.plot(df_plot, color=\"red\", lw=\"0.5\", label=\"az\")\n", + " ax2.set_ylabel(\"TMA El & Azimuth \\nPosition\\n[deg]\")\n", + "\n", + " ax2.set_xlabel(\"UTC\")\n", + " ax2.legend()\n", + " fig.autofmt_xdate()\n", + " fig.subplots_adjust(hspace=1)\n", + " fig.suptitle(t0)\n", + " fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4e64e96d-bc23-457e-82c6-6a845bd6283a", + "metadata": {}, + "outputs": [], + "source": [ + "df_hp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ae422435-89a2-45da-b999-15ae781d02c7", + "metadata": {}, + "outputs": [], + "source": [ + "df_hp.index" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffebe162-eef2-4160-a43d-5ed4f1fe4b30", + "metadata": {}, + "outputs": [], + "source": [ + "def compare_mount_hardpoints(\n", + " df_mtmount_ele,\n", + " df_mtmount_azi,\n", + " df_hp,\n", + " begin,\n", + " end\n", + "):\n", + " fig, axs = plt.subplots(1, 1, dpi=125, figsize=(15, 8))\n", + " ax = axs # [0]\n", + " df_plot = df_hp[\"measuredForce0\"][begin:end]\n", + " ax.plot(df_plot, color=\"red\", lw=\"0.5\", label=\"HP 0\")\n", + " df_plot = df_hp[\"measuredForce1\"][begin:end]\n", + " ax.plot(df_plot, color=\"blue\", lw=\"0.5\", label=\"HP 1\")\n", + " df_plot = df_hp[\"measuredForce2\"][begin:end]\n", + " ax.plot(df_plot, color=\"black\", lw=\"0.5\", label=\"HP 2\")\n", + " df_plot = df_hp[\"measuredForce3\"][begin:end]\n", + " ax.plot(df_plot, color=\"green\", lw=\"0.5\", label=\"HP 3\")\n", + " df_plot = df_hp[\"measuredForce4\"][begin:end]\n", + " ax.plot(df_plot, color=\"orange\", lw=\"0.5\", label=\"HP 4\")\n", + " df_plot = df_hp[\"measuredForce5\"][begin:end]\n", + " ax.plot(df_plot, color=\"yellow\", lw=\"0.5\", label=\"HP 5\")\n", + " ax.set_ylabel(\"HP Force \\n[N]\")\n", + " # ax.legend()\n", + " ax2 = ax.twinx()\n", + " # ax = axs[1]\n", + " df_plot = df_mtmount_ele[\"actualPosition\"][begin:end]\n", + " ax2.plot(df_plot, color=\"green\", lw=\"0.5\", label=\"el\")\n", + " # ax.axvline(begin, lw=\"0.5\", c=\"k\", label=\"Slew start\")\n", + " # ax.axvline(end, lw=\"0.5\", c=\"b\", label=\"Slew stop\")\n", + "\n", + " # ax = axs[2]\n", + " df_plot = df_mtmount_azi[\"actualPosition\"][begin:end]\n", + " ax2.plot(df_plot, color=\"red\", lw=\"0.5\", label=\"az\")\n", + " ax2.set_ylabel(\"TMA El & Azimuth \\nPosition\\n[deg]\")\n", + "\n", + " ax2.set_xlabel(\"UTC\")\n", + " ax2.legend()\n", + " fig.autofmt_xdate()\n", + " fig.subplots_adjust(hspace=1)\n", + " fig.suptitle(f\"Earthquake detected at {data.index[0]}\")\n", + " fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "060d5cf1-fed9-408c-a053-48b8d36f42c3", + "metadata": {}, + "outputs": [], + "source": [ + "df_mtmount_ele = getEfdData(\n", + " client,\n", + " \"lsst.sal.MTMount.elevation\", begin=begin_time, end=end_time\n", + ")\n", + "df_mtmount_azi = getEfdData(\n", + " client,\n", + " \"lsst.sal.MTMount.azimuth\", begin=begin_time, end=end_time\n", + ")\n", + "df_hp = getEfdData(\n", + " client, \"lsst.sal.MTM1M3.hardpointActuatorData\", begin=begin_time, end=end_time)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8959321e-a0b4-490c-bf89-b4f14fbf827b", + "metadata": {}, + "outputs": [], + "source": [ + "start_slew=data.index[0]\n", + "t0 = pd.to_datetime(start_slew.value, utc=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0060d06f-1614-4f77-8561-c021417c2d78", + "metadata": {}, + "outputs": [], + "source": [ + "import datetime\n", + "t1=(t0-datetime.timedelta(seconds=100)).strftime('%Y-%m-%dT%H:%M:%SZ')\n", + "t2=(t0+datetime.timedelta(seconds=260)).strftime('%Y-%m-%dT%H:%M:%SZ')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "648a4f85-5b97-4003-ba38-ea1ac28cd6c2", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "compare_mount_hardpoints(\n", + " df_mtmount_ele,\n", + " df_mtmount_azi,\n", + " df_hp,\n", + " t1,\n", + " t2\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf18328e-14c8-4bb0-a8df-f47e8f565be1", + "metadata": {}, + "outputs": [], + "source": [ + "def compare_mount_hardpoints_separate_plots(\n", + " df_mtmount_ele,\n", + " df_mtmount_azi,\n", + " df_hp,\n", + " begin,\n", + " end,\n", + "):\n", + " fig, axs = plt.subplots(3, 1, dpi=125, figsize=(15, 8))\n", + " ax = axs[0]\n", + " df_plot = df_hp[\"measuredForce0\"][begin:end]\n", + " ax.plot(df_plot, color=\"red\", lw=\"0.5\", label=\"HP 0\")\n", + " df_plot = df_hp[\"measuredForce1\"][begin:end]\n", + " ax.plot(df_plot, color=\"blue\", lw=\"0.5\", label=\"HP 1\")\n", + " df_plot = df_hp[\"measuredForce2\"][begin:end]\n", + " ax.plot(df_plot, color=\"black\", lw=\"0.5\", label=\"HP 2\")\n", + " df_plot = df_hp[\"measuredForce3\"][begin:end]\n", + " ax.plot(df_plot, color=\"green\", lw=\"0.5\", label=\"HP 3\")\n", + " df_plot = df_hp[\"measuredForce4\"][begin:end]\n", + " ax.plot(df_plot, color=\"orange\", lw=\"0.5\", label=\"HP 4\")\n", + " df_plot = df_hp[\"measuredForce5\"][begin:end]\n", + " ax.plot(df_plot, color=\"yellow\", lw=\"0.5\", label=\"HP 5\")\n", + " ax.set_ylabel(\"HP Force \\n[N]\")\n", + "\n", + " ax = axs[1]\n", + " df_plot = df_mtmount_ele[\"actualPosition\"][begin:end]\n", + " ax.plot(df_plot, color=\"green\", lw=\"0.5\")\n", + " ax.set_ylabel(\"TMA Elevation \\nPosition\\n[deg]\")\n", + " # ax.axvline(begin, lw=\"0.5\", c=\"k\", label=\"Slew start\")\n", + " # ax.axvline(end, lw=\"0.5\", c=\"b\", label=\"Slew stop\")\n", + "\n", + " ax = axs[2]\n", + " df_plot = df_mtmount_azi[\"actualPosition\"][begin:end]\n", + " ax.plot(df_plot, color=\"red\", lw=\"0.5\")\n", + " ax.set_ylabel(\"TMA Azimuth \\nPosition\\n[deg]\")\n", + "\n", + " ax.set_xlabel(\"UTC\")\n", + " fig.autofmt_xdate()\n", + " fig.subplots_adjust(hspace=1)\n", + " fig.suptitle(t0)\n", + " fig.legend()\n", + " fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0208ff6-9ba6-415b-9677-af20c08dbb26", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "compare_mount_hardpoints_separate_plots(\n", + " df_mtmount_ele,\n", + " df_mtmount_azi,\n", + " df_hp,\n", + " t1,\n", + " t2,\n", + ")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "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.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}