Coverage for python/lsst/summit/extras/ringssSeeing.py: 26%

1# This file is part of summit_extras.

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.

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 numpy as np

30import pandas as pd

31from astropy.time import Time, TimeDelta

32from matplotlib.dates import DateFormatter, num2date

34from lsst.summit.utils.dateTime import getDayObsEndTime, getDayObsStartTime

35from lsst.summit.utils.efdUtils import getEfdData, getMostRecentRowWithDataBefore

36from lsst.utils.plotting.figures import make_figure

38HAS_EFD_CLIENT = True

39try:

40 from lsst_efd_client import EfdClient

41except ImportError:

42 HAS_EFD_CLIENT = False

44if TYPE_CHECKING:

45 from datetime import datetime

47 from matplotlib.figure import Figure

48 from pandas import DataFrame, Series

50 from lsst.daf.butler import Butler, DataCoordinate, DimensionRecord

52RINGSS_TOPIC = "lsst.sal.ESS.logevent_ringssMeasurement"

55@dataclass

56class SeeingConditions:

57 """A single snapshot of RINGSS seeing conditions.

59 Values are either taken directly from a RINGSS EFD row or linearly

60 interpolated between two rows bracketing a requested time. All

61 numeric attributes default to NaN when the underlying column was

62 missing.

64 Attributes

65 ----------

66 timestamp : `astropy.time.Time`

67 The time to which these conditions apply. For an interpolated

68 instance this is the midpoint between the two source rows.

69 eRMS : `float`

70 Residual defocus signature RMS reported by RINGSS.

71 flux : `float`

72 Measured flux from the RINGSS target star, in instrumental

73 units.

74 fwhmFree : `float`

75 Free-atmosphere seeing FWHM, in arcsec.

76 fwhmScintillation : `float`

77 Scintillation-derived total seeing FWHM, in arcsec.

78 fwhmSector : `float`

79 Profile-weighted sector seeing FWHM, in arcsec.

80 hrNum : `float`

81 RINGSS star HR number.

82 tau0 : `float`

83 Atmospheric coherence time, in ms.

84 theta0 : `float`

85 Isoplanatic angle, in arcsec.

86 totalVariance : `float`

87 Total variance of the wavefront signal.

88 profile0m : `float`

89 Turbulence integral at the ground layer (0 m).

90 profile250m : `float`

91 Turbulence integral at 250 m.

92 profile500m : `float`

93 Turbulence integral at 500 m.

94 profile1000m : `float`

95 Turbulence integral at 1 km.

96 profile2000m : `float`

97 Turbulence integral at 2 km.

98 profile4000m : `float`

99 Turbulence integral at 4 km.

100 profile8000m : `float`

101 Turbulence integral at 8 km.

102 profile16000m : `float`

103 Turbulence integral at 16 km.

104 windSpeed : `float`

105 Wind speed at the turbulence layer, in m/s.

106 zenithDistance : `float`

107 Zenith distance of the RINGSS target at the measurement time,

108 in degrees.

109 """

110

111 timestamp: Time

112

113 eRMS: float = float("nan")

114 flux: float = float("nan")

115 fwhmFree: float = float("nan")

116 fwhmScintillation: float = float("nan")

117 fwhmSector: float = float("nan")

118 hrNum: float = float("nan")

119 tau0: float = float("nan")

120 theta0: float = float("nan")

121 totalVariance: float = float("nan")

122 profile0m: float = float("nan")

123 profile250m: float = float("nan")

124 profile500m: float = float("nan")

125 profile1000m: float = float("nan")

126 profile2000m: float = float("nan")

127 profile4000m: float = float("nan")

128 profile8000m: float = float("nan")

129 profile16000m: float = float("nan")

130 windSpeed: float = float("nan")

131 zenithDistance: float = float("nan")

132

133 @property

134 def isoparalacticAngle(self) -> float:

135 """Isoparalactic angle, in arcsec.

136

137 Returns

138 -------

139 angle : `float`

140 The isoparalactic angle.

141

142 Raises

143 ------

144 NotImplementedError

145 Always raised; there is no corresponding RINGSS field

146 available in the EFD topic yet.

147 """

148 raise NotImplementedError("isoparalacticAngle is not available from the RINGSS EFD topic")

149

150 @property

151 def starName(self) -> str:

