Initial push of long term tracking notebook

notebooks/SITCOM-1310_Long_Tracking_Drift_06Apr24.ipynb

@@ -0,0 +1,398 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "d9e95b1e-5d96-4ceb-a634-45d830a2e26e",
+   "metadata": {},
+   "source": [
+    "# TMA Analysis code supporting SITCOM-1310\n",
+    "Craig Lage - 07-Apr-24\n",
+    "\n",
+    "During several nights, the TMA was set to track for long period, form 10 minutes to over an hour.  Then the StarTrackers were used to determine how well the TMA was tracking the sky.  This long term drift could be due to TMA deficiencies, or to deficiencies in the pointing model."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5d4f88a6-0f11-4fee-bd62-e1ab33acaf08",
+   "metadata": {},
+   "source": [
+    "# Prepare the notebook"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ce557368-629b-4b0b-b3a5-e6bba79813c7",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "# Directory to store the data\n",
+    "from pathlib import Path\n",
+    "dataDir = Path(\"./plots\")\n",
+    "dataDir.mkdir(exist_ok=True, parents=True)\n",
+    "\n",
+    "# These are the nights when the tests were run\n",
+    "dayObss = [20240403, 20240404, 20240405]\n",
+    "#dayObss = [20240417,20240418]\n",
+    "\n",
+    "# Some sequences give NaNs for some reason.  Skipping those.\n",
+    "# This is the dayObs, and the first seqNum of the exposure series.\n",
+    "# These seem to be ones at the beginning or end of the night\n",
+    "onesToSkip = [[20240403, 98], [20240404, 15], [20240404, 17], [20240404, 1042]]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "647883db-3024-479d-80d6-3ebd3a8b9d29",
+   "metadata": {},
+   "source": [
+    "# Note that the fast camera code requires using the \n",
+    "# tickets/DM-43774 branch of summit_extras."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1e44aafe-c1bf-4aff-b979-510ab19b71a7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys, time, os, asyncio, glob\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
+    "from matplotlib.backends.backend_pdf import PdfPages\n",
+    "from astropy.time import Time, TimeDelta\n",
+    "from lsst.summit.utils.utils import dayObsIntToString\n",
+    "from lsst_efd_client import EfdClient\n",
+    "from lsst.summit.utils.tmaUtils import TMAEventMaker, getAzimuthElevationDataForEvent, filterBadValues\n",
+    "\n",
+    "from lsst.summit.extras.fastStarTrackerAnalysis import getRegularSequences\n",
+    "from lsst.summit.extras.fastStarTrackerAnalysis import findFastStarTrackerImageSources\n",
+    "\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dbc206d0-33b3-492c-8e3f-2cb21ea86c03",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "client = EfdClient(\"usdf_efd\")\n",
+    "eventMaker = TMAEventMaker()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4b6069aa-8ba1-4145-b32a-f7d4dd170ff9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def findMatchingImages(event, jsonData, camera):\n",
+    "    # This finds StarTracker images that were\n",
+    "    # taken during a long slew by matching RA and Dec\n",
+    "    prePad = 30.0\n",
+    "    postPad = 30.0\n",
+    "    begin = track.begin.value \n",
+    "    end = track.end.value\n",
+    "    midTime = (begin + end) / 2.0\n",
+    "    for i in range(len(jsonData)):\n",
+    "        # Find the exposure closest to the midpoint\n",
+    "        if np.isnan(jsonData.iloc[i][f'MJD {camera}']) or np.isnan(jsonData.iloc[i][f'Ra {camera}']):\n",
+    "            continue\n",
+    "        thisTime = Time(jsonData.iloc[i][f'MJD {camera}'], format='mjd').unix_tai - 37.0\n",
+    "        if abs(thisTime - midTime) < 30.0:\n",
+    "            break\n",
+    "    ra = jsonData.iloc[i][f'Ra {camera}']\n",
+    "    dec = jsonData.iloc[i][f'Dec {camera}']\n",
+    "    seqNums = []\n",
+    "    for i in range(len(jsonData)):\n",
+    "        if np.isnan(jsonData.iloc[i][f'Ra {camera}']):\n",
+    "            continue\n",
+    "        thisTime = Time(jsonData.iloc[i][f'MJD {camera}'], format='mjd').unix_tai - 37.0\n",
+    "        if (thisTime > begin - prePad) and (thisTime < end + postPad):\n",
+    "            dra = abs(ra - jsonData.iloc[i][f'Ra {camera}'])\n",
+    "            ddec = abs(dec - jsonData.iloc[i][f'Dec {camera}'])\n",
+    "            #print(i, dra, ddec)\n",
+    "            if (dra < 0.1) and (ddec < 0.1):\n",
+    "                seqNums.append(jsonData.index[i])\n",
+    "    return seqNums\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d2f7ee8f-f171-49f8-8f10-69f11e5d2685",
+   "metadata": {},
+   "source": [
+    "# This uses the json metadata from RubinTV, which contains the astrometric solutions for the StarTracker exposures.  You need to download the metadata for each dayObs and put it in a location where the program can find it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b53c950c-2a61-4031-8414-1c62e1275a60",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# This disables the annoying logging associated with filterBadValues\n",
+    "import logging.config\n",
+    "logging.config.dictConfig({\n",
+    "    'version': 1,\n",
+    "    'disable_existing_loggers': True,\n",
+    "})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d503c178-f472-4c55-a3e2-3a4b134ea56a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fastPlateScale = 0.67 # arcsec/pixel\n",
+    "pdf = PdfPages(dataDir / \"Long_Term_Tracking_Events_19Apr24.pdf\")\n",
+    "fig = plt.figure(figsize=(15,10))\n",
+    "\n",
+    "# The lists below will be used to plot the \n",
+    "# tracking drift vs AltAz\n",
+    "alts = []\n",
+    "azs = []\n",
+    "raDrifts = []\n",
+    "decDrifts = []\n",
+    "\n",
+    "for dayObs in dayObss:\n",
+    "    dayObsString = dayObsIntToString(dayObs)\n",
+    "    filename = f'../json_metadata/startracker_narrow_{dayObsString}.json'\n",
+    "    jsonData = pd.read_json(filename).T\n",
+    "    jsonData = jsonData.sort_index()\n",
+    "    events = eventMaker.getEvents(dayObs)\n",
+    "    longTracks = [e for e in events if (e.type.name == 'TRACKING' and e.duration > 500)]\n",
+    "    files = getRegularSequences(dayObs) # This is for the fast camera data\n",
+    "\n",
+    "    for n, track in enumerate(longTracks):\n",
+    "        seqNums = []\n",
+    "        seqNums = findMatchingImages(track, jsonData, 'narrow')\n",
+    "        print(f\"For long track beginning at {track.begin.isot},there are {len(seqNums)} narrow images that match\")\n",
+    "        if len(seqNums) == 0:\n",
+    "            continue\n",
+    "        # This is for the encoder data        \n",
+    "        azimuthData, elevationData = getAzimuthElevationDataForEvent(client, track)\n",
+    "        azTimes = azimuthData['timestamp'].values\n",
+    "        azTimes -= azTimes[0]\n",
+    "        azErrors = azimuthData['azError'].values\n",
+    "        filterBadValues(azErrors, maxConsecutiveValues=10)\n",
+    "        elTimes = elevationData['timestamp'].values\n",
+    "        elTimes -= elTimes[0]\n",
+    "        elErrors = elevationData['elError'].values\n",
+    "        filterBadValues(elErrors, maxConsecutiveValues=10)\n",
+    "\n",
+    "        # Now this is for the narrow data\n",
+    "        times = []\n",
+    "        ras = []\n",
+    "        decs = []\n",
+    "        for seqNum in seqNums:\n",
+    "            try:\n",
+    "                time = Time(jsonData[jsonData.index == seqNum]['MJD narrow'].values[0], format='mjd').unix_tai\n",
+    "                ra = jsonData[jsonData.index == seqNum]['Calculated Ra narrow'].values[0]\n",
+    "                dec = jsonData[jsonData.index == seqNum]['Calculated Dec narrow'].values[0]\n",
+    "                ras.append(ra)\n",
+    "                decs.append(dec)\n",
+    "                times.append(time)\n",
+    "            except:\n",
+    "                continue\n",
+    "\n",
+    "        if [dayObs, seqNums[0]] in onesToSkip:\n",
+    "            continue\n",
+    "        azs.append(np.median(azimuthData['actualPosition'].values))\n",
+    "        alts.append(np.median(elevationData['actualPosition'].values))\n",
+    "        \n",
+    "        times = np.array(times)\n",
+    "        times -= times[0]\n",
+    "        ras = np.array(ras)\n",
+    "        ras = (ras - ras[0]) * 3600.0\n",
+    "        decs = np.array(decs)\n",
+    "        decs = (decs - decs[0]) * 3600.0\n",
+    "        raDrift = (ras[-1] - ras[0]) / times[-1] * 30.0\n",
+    "        decDrift = (decs[-1] - decs[0]) / times[-1] * 30.0\n",
+    "        if np.isnan(raDrift) or np.isnan(decDrift):\n",
+    "            continue\n",
+    "        raDrifts.append(raDrift)\n",
+    "        decDrifts.append(decDrift)\n",
+    "        axs = fig.subplots(2,3)\n",
+    "        plt.suptitle(f\"Long term tracking drift {dayObsString}: {seqNums[0]} - {seqNums[-1]}  RA = {ra:.2f}, Dec = {dec:.2f}\", \\\n",
+    "                     fontsize=18)\n",
+    "        plt.subplots_adjust(wspace=0.3)\n",
+    "        axs[0][0].set_title(f\"Narrow RA drift = {raDrift:.2f} arcseconds/30 seconds\")\n",
+    "        axs[0][0].scatter(times, ras)\n",
+    "        axs[0][0].set_ylabel(\"RA Drift (arcseconds)\")\n",
+    "        axs[0][0].set_xlabel(\"Time (seconds)\")\n",
+    "        axs[0][1].set_title(f\"Narrow Dec drift = {decDrift:.2f} arcseconds/30 seconds\")\n",
+    "        axs[0][1].scatter(times, decs)\n",
+    "        axs[0][1].set_ylabel(\"Dec Drift (arcseconds)\")\n",
+    "        axs[0][1].set_xlabel(\"Time (seconds)\")\n",
+    "        axs[0][2].set_title(\"Encoder tracking errors\")\n",
+    "        axs[0][2].plot(azTimes, azErrors, color='green', label='Az')\n",
+    "        axs[0][2].plot(elTimes, elErrors, color='red', label='El')\n",
+    "        axs[0][2].set_ylabel(\"Errors (arcsec)\")\n",
+    "        axs[0][2].set_xlabel(\"Time (sec)\")\n",
+    "        axs[0][2].set_ylim(-1.0,1.0)\n",
+    "        axs[0][2].legend()\n",
+    "\n",
+    "        # Now for the fast data\n",
+    "        # These are needed to make sure you don't grab an extra \n",
+    "        # seqNum before or after the long track.\n",
+    "        seqNums = findMatchingImages(track, jsonData, 'fast')\n",
+    "        print(f\"For long track beginning at {track.begin.isot},there are {len(seqNums)} fast images that match\")\n",
+    "        if len(seqNums) == 0:\n",
+    "            continue\n",
+    "        cenXs = []\n",
+    "        cenYs = []\n",
+    "        fastTimes = []\n",
+    "\n",
+    "        for seqNum in seqNums:\n",
+    "            try:\n",
+    "                filename = [b for (a,b) in files if a == seqNum][0]\n",
+    "                result = findFastStarTrackerImageSources(filename, boxSize=50, streaming=False)\n",
+    "                fastTime = Time(jsonData[jsonData.index == seqNum]['MJD fast'].values[0], format='mjd').unix_tai\n",
+    "                cenXs.append(result[0].centroidX)\n",
+    "                cenYs.append(result[0].centroidY)\n",
+    "                fastTimes.append(fastTime)\n",
+    "            except:\n",
+    "                continue\n",
+    "        \n",
+    "        if len(cenXs) == 0:\n",
+    "            print(\"Fast images failed to solve\")\n",
+    "            axs[1][0].set_title(\"No valid fast data\")\n",
+    "            axs[1][1].set_title(\"No valid fast data\")\n",
+    "            axs[1][2].axis('off')\n",
+    "        else:\n",
+    "            fastTimes = np.array(fastTimes)\n",
+    "            fastTimes -= fastTimes[0]\n",
+    "            cenXs = np.array(cenXs) * fastPlateScale\n",
+    "            cenYs = np.array(cenYs) * fastPlateScale\n",
+    "            xFit = np.polyfit(fastTimes, cenXs, 1)\n",
+    "            yFit = np.polyfit(fastTimes, cenYs, 1)\n",
+    "            xPlot = np.linspace(fastTimes[0], fastTimes[-1], 100)\n",
+    "            yPlotXcen = np.polyval(xFit, xPlot)\n",
+    "            yPlotYcen = np.polyval(yFit, xPlot)\n",
+    "            xDrift = (yPlotXcen[-1] - yPlotXcen[0]) / (xPlot[-1] - xPlot[0]) * 30.0\n",
+    "            yDrift = (yPlotYcen[-1] - yPlotYcen[0]) / (xPlot[-1] - xPlot[0]) * 30.0\n",
+    "            axs[1][0].set_title(f\"Fast X drift = {xDrift:.2f} arcseconds/30 seconds\")\n",
+    "            axs[1][0].scatter(fastTimes, cenXs)\n",
+    "            axs[1][0].plot(xPlot, yPlotXcen, color='red', ls = '--')\n",
+    "            axs[1][0].set_ylabel(\"X Drift (arcseconds)\")\n",
+    "            axs[1][0].set_xlabel(\"Time (seconds)\")\n",
+    "            axs[1][1].set_title(f\"Fast Y drift = {yDrift:.2f} arcseconds/30 seconds\")\n",
+    "            axs[1][1].scatter(fastTimes, cenYs)\n",
+    "            axs[1][1].plot(xPlot, yPlotYcen, color='red', ls = '--')\n",
+    "            axs[1][1].set_ylabel(\"Y Drift (arcseconds)\")\n",
+    "            axs[1][1].set_xlabel(\"Time (seconds)\")\n",
+    "            axs[1][2].axis('off')\n",
+    "        pdf.savefig(fig)  # saves the current figure into a pdf page\n",
+    "        plt.clf()\n",
+    "pdf.close()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cd75064d-23e4-4377-b7e9-6d254dd90979",
+   "metadata": {},
+   "source": [
+    "# This cell plots the drift errors vs AltAz, using the size of the plot marker to indicate the tracking drift.\n",
+    "# This only applies to the narrow camera."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c749da28-fefe-4687-9583-cfb13376b578",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "radAzs = np.radians(azs)\n",
+    "fig=plt.figure(figsize=(16,8))\n",
+    "ax1 = plt.subplot(121, projection='polar')\n",
+    "ax1.invert_yaxis()\n",
+    "ax1.set_title(\"RA drift - Size proportional to errors\\n Blue = pos, Red = Neg\", fontsize=18)\n",
+    "radAzsPos = []\n",
+    "radAzsNeg = []\n",
+    "altsPos = []\n",
+    "altsNeg = []\n",
+    "raDriftsPos = []\n",
+    "raDriftsNeg = []\n",
+    "for n, raDrift in enumerate(raDrifts):\n",
+    "    if raDrift > 0.0:\n",
+    "        radAzsPos.append(radAzs[n])\n",
+    "        altsPos.append(alts[n])\n",
+    "        raDriftsPos.append(raDrift)\n",
+    "    else:\n",
+    "        radAzsNeg.append(radAzs[n])\n",
+    "        altsNeg.append(alts[n])\n",
+    "        raDriftsNeg.append(raDrift)\n",
+    "raSizePos = abs(np.array(raDriftsPos)) * 500.0\n",
+    "raSizeNeg = abs(np.array(raDriftsNeg)) * 500.0\n",
+    "\n",
+    "ax1.scatter(radAzsPos, altsPos, color='blue', marker='o', s=raSizePos)\n",
+    "ax1.scatter(radAzsNeg, altsNeg, color='red', marker='o', s=raSizeNeg)\n",
+    "\n",
+    "ax2 = plt.subplot(122, projection='polar')\n",
+    "ax2.invert_yaxis()\n",
+    "ax2.set_title(\"Dec drift - Size proportional to errors\\n Blue = pos, Red = Neg\", fontsize=18)\n",
+    "radAzsPos = []\n",
+    "radAzsNeg = []\n",
+    "altsPos = []\n",
+    "altsNeg = []\n",
+    "decDriftsPos = []\n",
+    "decDriftsNeg = []\n",
+    "for n, decDrift in enumerate(decDrifts):\n",
+    "    if decDrift > 0.0:\n",
+    "        radAzsPos.append(radAzs[n])\n",
+    "        altsPos.append(alts[n])\n",
+    "        decDriftsPos.append(decDrift)\n",
+    "    else:\n",
+    "        radAzsNeg.append(radAzs[n])\n",
+    "        altsNeg.append(alts[n])\n",
+    "        decDriftsNeg.append(decDrift)\n",
+    "decSizePos = abs(np.array(decDriftsPos)) * 500.0\n",
+    "decSizeNeg = abs(np.array(decDriftsNeg)) * 500.0\n",
+    "ax2.scatter(radAzsPos, altsPos, color='blue', marker='o', s=decSizePos)\n",
+    "ax2.scatter(radAzsNeg, altsNeg, color='red', marker='o', s=decSizeNeg)\n",
+    "plt.savefig(dataDir / \"Long_Term_Tracking_Drifts_AltAz_06Apr24.pdf\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "150f438b-77d4-4db9-b026-e9e6e2f18864",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}