diff --git a/python/lsst/summit/utils/efdUtils.py b/python/lsst/summit/utils/efdUtils.py index d6d3512c..0f46250a 100644 --- a/python/lsst/summit/utils/efdUtils.py +++ b/python/lsst/summit/utils/efdUtils.py @@ -832,6 +832,8 @@ def getCommands( timeKey = time.to_pydatetime() if timeKey in commandTimes: - raise ValueError(f"There is already a command at {timeKey=} - make a better data structure!") + msg = f"There is already a command at {timeKey=} - make a better data structure!" + msg += f"Colliding commands = {commandTimes[timeKey]} and {command}" + raise ValueError(msg) commandTimes[timeKey] = command return commandTimes diff --git a/python/lsst/summit/utils/simonyi/__init__.py b/python/lsst/summit/utils/simonyi/__init__.py new file mode 100644 index 00000000..83c18f36 --- /dev/null +++ b/python/lsst/summit/utils/simonyi/__init__.py @@ -0,0 +1,22 @@ +# This file is part of summit_utils. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +from .mountData import * diff --git a/python/lsst/summit/utils/simonyi/mountAnalysis.py b/python/lsst/summit/utils/simonyi/mountAnalysis.py new file mode 100644 index 00000000..5fad628a --- /dev/null +++ b/python/lsst/summit/utils/simonyi/mountAnalysis.py @@ -0,0 +1,359 @@ +# This file is part of summit_utils. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +from __future__ import annotations + +__all__ = ["calculateMountErrors", "plotMountErrors"] + +import logging +from dataclasses import dataclass +from typing import TYPE_CHECKING +from zoneinfo import ZoneInfo + +import matplotlib.dates as mdates +import matplotlib.pyplot as plt +import numpy as np +from matplotlib.dates import num2date +from matplotlib.ticker import FuncFormatter + +from lsst.summit.utils.tmaUtils import filterBadValues +from lsst.summit.utils.utils import dayObsIntToString + +from .mountData import getAzElRotDataForExposure + +if TYPE_CHECKING: + from lsst_efd_client import EfdClient + from matplotlib.figure import Figure + + from lsst.daf.butler import DimensionRecord + + from .mountData import MountData + + +NON_TRACKING_IMAGE_TYPES = ["BIAS", "FLAT"] + +COMCAM_ANGLE_TO_EDGE_OF_FIELD_ARCSEC = 1800.0 +LSSTCAM_ANGLE_TO_EDGE_OF_FIELD_ARCSEC = 8500.0 + +# These levels determine the colouring of the cells in the RubinTV. +# Yellow for warning level, red for bad level +MOUNT_IMAGE_WARNING_LEVEL = 0.05 +MOUNT_IMAGE_BAD_LEVEL = 0.10 # and red for this + +N_REPLACED_WARNING_LEVEL = 999999 # fill these values in once you've spoken to Craig and Brian +N_REPLACED_BAD_LEVEL = 999999 # fill these values in once you've spoken to Craig and Brian + + +@dataclass +class MountErrors: + azRms: float + elRms: float + rotRms: float + imageAzRms: float + imageElRms: float + imageRotRms: float + imageImpactRms: float + residualFiltering: bool + nReplacedAz: int + nReplacedEl: int + + +def tickFormatter(value: float, tick_number: float) -> str: + # Convert the value to a string without subtracting large numbers + # tick_number is unused. + return f"{value:.2f}" + + +def calculateMountErrors( + expRecord: DimensionRecord, + client: EfdClient, + maxDelta=0.1, + doFilterResiduals=True, +) -> tuple[MountErrors, MountData] | tuple[None, None]: + """Queries EFD for a given exposure and calculates the RMS errors in the + axes during the exposure, optionally plotting and saving the data. + """ + logger = logging.getLogger(__name__) + + imgType = expRecord.observation_type.upper() + if imgType in NON_TRACKING_IMAGE_TYPES: + logger.info(f"Skipping mount torques for non-tracking image type {imgType} for {expRecord.id}") + return None, None + + mountData = getAzElRotDataForExposure(client, expRecord) + + elevation = 90 - expRecord.zenith_angle + + azError = mountData.azimuthData["azError"].to_numpy() + elError = mountData.elevationData["elError"].to_numpy() + rotError = mountData.rotationData["rotError"].to_numpy() + nReplacedAz = 0 + nReplacedEl = 0 + if doFilterResiduals: + # Filtering out bad values + nReplacedAz = filterBadValues(azError, maxDelta) + nReplacedEl = filterBadValues(elError, maxDelta) + mountData.azimuthData["azError"] = azError + mountData.elevationData["elError"] = elError + azRms = np.sqrt(np.mean(azError * azError)) + elRms = np.sqrt(np.mean(elError * elError)) + rotRms = np.sqrt(np.mean(rotError * rotError)) + + # Calculate Image impact RMS + imageAzRms = azRms * np.cos(elevation * np.pi / 180.0) + imageElRms = elRms + imageRotRms = rotRms * COMCAM_ANGLE_TO_EDGE_OF_FIELD_ARCSEC * np.pi / 180.0 / 3600.0 + imageImpactRms = np.sqrt(imageAzRms**2 + imageElRms**2 + imageRotRms**2) + + mountErrors = MountErrors( + azRms=azRms, + elRms=elRms, + rotRms=rotRms, + imageAzRms=imageAzRms, + imageElRms=imageElRms, + imageRotRms=imageRotRms, + imageImpactRms=imageImpactRms, + residualFiltering=doFilterResiduals, + nReplacedAz=nReplacedAz, + nReplacedEl=nReplacedEl, + ) + + return (mountErrors, mountData) + + +def plotMountErrors( + mountData: MountData, + mountErrors: MountErrors, + figure=None, + saveFilename: str = "", +) -> Figure: + imageImpactRms = mountErrors.imageImpactRms + expRecord = mountData.expRecord + if expRecord is not None: + dayObsString = dayObsIntToString(expRecord.day_obs) + dataIdString = f"{expRecord.instrument} {dayObsString} - seqNum {expRecord.seq_num}" + title = f"{dataIdString} - Exposure time = {expRecord.exposure_time:.1f}s" + + if figure is None: + figure = plt.figure(figsize=(12, 8)) + else: + figure.clear() + ax = figure.gca() + ax.clear() + + utc = ZoneInfo("UTC") + chile_tz = ZoneInfo("America/Santiago") + + # Function to convert UTC to Chilean time + def offset_time_aware(utc_time): + # Ensure the time is timezone-aware in UTC + if utc_time.tzinfo is None: + utc_time = utc.localize(utc_time) + return utc_time.astimezone(chile_tz) + + [[ax1, ax4], [ax2, ax5], [ax3, ax6]] = figure.subplots( + 3, + 2, + sharex="col", + sharey=False, + gridspec_kw={"wspace": 0.25, "hspace": 0, "height_ratios": [2.5, 1, 1], "width_ratios": [1.5, 1]}, + ) + # [ax1, ax4] = [azimuth, rotator] + # [ax2, ax5] = [azError, rotError] + # [ax3, ax6] = [azTorque, rotTorque] + + # Use the native color cycle for the lines. Because they're on + # different axes they don't cycle by themselves + axs = [ax1, ax2, ax3, ax4, ax5, ax6] + lineColors = [p["color"] for p in plt.rcParams["axes.prop_cycle"]] + nColors = len(lineColors) + colorCounter = 0 + + ax1.plot( + mountData.azimuthData["actualPosition"], + label="Azimuth position", + c=lineColors[colorCounter % nColors], + ) + colorCounter += 1 + ax1.yaxis.set_major_formatter(FuncFormatter(tickFormatter)) + ax1.set_ylabel("Azimuth (degrees)") + + ax1_twin = ax1.twinx() + ax1_twin.plot( + mountData.elevationData["actualPosition"], + label="Elevation position", + c=lineColors[colorCounter % nColors], + ) + colorCounter += 1 + + ax2.plot( + mountData.azimuthData["azError"], + label="Azimuth tracking error", + c=lineColors[colorCounter % nColors], + ) + colorCounter += 1 + ax2.plot( + mountData.elevationData["elError"], + label="Elevation tracking error", + c=lineColors[colorCounter % nColors], + ) + colorCounter += 1 + ax2.axhline(0.01, ls="-.", color="black") + ax2.axhline(-0.01, ls="-.", color="black") + ax2.yaxis.set_major_formatter(FuncFormatter(tickFormatter)) + ax2.set_ylabel("Tracking error (arcsec)") + ax2.set_xticks([]) # remove x tick labels on the hidden upper x-axis + ax2.set_ylim(-0.05, 0.05) + ax2.set_yticks([-0.04, -0.02, 0.0, 0.02, 0.04]) + ax2.legend() + ax2.text(0.1, 0.9, f"Image impact RMS = {imageImpactRms:.3f} arcsec (with rot).", transform=ax2.transAxes) + if mountErrors.residualFiltering: + ax2.text( + 0.1, + 0.8, + ( + f"{mountErrors.nReplacedAz} bad az values and " + f"{mountErrors.nReplacedEl} bad el values were replaced" + ), + transform=ax2.transAxes, + ) + ax3_twin = ax3.twinx() + ax3.plot( + mountData.azimuthData["actualTorque"], + label="Azimuth torque", + c=lineColors[colorCounter % nColors], + ) + colorCounter += 1 + ax3_twin.plot( + mountData.elevationData["actualTorque"], + label="Elevation torque", + c=lineColors[colorCounter % nColors], + ) + colorCounter += 1 + ax3.set_ylabel("Azimuth torque (Nm)") + ax3_twin.set_ylabel("Elevation torque (Nm)") + ax3.set_xlabel("Time (UTC)") # yes, it really is UTC, matplotlib converts this automatically! + + # put the ticks at an angle, and right align with the tick marks + ax3.set_xticks(ax3.get_xticks()) # needed to supress a user warning + xlabels = ax3.get_xticks() + ax3.set_xticklabels(xlabels) + ax3.tick_params(axis="x", rotation=45) + ax3.xaxis.set_major_locator(mdates.AutoDateLocator()) + ax3.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S")) + + ax4.plot( + mountData.rotationData["actualPosition"], + label="Rotator position", + c=lineColors[colorCounter % nColors], + ) + colorCounter += 1 + ax4.yaxis.set_major_formatter(FuncFormatter(tickFormatter)) + ax4.yaxis.tick_right() + ax4.set_ylabel("Rotator angle (degrees)") + ax4.yaxis.set_label_position("right") + ax5.plot( + mountData.rotationData["rotError"], + c=lineColors[colorCounter % nColors], + ) + + colorCounter += 1 + ax5.axhline(0.1, ls="-.", color="black") + ax5.axhline(-0.1, ls="-.", color="black") + ax5.yaxis.set_major_formatter(FuncFormatter(tickFormatter)) + ax5.set_ylabel("Tracking error (arcsec)") + ax5.tick_params(labelbottom=False) # Hide x-axis tick labels without removing ticks + ax5.set_ylim(-0.5, 0.5) + ax5.set_yticks([-0.4, -0.2, 0.0, 0.2, 0.4]) + ax5.yaxis.tick_right() + ax5.yaxis.set_label_position("right") + + ax6.plot(mountData.rotationTorques["torque0"], label="Torque0", c=lineColors[colorCounter % nColors]) + colorCounter += 1 + ax6.plot(mountData.rotationTorques["torque1"], label="Torque1", c=lineColors[colorCounter % nColors]) + ax6.set_ylabel("Rotator torque (Nm)") + ax6.set_xlabel("Time (UTC)") # yes, it really is UTC, matplotlib converts this automatically! + # put the ticks at an angle, and right align with the tick marks + ax6.set_xticks(ax6.get_xticks()) # needed to supress a user warning + xlabels = ax6.get_xticks() + ax6.set_xticklabels(xlabels) + ax6.tick_params(axis="x", rotation=45) + ax6.xaxis.set_major_locator(mdates.AutoDateLocator()) + ax6.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S")) + ax6.yaxis.tick_right() + ax6.yaxis.set_label_position("right") + ax6.legend() + + ax1_twin.yaxis.set_major_formatter(FuncFormatter(tickFormatter)) + ax1_twin.set_ylabel("Elevation (degrees)") + ax1.tick_params(labelbottom=False) # Hide x-axis tick labels without removing ticks + # combine the legends and put inside the plot + handles1a, labels1a = ax1.get_legend_handles_labels() + handles1b, labels1b = ax1_twin.get_legend_handles_labels() + handles2a, labels2a = ax3.get_legend_handles_labels() + handles2b, labels2b = ax3_twin.get_legend_handles_labels() + handles = handles1a + handles1b + handles2a + handles2b + labels = labels1a + labels1b + labels2a + labels2b + # ax2 is "in front" of ax1 because it has the vlines plotted on it, and + # vlines are on ax2 so that they appear at the bottom of the legend, so + # make sure to plot the legend on ax2, otherwise the vlines will go on + # top of the otherwise-opaque legend. + ax1_twin.legend(handles, labels, facecolor="white", framealpha=1) + + ax1.set_title("Azimuth and Elevation") + ax4.set_title("Rotator") + figure.subplots_adjust(top=0.85) # Adjust the top margin to make room for the suptitle + figure.suptitle(title, fontsize=14, y=1.04) # Adjust y to move the title up + + # Create the upper axis for Chilean time + ax1_twiny = ax1.twiny() + ax1_twiny.set_xlim(ax1.get_xlim()) # Set the limits of the upper axis to match the lower axis + utcTicks = ax1.get_xticks() # Use the same ticks as the lower UTC axis + utcTickLabels = [num2date(tick, tz=utc) for tick in utcTicks] + chileTickLabels = [offset_time_aware(label) for label in utcTickLabels] + # Set the same tick positions but with Chilean time labels + ax1_twiny.set_xticks(utcTicks) + ax1_twiny.set_xticklabels([tick.strftime("%H:%M:%S") for tick in chileTickLabels]) + ax1_twiny.tick_params(axis="x", rotation=45) + ax1_twiny.set_xlabel("Time (Chilean Time)") + + ax4_twiny = ax4.twiny() + ax4_twiny.set_xlim(ax4.get_xlim()) # Set the limits of the upper axis to match the lower axis + utcTicks = ax4.get_xticks() # Use the same ticks as the lower UTC axis + utcTickLabels = [num2date(tick, tz=utc) for tick in utcTicks] + chileTickLabels = [offset_time_aware(label) for label in utcTickLabels] + # Set the same tick positions but with Chilean time labels + ax4_twiny.set_xticks(utcTicks) + ax4_twiny.set_xticklabels([tick.strftime("%H:%M:%S") for tick in chileTickLabels]) + ax4_twiny.tick_params(axis="x", rotation=45) + ax4_twiny.set_xlabel("Time (Chilean Time)") + + # Add exposure start and end: + for ax in axs: + if expRecord is not None: + # assert expRecord is not None, "expRecord is None" + ax.axvline(expRecord.timespan.begin.utc.datetime, ls="--", color="green") + ax.axvline(expRecord.timespan.end.utc.datetime, ls="--", color="red") + + if saveFilename: + figure.savefig(saveFilename, bbox_inches="tight") + + return figure diff --git a/python/lsst/summit/utils/simonyi/mountData.py b/python/lsst/summit/utils/simonyi/mountData.py new file mode 100644 index 00000000..df75f4d3 --- /dev/null +++ b/python/lsst/summit/utils/simonyi/mountData.py @@ -0,0 +1,115 @@ +# This file is part of summit_utils. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING + +from ..efdUtils import getEfdData + +if TYPE_CHECKING: + from astropy.time import Time + from lsst_efd_client import EfdClient + from pandas import DataFrame + + from lsst.daf.butler import DimensionRecord + + +@dataclass +class MountData: + begin: Time + end: Time + azimuthData: DataFrame + elevationData: DataFrame + rotationData: DataFrame + rotationTorques: DataFrame + includedPrePadding: float + includedPostPadding: float + expRecord: DimensionRecord | None + + +def getAzElRotDataForPeriod( + client: EfdClient, begin: Time, end: Time, prePadding: float = 0, postPadding: float = 0 +) -> MountData: + azimuthData = getEfdData( + client, + "lsst.sal.MTMount.azimuth", + begin=begin, + end=end, + prePadding=prePadding, + postPadding=postPadding, + ) + elevationData = getEfdData( + client, + "lsst.sal.MTMount.elevation", + begin=begin, + end=end, + prePadding=prePadding, + postPadding=postPadding, + ) + rotationData = getEfdData( + client, + "lsst.sal.MTRotator.rotation", + begin=begin, + end=end, + prePadding=prePadding, + postPadding=postPadding, + ) + rotationTorques = getEfdData( + client, + "lsst.sal.MTRotator.motors", + begin=begin, + end=end, + prePadding=prePadding, + postPadding=postPadding, + ) + + azValues = azimuthData["actualPosition"].to_numpy() + elValues = elevationData["actualPosition"].to_numpy() + rotValues = rotationData["actualPosition"].to_numpy() + azDemand = azimuthData["demandPosition"].to_numpy() + elDemand = elevationData["demandPosition"].to_numpy() + rotDemand = rotationData["demandPosition"].to_numpy() + + azError = (azValues - azDemand) * 3600 + elError = (elValues - elDemand) * 3600 + rotError = (rotValues - rotDemand) * 3600 + + azimuthData["azError"] = azError + elevationData["elError"] = elError + rotationData["rotError"] = rotError + + mountData = MountData( + begin, end, azimuthData, elevationData, rotationData, rotationTorques, prePadding, postPadding, None + ) + return mountData + + +def getAzElRotDataForExposure( + client: EfdClient, expRecord: DimensionRecord, prePadding: float = 0, postPadding: float = 0 +) -> MountData: + + begin = expRecord.timespan.begin + end = expRecord.timespan.end + mountData = getAzElRotDataForPeriod(client, begin, end, prePadding, postPadding) + mountData.expRecord = expRecord + return mountData