152 """HR-prefixed name of the RINGSS target star.

153

154 Returns

155 -------

156 name : `str`

157 The star name, e.g. ``"HR1234"``.

158 """

159 return f"HR{int(self.hrNum)}"

160

161 @property

162 def seeing(self) -> float:

163 """Scintillation-derived total seeing FWHM, in arcsec.

164

165 Returns

166 -------

167 seeing : `float`

168 Alias for `fwhmScintillation`.

169 """

170 return self.fwhmScintillation

171

172 @property

173 def seeing2(self) -> float:

174 """Profile-weighted seeing FWHM, in arcsec.

175

176 Returns

177 -------

178 seeing : `float`

179 Alias for `fwhmSector`.

180 """

181 return self.fwhmSector

182

183 @property

184 def freeAtmosphericSeeing(self) -> float:

185 """Free-atmosphere seeing FWHM, in arcsec.

186

187 Returns

188 -------

189 seeing : `float`

190 Alias for `fwhmFree`.

191 """

192 return self.fwhmFree

193

194 @property

195 def groundLayerSeeing(self) -> float:

196 """Ground-layer seeing FWHM in arcsec.

197

198 Returns

199 -------

200 seeing : `float`

201 The seeing contribution from the ground layer.

202

203 Raises

204 ------

205 NotImplementedError

206 Always raised; the conversion from the turbulence

207 integral at 0 m to arcsec is not yet implemented.

208 """

209 raise NotImplementedError("I need the transormation to get this value in arcsec")

210

211 def __init__(self, rows: list[Series], targetTime: Time | None = None) -> None:

212 # Ensure we have valid data

213 assert len(rows) >= 1, "Must provide some data!"

214 assert len(rows) <= 2, "Provided data must be either or two rows."

215

216 # Extract timestamp from index

217 timestamps = [row.name for row in rows]

218 self.timestamp = Time(timestamps[0]) # Use the first timestamp as default

219

220 if len(rows) == 1:

221 # If only one row, use its values directly

222 row = rows[0]

223

224 self.eRMS = row.get("eRMS", float("nan"))

225 self.flux = row.get("flux", float("nan"))

226 self.fwhmFree = row.get("fwhmFree", float("nan"))

227 self.fwhmScintillation = row.get("fwhmScintillation", float("nan"))

228 self.fwhmSector = row.get("fwhmSector", float("nan"))

229 self.hrNum = row.get("hrNum", float("nan"))

230 self.tau0 = row.get("tau0", float("nan"))

231 self.theta0 = row.get("theta0", float("nan"))

232 self.totalVariance = row.get("totalVariance", float("nan"))

233 self.profile0m = row.get("turbulenceProfiles0", float("nan"))

234 self.profile250m = row.get("turbulenceProfiles1", float("nan"))

235 self.profile500m = row.get("turbulenceProfiles2", float("nan"))

236 self.profile1000m = row.get("turbulenceProfiles3", float("nan"))

237 self.profile2000m = row.get("turbulenceProfiles4", float("nan"))

238 self.profile4000m = row.get("turbulenceProfiles5", float("nan"))

239 self.profile8000m = row.get("turbulenceProfiles6", float("nan"))

240 self.profile16000m = row.get("turbulenceProfiles7", float("nan"))

241 self.windSpeed = row.get("wind", float("nan"))

242 self.zenithDistance = row.get("zenithDistance", float("nan"))

243 else:

244 # Interpolate between two rows at the requested target time.

245 # If no target is given, fall back to the midpoint.

246 t1, t2 = timestamps

247 if targetTime is None:

248 t = t1 + (t2 - t1) / 2

249 else:

250 t = pd.Timestamp(targetTime.datetime).tz_localize("UTC")

251 self.timestamp = Time(t)

252

253 # Weight for linear interpolation (0 to 1), clamped to the

254 # bracket so extrapolation beyond the two rows can't occur.

255 w = (t - t1) / (t2 - t1)

256 w = max(0.0, min(1.0, float(w)))

257

258 row1, row2 = rows[0], rows[1]

259

260 # Helper function for interpolation

261 def interpolate(col: str) -> float:

262 if col in row1 and col in row2:

263 v1, v2 = row1.get(col, float("nan")), row2.get(col, float("nan"))

264 return v1 + (v2 - v1) * w

265 return float("nan")

266

