Coverage for python / lsst / summit / extras / ringssSeeing.py: 0%
211 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 09:28 +0000
1# This file is part of summit_extras.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22from __future__ import annotations
24import logging
25from dataclasses import dataclass
26from typing import TYPE_CHECKING
27from zoneinfo import ZoneInfo
29import pandas as pd
30from astropy.time import Time, TimeDelta
31from matplotlib.dates import DateFormatter, num2date
33from lsst.summit.utils.dateTime import getDayObsEndTime, getDayObsStartTime
34from lsst.summit.utils.efdUtils import getEfdData, getMostRecentRowWithDataBefore
35from lsst.utils.plotting.figures import make_figure
37HAS_EFD_CLIENT = True
38try:
39 from lsst_efd_client import EfdClient
40except ImportError:
41 HAS_EFD_CLIENT = False
43if TYPE_CHECKING:
44 from matplotlib.figure import Figure
45 from pandas import DataFrame, Series
47 from lsst.daf.butler import Butler, DataCoordinate, DimensionRecord
49RINGSS_TOPIC = "lsst.sal.ESS.logevent_ringssMeasurement"
52@dataclass
53class SeeingConditions:
54 """Class to hold the seeing conditions from the RINGSS instrument.
56 Attributes
57 ----------
58 timestamp : `Time`
59 The time of the seeing conditions.
60 XXX : `timeDelta`
61 Distance to the nearest actual reading, given to give a guide for
62 how accurate the data is likely to be.
63 seeing : `float`
64 XXX fill this in
65 freeAtmSeeing : `float`
66 XXX fill this in
67 groundLayer : `float`
68 XXX fill this in
69 eRMS : `float`
70 XXX fill this in
71 flux : `float`
72 XXX fill this in
73 fwhmFree : `float`
74 XXX fill this in
75 fwhmScintillation : `float`
76 XXX fill this in
77 fwhmSector : `float`
78 XXX fill this in
79 hrNum : `float`
80 XXX fill this in
81 tau0 : `float`
82 XXX fill this in
83 theta0 : `float`
84 XXX fill this in
85 totalVariance : `float`
86 XXX fill this in
87 turbulenceProfiles0 : `float`
88 XXX fill this in
89 turbulenceProfiles1 : `float`
90 XXX fill this in
91 turbulenceProfiles2 : `float`
92 XXX fill this in
93 turbulenceProfiles3 : `float`
94 XXX fill this in
95 turbulenceProfiles4 : `float`
96 XXX fill this in
97 turbulenceProfiles5 : `float`
98 XXX fill this in
99 turbulenceProfiles6 : `float`
100 XXX fill this in
101 turbulenceProfiles7 : `float`
102 XXX fill this in
103 wind : `float`
104 XXX fill this in
105 zenithDistance : `float`
106 XXX fill this in
107 """
109 timestamp: Time
111 eRMS: float = float("nan")
112 flux: float = float("nan")
113 fwhmFree: float = float("nan")
114 fwhmScintillation: float = float("nan")
115 fwhmSector: float = float("nan")
116 hrNum: float = float("nan")
117 tau0: float = float("nan")
118 theta0: float = float("nan")
119 totalVariance: float = float("nan")
120 profile0m: float = float("nan")
121 profile250m: float = float("nan")
122 profile500m: float = float("nan")
123 profile1000m: float = float("nan")
124 profile2000m: float = float("nan")
125 profile4000m: float = float("nan")
126 profile8000m: float = float("nan")
127 profile16000m: float = float("nan")
128 windSpeed: float = float("nan")
129 zenithDistance: float = float("nan")
131 @property
132 def isoparalacticAngle(self) -> float:
133 """Alias for isoparalacticAngle."""
134 return self.isoparalacticAngle
136 @property
137 def starName(self) -> str:
138 """Alias for starName including the HD part."""
139 return f"HD{self.starName}"
141 @property
142 def seeing(self) -> float:
143 """Alias for fwhmScintillation - the seeing in arcsec."""
144 return self.fwhmScintillation
146 @property
147 def seeing2(self) -> float:
148 """The seeing profile adjusted weight in arcsec."""
149 return self.fwhmSector
151 @property
152 def freeAtmosphericSeeing(self) -> float:
153 """Alias for fwhmFree - the free atmospheric seeing in arcsec."""
154 return self.fwhmFree
156 @property
157 def groundLayerSeeing(self) -> float:
158 """The transformation of turbulenceProfiles0 to the ground layer seeing
159 in arcsec.
160 """
161 raise NotImplementedError("I need the transormation to get this value in arcsec")
163 def __init__(self, rows: list[Series]) -> None:
164 # Ensure we have valid data
165 assert len(rows) >= 1, "Must provide some data!"
166 assert len(rows) <= 2, "Provided data must be either or two rows."
168 # Extract timestamp from index
169 timestamps = [row.name for row in rows]
170 self.timestamp = Time(timestamps[0]) # Use the first timestamp as default
172 if len(rows) == 1:
173 # If only one row, use its values directly
174 row = rows[0]
176 self.eRMS = row.get("eRMS", float("nan"))
177 self.flux = row.get("flux", float("nan"))
178 self.fwhmFree = row.get("fwhmFree", float("nan"))
179 self.fwhmScintillation = row.get("fwhmScintillation", float("nan"))
180 self.fwhmSector = row.get("fwhmSector", float("nan"))
181 self.hrNum = row.get("hrNum", float("nan"))
182 self.tau0 = row.get("tau0", float("nan"))
183 self.theta0 = row.get("theta0", float("nan"))
184 self.totalVariance = row.get("totalVariance", float("nan"))
185 self.profile0m = row.get("turbulenceProfiles0", float("nan"))
186 self.profile250m = row.get("turbulenceProfiles1", float("nan"))
187 self.profile500m = row.get("turbulenceProfiles2", float("nan"))
188 self.profile1000m = row.get("turbulenceProfiles3", float("nan"))
189 self.profile2000m = row.get("turbulenceProfiles4", float("nan"))
190 self.profile4000m = row.get("turbulenceProfiles5", float("nan"))
191 self.profile8000m = row.get("turbulenceProfiles6", float("nan"))
192 self.profile16000m = row.get("turbulenceProfiles7", float("nan"))
193 self.windSpeed = row.get("wind", float("nan"))
194 self.zenithDistance = row.get("zenithDistance", float("nan"))
195 else:
196 # Interpolate between two rows
197 t1, t2 = timestamps
198 # Calculate the midpoint timestamp for interpolation
199 t = t1 + (t2 - t1) / 2
200 self.timestamp = Time(t)
202 # Weight for interpolation (0 to 1)
203 w = (t - t1) / (t2 - t1)
205 row1, row2 = rows[0], rows[1]
207 # Helper function for interpolation
208 def interpolate(col):
209 if col in row1 and col in row2:
210 v1, v2 = row1.get(col, float("nan")), row2.get(col, float("nan"))
211 return v1 + (v2 - v1) * w
212 return float("nan")
214 self.eRMS = interpolate("eRMS")
215 self.flux = interpolate("flux")
216 self.fwhmFree = interpolate("fwhmFree")
217 self.fwhmScintillation = interpolate("fwhmScintillation")
218 self.fwhmSector = interpolate("fwhmSector")
219 self.hrNum = interpolate("hrNum")
220 self.tau0 = interpolate("tau0")
221 self.theta0 = interpolate("theta0")
222 self.totalVariance = interpolate("totalVariance")
223 self.profile0m = interpolate("turbulenceProfiles0")
224 self.profile250m = interpolate("turbulenceProfiles1")
225 self.profile500m = interpolate("turbulenceProfiles2")
226 self.profile1000m = interpolate("turbulenceProfiles3")
227 self.profile2000m = interpolate("turbulenceProfiles4")
228 self.profile4000m = interpolate("turbulenceProfiles5")
229 self.profile8000m = interpolate("turbulenceProfiles6")
230 self.profile16000m = interpolate("turbulenceProfiles7")
231 self.windSpeed = interpolate("wind")
232 self.zenithDistance = interpolate("zenithDistance")
234 def __repr__(self):
235 return (
236 f"SeeingConditions @ {self.timestamp.isot}\n"
237 f" Seeing = XXX define me!\n"
238 f' Free Atm Seeing = {self.fwhmFree:.2f}"\n'
239 # XXX replace this with self.groundLayerSeeing
240 f" Ground layer = {self.profile0m}\n"
241 f" Wind speed = {self.windSpeed:.2f}m/s\n"
242 )
245class RingssSeeingMonitor:
246 def __init__(
247 self, efdClient: EfdClient, warningThreshold: float = 300, errorThreshold: float = 600
248 ) -> None:
249 self.warningThreshold = warningThreshold
250 self.errorThreshold = errorThreshold
251 self.log = logging.getLogger(__name__)
252 self.client = efdClient
254 def getSeeingAtTime(self, time: Time) -> SeeingConditions:
255 """Get the seeing conditions at a specific time.
257 Parameters
258 ----------
259 time : `Time`
260 The time at which to get the seeing conditions.
262 Returns
263 -------
264 seeing : `SeeingConditions`
265 The seeing conditions at the specified time.
266 """
267 begin = time - TimeDelta(self.errorThreshold, format="sec")
268 end = time + TimeDelta(self.errorThreshold, format="sec")
269 data = getEfdData(self.client, RINGSS_TOPIC, begin=begin, end=end, warn=False)
270 if data.empty:
271 raise ValueError(
272 f"Failed to find a RINGSS seeing measurement within"
273 f" {self.errorThreshold / 60:.1f} minutes of {time.isot}"
274 )
276 # Ensure the timestamp is sorted
277 data.sort_index(inplace=True)
279 # Check if the *exact* time exists - seems unlikely, but need to check
280 if time in data.index:
281 row = data.loc[time]
282 return SeeingConditions(
283 rows=[row],
284 )
286 # Convert the astropy Time object to a timezone-aware pandas Timestamp
287 time_datetime = pd.Timestamp(time.datetime).tz_localize("UTC")
289 # Use the timezone-aware timestamp for comparison
290 earlier = (
291 data[data.index < time_datetime].iloc[-1] if not data[data.index < time_datetime].empty else None
292 )
293 later = (
294 data[data.index > time_datetime].iloc[0] if not data[data.index > time_datetime].empty else None
295 )
297 if later is None and earlier is not None and (time - earlier.name).sec < self.errorThreshold:
298 self.log.info("Returning the last available value.")
299 return SeeingConditions(
300 rows=[earlier],
301 )
303 if earlier is None or later is None:
304 raise ValueError("Cannot interpolate: insufficient data before or after the requested time.")
306 # Check time difference: to log warnings/raise as necessary
307 earlierTime = earlier.name
308 laterTime = later.name
309 interval = (laterTime - earlierTime).seconds
311 if interval > self.errorThreshold:
312 raise ValueError(
313 f"Requested time {time.isot} would require interpolating between values more "
314 f"than {self.errorThreshold} apart: {interval:.2f} seconds."
315 )
316 if interval > self.warningThreshold:
317 self.log.warning(
318 f"Interpolating between values more than {self.warningThreshold / 60:.1f} mins apart."
319 )
321 return SeeingConditions(
322 rows=[earlier, later],
323 )
325 def getMostRecentTimestamp(self) -> Time:
326 """Get the most recent timestamp for which seeing data is available.
328 Returns
329 -------
330 time : `Time`
331 The most recent timestamp with seeing data.
332 """
333 now = Time.now() + TimeDelta(10, format="sec")
334 row = getMostRecentRowWithDataBefore(
335 self.client,
336 RINGSS_TOPIC,
337 now,
338 maxSearchNMinutes=self.errorThreshold / 60,
339 )
340 return Time(row.name)
342 def getMostRecentSeeing(self) -> SeeingConditions:
343 """Get the most recent seeing conditions.
345 Returns
346 -------
347 seeing : `SeeingConditions`
348 The most recent seeing conditions.
350 Raises
351 ------
352 ValueError
353 If no data is available in the EFD within the error threshold.
354 """
355 now = Time.now() + TimeDelta(10, format="sec")
356 try:
357 row = getMostRecentRowWithDataBefore(
358 self.client,
359 RINGSS_TOPIC,
360 now,
361 maxSearchNMinutes=self.errorThreshold / 60,
362 )
363 except ValueError as e:
364 raise ValueError("Could not get SeeingConditions - no data available in the EFD.") from e
366 # XXX add validity checks here
368 return SeeingConditions([row])
370 def getSeeingForDataId(self, butler: Butler, dataId: DataCoordinate) -> SeeingConditions:
371 """Get the seeing conditions for a specific data ID.
373 Parameters
374 ----------
375 butler : `Butler`
376 The Butler instance to query.
377 dataId : `DataCoordinate`
378 The data ID for which to get the seeing conditions.
380 Returns
381 -------
382 seeing : `SeeingConditions`
383 The seeing conditions for the specified data ID.
384 """
385 (expRecord,) = butler.registry.queryDimensionRecords("exposure", dataId=dataId)
386 return self.getSeeingForExpRecord(expRecord)
388 def getSeeingForExpRecord(self, expRecord: DimensionRecord) -> SeeingConditions:
389 """Get the seeing conditions for a specific exposure record.
391 Parameters
392 ----------
393 expRecord : `DimensionRecord`
394 The exposure record for which to get the seeing conditions.
396 Returns
397 -------
398 seeing : `SeeingConditions`
399 The seeing conditions for the specified exposure record.
400 """
401 midPoint = expRecord.timespan.begin + TimeDelta(expRecord.exposure_time / 2, format="sec")
402 return self.getSeeingAtTime(midPoint)
404 def plotSeeingForDayObs(
405 self, dayObs: int, addMostRecentBox: bool = True, fig: Figure | None = None
406 ) -> Figure:
407 """Plot the seeing conditions for a specific day observation.
409 Parameters
410 ----------
411 dayObs : `int`
412 The dayObs to plot the seeing for, in YYYYMMDD format.
413 addMostRecentBox : `bool`, optional
414 Whether to add a box with the most recent seeing conditions.
415 fig : `Figure`, optional
416 The figure to plot on. If None, a new figure will be created.
418 Returns
419 -------
420 fig : `Figure`
421 The figure with the plotted seeing conditions.
422 """
423 startTime = getDayObsStartTime(dayObs)
424 endTime = getDayObsEndTime(dayObs)
425 data = getEfdData(self.client, RINGSS_TOPIC, begin=startTime, end=endTime, warn=False)
426 fig = self.plotSeeing(data, addMostRecentBox=addMostRecentBox, fig=fig)
427 return fig
429 def plotSeeing(
430 self, dataframe: DataFrame, addMostRecentBox: bool = True, fig: Figure | None = None
431 ) -> Figure:
432 """Plot the seeing conditions from a DataFrame.
434 Parameters
435 ----------
436 dataframe : `DataFrame`
437 The DataFrame containing the seeing conditions data.
438 addMostRecentBox : `bool`, optional
439 Whether to add a box with the most recent seeing conditions.
440 fig : `Figure`, optional
441 The figure to plot on. If None, a new figure will be created.
443 Returns
444 -------
445 fig : `Figure`
446 The figure with the plotted seeing conditions.
447 """
448 ls = "-"
449 ms = "o"
450 df = dataframe
452 if df.empty:
453 raise ValueError("No data to plot for the given time range.")
455 seeings = [SeeingConditions([row]) for _, row in df.iterrows()]
457 if fig is None:
458 fig = make_figure(figsize=(18, 10))
459 ax1 = fig.add_subplot(111)
460 else:
461 fig.clear()
462 ax1 = fig.add_subplot(111)
464 utc = ZoneInfo("UTC")
465 chile_tz = ZoneInfo("America/Santiago")
467 # Function to convert UTC to Chilean time
468 def offset_time_aware(utc_time):
469 # Ensure the time is timezone-aware in UTC
470 if utc_time.tzinfo is None:
471 utc_time = utc.localize(utc_time)
472 return utc_time.astimezone(chile_tz)
474 df.index = pd.DatetimeIndex([t for t in df.index])
476 ax1.plot([seeing.fwhmFree for seeing in seeings], "b", label='Free atmos. seeing"', ls=ls, marker=ms)
477 ax1.plot([seeing.seeing for seeing in seeings], "r", label='Seeing "', ls=ls, marker=ms)
478 ax1.plot(
479 [seeing.seeing2 for seeing in seeings], "g", label='Profile adjusted seeing "', ls=ls, marker=ms
480 )
482 ax2 = ax1.twiny()
483 ax2.set_xlim(ax1.get_xlim())
485 # Format both axes to show only time
486 ax1.xaxis.set_major_formatter(DateFormatter("%H:%M:%S"))
487 ax2.xaxis.set_major_formatter(DateFormatter("%H:%M:%S"))
489 # Apply the timezone-aware offset to the top axis ticks
490 ax2.set_xticks(ax1.get_xticks())
491 offset_ticks = [offset_time_aware(num2date(tick)) for tick in ax1.get_xticks()]
492 ax2.set_xticklabels([tick.strftime("%H:%M:%S") for tick in offset_ticks])
494 ax1.set_ylim(0, 1.1 * max([s.seeing2 for s in seeings]))
495 ax1.set_xlabel("Time (UTC)")
496 ax2.set_xlabel("Time (Chilean Time)")
497 ax1.set_ylabel("Seeing (arcsec)")
498 ax1.yaxis.grid(True, which="major", linestyle="--", linewidth=0.5, alpha=0.7)
499 ax1.xaxis.grid(True, which="major", linestyle="--", linewidth=0.5, alpha=0.7)
501 # Update legend with larger font size
502 ax1.legend(loc="lower left", fontsize=14)
504 # Calculate current seeing and age of data
505 if addMostRecentBox:
506 currentSeeing = seeings[-1]
507 justTime = currentSeeing.timestamp.isot.split("T")[1].split(".")[0]
509 text = f'Current Seeing: {currentSeeing.seeing:.2f}"\n' f"Last updated @ {justTime} UTC"
510 ax1.text(
511 0.05,
512 0.95,
513 text,
514 transform=ax1.transAxes,
515 fontsize=14,
516 verticalalignment="top",
517 horizontalalignment="left",
518 bbox=dict(boxstyle="round", facecolor="white", alpha=0.7),
519 )
521 return fig