267 self.eRMS = interpolate("eRMS")

268 self.flux = interpolate("flux")

269 self.fwhmFree = interpolate("fwhmFree")

270 self.fwhmScintillation = interpolate("fwhmScintillation")

271 self.fwhmSector = interpolate("fwhmSector")

272 self.hrNum = interpolate("hrNum")

273 self.tau0 = interpolate("tau0")

274 self.theta0 = interpolate("theta0")

275 self.totalVariance = interpolate("totalVariance")

276 self.profile0m = interpolate("turbulenceProfiles0")

277 self.profile250m = interpolate("turbulenceProfiles1")

278 self.profile500m = interpolate("turbulenceProfiles2")

279 self.profile1000m = interpolate("turbulenceProfiles3")

280 self.profile2000m = interpolate("turbulenceProfiles4")

281 self.profile4000m = interpolate("turbulenceProfiles5")

282 self.profile8000m = interpolate("turbulenceProfiles6")

283 self.profile16000m = interpolate("turbulenceProfiles7")

284 self.windSpeed = interpolate("wind")

285 self.zenithDistance = interpolate("zenithDistance")

286

287 def __repr__(self) -> str:

288 return (

289 f"SeeingConditions @ {self.timestamp.isot}\n"

290 f" Seeing = XXX define me!\n"

291 f' Free Atm Seeing = {self.fwhmFree:.2f}"\n'

292 # XXX replace this with self.groundLayerSeeing

293 f" Ground layer = {self.profile0m}\n"

294 f" Wind speed = {self.windSpeed:.2f}m/s\n"

295 )

296

297

298class RingssSeeingMonitor:

299 """Query and plot RINGSS seeing data from the EFD.

300

301 Parameters

302 ----------

303 efdClient : `lsst_efd_client.EfdClient`

304 Client used to query the EFD for RINGSS data.

305 warningThreshold : `float`, optional

306 Interpolation interval, in seconds, above which a warning is

307 logged when returning interpolated seeing conditions.

308 errorThreshold : `float`, optional

309 Interpolation interval, in seconds, above which an interpolated

310 value is refused outright. Also used as the search half-window

311 when looking for the nearest RINGSS measurement.

312 """

313

314 def __init__(

315 self, efdClient: EfdClient, warningThreshold: float = 300, errorThreshold: float = 600

316 ) -> None:

317 self.warningThreshold = warningThreshold

318 self.errorThreshold = errorThreshold

319 self.log = logging.getLogger(__name__)

320 self.client = efdClient

321

322 def getSeeingAtTime(self, time: Time) -> SeeingConditions:

323 """Get the seeing conditions at a specific time.

324

325 Parameters

326 ----------

327 time : `Time`

328 The time at which to get the seeing conditions.

329

330 Returns

331 -------

332 seeing : `SeeingConditions`

333 The seeing conditions at the specified time.

334 """

335 begin = time - TimeDelta(self.errorThreshold, format="sec")

336 end = time + TimeDelta(self.errorThreshold, format="sec")

337 data = getEfdData(self.client, RINGSS_TOPIC, begin=begin, end=end, warn=False)

338 if data.empty:

339 raise ValueError(

340 f"Failed to find a RINGSS seeing measurement within"

341 f" {self.errorThreshold / 60:.1f} minutes of {time.isot}"

342 )

343

344 # Ensure the timestamp is sorted

345 data.sort_index(inplace=True)

346

347 # Convert the astropy Time object to a tz-aware pandas Timestamp

348 # so comparisons against ``data.index`` work.

349 time_datetime = pd.Timestamp(time.datetime).tz_localize("UTC")

350

351 # Check if the *exact* time exists - seems unlikely, but need to check

352 if time_datetime in data.index:

353 row = data.loc[time_datetime]

354 return SeeingConditions(

355 rows=[row],

356 )

357

358 # Use the timezone-aware timestamp for comparison

359 earlier = (

360 data[data.index < time_datetime].iloc[-1] if not data[data.index < time_datetime].empty else None

361 )

362 later = (

363 data[data.index > time_datetime].iloc[0] if not data[data.index > time_datetime].empty else None

364 )

365

366 if (

367 later is None

368 and earlier is not None

369 and (time_datetime - earlier.name).total_seconds() < self.errorThreshold

370 ):

371 self.log.info("Returning the last available value.")

372 return SeeingConditions(

373 rows=[earlier],

374 )

375

376 if earlier is None or later is None:

377 raise ValueError("Cannot interpolate: insufficient data before or after the requested time.")

378

379 # Check time difference: to log warnings/raise as necessary

380 earlierTime = earlier.name

381 laterTime = later.name

382 interval = (laterTime - earlierTime).total_seconds()

383

384 if interval > self.errorThreshold:

385 raise ValueError(

386 f"Requested time {time.isot} would require interpolating between values more "

387 f"than {self.errorThreshold} apart: {interval:.2f} seconds."

388 )

389 if interval > self.warningThreshold:

390 self.log.warning(

391 f"Interpolating between values more than {self.warningThreshold / 60:.1f} mins apart."

392 )

393

394 return SeeingConditions(

395 rows=[earlier, later],

396 targetTime=time,

397 )

398

399 def getMostRecentTimestamp(self) -> Time:

400 """Get the most recent timestamp for which seeing data is available.

401

402 Returns

403 -------

404 time : `Time`

405 The most recent timestamp with seeing data.

406 """

407 now = Time.now() + TimeDelta(10, format="sec")

408 row = getMostRecentRowWithDataBefore(

409 self.client,

410 RINGSS_TOPIC,

411 now,

412 maxSearchNMinutes=self.errorThreshold / 60,

413 )

414 return Time(row.name)

415

416 def getMostRecentSeeing(self) -> SeeingConditions:

417 """Get the most recent seeing conditions.

418

419 Returns

420 -------

421 seeing : `SeeingConditions`

422 The most recent seeing conditions.

423

424 Raises

425 ------

426 ValueError

427 If no data is available in the EFD within the error threshold.

428 """

429 now = Time.now() + TimeDelta(10, format="sec")

430 try:

431 row = getMostRecentRowWithDataBefore(

432 self.client,

433 RINGSS_TOPIC,

434 now,

435 maxSearchNMinutes=self.errorThreshold / 60,

436 )

437 except ValueError as e:

438 raise ValueError("Could not get SeeingConditions - no data available in the EFD.") from e

439

440 # XXX add validity checks here

441

442 return SeeingConditions([row])

443

444 def getSeeingForDataId(self, butler: Butler, dataId: DataCoordinate) -> SeeingConditions:

445 """Get the seeing conditions for a specific data ID.

446

447 Parameters

448 ----------

449 butler : `Butler`

450 The Butler instance to query.

451 dataId : `DataCoordinate`

452 The data ID for which to get the seeing conditions.

453

454 Returns

455 -------

456 seeing : `SeeingConditions`

457 The seeing conditions for the specified data ID.

458 """

459 (expRecord,) = butler.registry.queryDimensionRecords("exposure", dataId=dataId)

460 return self.getSeeingForExpRecord(expRecord)

461

462 def getSeeingForExpRecord(self, expRecord: DimensionRecord) -> SeeingConditions:

463 """Get the seeing conditions for a specific exposure record.

464

465 Parameters

466 ----------

467 expRecord : `DimensionRecord`

468 The exposure record for which to get the seeing conditions.

469

470 Returns

471 -------

472 seeing : `SeeingConditions`

473 The seeing conditions for the specified exposure record.

474 """

475 midPoint = expRecord.timespan.begin + TimeDelta(expRecord.exposure_time / 2, format="sec")

476 return self.getSeeingAtTime(midPoint)

477

478 def plotSeeingForDayObs(

479 self, dayObs: int, addMostRecentBox: bool = True, fig: Figure | None = None

480 ) -> Figure:

481 """Plot the seeing conditions for a specific day observation.

482

483 Parameters

484 ----------

485 dayObs : `int`

486 The dayObs to plot the seeing for, in YYYYMMDD format.

487 addMostRecentBox : `bool`, optional

488 Whether to add a box with the most recent seeing conditions.

489 fig : `Figure`, optional

490 The figure to plot on. If None, a new figure will be created.

491

492 Returns

493 -------

494 fig : `Figure`

495 The figure with the plotted seeing conditions.

496 """

497 startTime = getDayObsStartTime(dayObs)

498 endTime = getDayObsEndTime(dayObs)

499 data = getEfdData(self.client, RINGSS_TOPIC, begin=startTime, end=endTime, warn=False)

500 fig = self.plotSeeing(data, addMostRecentBox=addMostRecentBox, fig=fig)

501 return fig

502

503 def plotSeeing(

504 self, dataframe: DataFrame, addMostRecentBox: bool = True, fig: Figure | None = None

505 ) -> Figure:

506 """Plot the seeing conditions from a DataFrame.

507

508 Parameters

509 ----------

510 dataframe : `DataFrame`

511 The DataFrame containing the seeing conditions data.

512 addMostRecentBox : `bool`, optional

513 Whether to add a box with the most recent seeing conditions.

514 fig : `Figure`, optional

515 The figure to plot on. If None, a new figure will be created.

516

517 Returns

518 -------

519 fig : `Figure`

520 The figure with the plotted seeing conditions.

521 """

522 ls = "-"

523 ms = "o"

524 df = dataframe

525

526 if df.empty:

527 raise ValueError("No data to plot for the given time range.")

528

529 seeings = [SeeingConditions([row]) for _, row in df.iterrows()]

530

531 if fig is None:

532 fig = make_figure(figsize=(18, 10))

533 ax1 = fig.add_subplot(111)

534 else:

535 fig.clear()

536 ax1 = fig.add_subplot(111)

537

538 utc = ZoneInfo("UTC")

539 chile_tz = ZoneInfo("America/Santiago")

540

541 # Function to convert UTC to Chilean time

542 def offset_time_aware(utc_time: datetime) -> datetime:

543 # Ensure the time is timezone-aware in UTC

544 if utc_time.tzinfo is None:

545 utc_time = utc_time.replace(tzinfo=utc)

546 return utc_time.astimezone(chile_tz)

547

548 times = pd.DatetimeIndex([t for t in df.index])

549

550 ax1.plot(

551 times,

552 [seeing.fwhmFree for seeing in seeings],

553 "b",

554 label='Free atmos. seeing"',

555 ls=ls,

556 marker=ms,

557 )

558 ax1.plot(times, [seeing.seeing for seeing in seeings], "r", label='Seeing "', ls=ls, marker=ms)

559 ax1.plot(

560 times,

561 [seeing.seeing2 for seeing in seeings],

562 "g",

563 label='Profile adjusted seeing "',

564 ls=ls,

565 marker=ms,

566 )

567

568 ax2 = ax1.twiny()

569 ax2.set_xlim(ax1.get_xlim())

570

571 # Format both axes to show only time

572 ax1.xaxis.set_major_formatter(DateFormatter("%H:%M:%S"))

573 ax2.xaxis.set_major_formatter(DateFormatter("%H:%M:%S"))

574

575 # Apply the timezone-aware offset to the top axis ticks

576 ax2.set_xticks(ax1.get_xticks())

577 offset_ticks = [offset_time_aware(num2date(tick)) for tick in ax1.get_xticks()]

578 ax2.set_xticklabels([tick.strftime("%H:%M:%S") for tick in offset_ticks])

579

580 ax1.set_ylim(0, 1.1 * np.nanmax([s.seeing2 for s in seeings]))

581 ax1.set_xlabel("Time (UTC)")

582 ax2.set_xlabel("Time (Chilean Time)")

583 ax1.set_ylabel("Seeing (arcsec)")

584 ax1.yaxis.grid(True, which="major", linestyle="--", linewidth=0.5, alpha=0.7)

585 ax1.xaxis.grid(True, which="major", linestyle="--", linewidth=0.5, alpha=0.7)

586

587 # Update legend with larger font size

588 ax1.legend(loc="lower left", fontsize=14)

589

590 # Calculate current seeing and age of data

591 if addMostRecentBox:

592 currentSeeing = seeings[-1]

593 justTime = currentSeeing.timestamp.isot.split("T")[1].split(".")[0]

594

595 text = f'Current Seeing: {currentSeeing.seeing:.2f}"\n' f"Last updated @ {justTime} UTC"

596 ax1.text(

597 0.05,

598 0.95,

599 text,

600 transform=ax1.transAxes,

601 fontsize=14,

602 verticalalignment="top",

603 horizontalalignment="left",

604 bbox=dict(boxstyle="round", facecolor="white", alpha=0.7),

605 )

606

607 return fig

Coverage for python / lsst / summit / extras / ringssSeeing.py: 26%

215 statements