Coverage for python / lsst / summit / utils / guiders / transformation.py: 18%

1# This file is part of summit_utils. 

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/>. 

21from __future__ import annotations 

22 

23__all__ = [ 

24 "makeRotationTransform", 

25 "makeCcdToDvcsTransform", 

26 "makeRoiBbox", 

27 "stampToCcd", 

28 "getAmpNameFromMetadata", 

29 "roiImageToDvcs", 

30 "focalToPixel", 

31 "pixelToFocal", 

32 "ampToCcdView", 

33 "convertFocalToAltaz", 

34 "convertPixelsToAltaz", 

35 "convertToFocalPlane", 

36 "convertRoiToCcd", 

37 "convertCcdToRoi", 

38 "convertCcdToDvcs", 

39 "convertToAltaz", 

40 "convertPixelToRadec", 

41 "makeInitGuiderWcs", 

42 "getCamRotAngle", 

43 "DriftResult", 

44] 

45 

46from dataclasses import dataclass 

47from typing import TYPE_CHECKING 

48 

49import astropy.units as u 

50import numpy as np 

51from astropy.coordinates import AltAz, SkyCoord 

52from astropy.time import Time 

53 

54from lsst.afw import cameraGeom 

55from lsst.afw.cameraGeom import Detector 

56from lsst.afw.image import ImageF 

57 

58# at top-level imports 

59from lsst.daf.base import PropertyList 

60from lsst.geom import AffineTransform, Angle, Box2D, Box2I, Extent2D, Point2D, degrees 

61from lsst.obs.base import createInitialSkyWcsFromBoresight 

62from lsst.obs.lsst import LsstCam 

63from lsst.obs.lsst.cameraTransforms import LsstCameraTransforms 

64from lsst.obs.lsst.translators.lsst import SIMONYI_LOCATION 

65 

66if TYPE_CHECKING: 

67 from lsst.geom import SkyWcs 

68 

69 from .reading import GuiderData 

70 

71# build integer rotation matrices 

72ROTATION_MATRICES = { 

73 0: np.array([[1.0, 0.0], [0.0, 1.0]]), 

74 1: np.array([[0.0, -1], [1.0, 0.0]]), 

75 2: np.array([[-1.0, 0.0], [0.0, -1.0]]), 

76 3: np.array([[0.0, 1.0], [-1.0, 0.0]]), 

77} 

78# build inverse rotation matrices 

79INVERSE_ROTATION_MATRICES = {k: v.transpose() for k, v in ROTATION_MATRICES.items()} 

80 

81 

82def convertPixelToRadec(wcs: SkyWcs, x_flat: np.ndarray, y_flat: np.ndarray) -> tuple[np.ndarray, np.ndarray]: 

83 """ 

84 Map detector-pixel coordinates to ICRS RA and DEC (in radians). 

85 

86 Parameters 

87 ---------- 

88 wcs : `SkyWcs` 

89 World Coordinate System object with pixelToSkyArray method. 

90 x_flat : `np.ndarray` 

91 Flattened array of x pixel coordinates. 

92 y_flat : `np.ndarray` 

93 Flattened array of y pixel coordinates. 

94 

95 Returns 

96 ------- 

97 raFlat : `np.ndarray` 

98 Array of right ascension values in radians. 

99 decFlat : `np.ndarray` 

100 Array of declination values in radians. 

101 """ 

102 return wcs.pixelToSkyArray(x_flat, y_flat) 

103 

104 

105def convertPixelsToAltaz( 

106 wcs: SkyWcs, time: Time, xPix: np.ndarray, yPix: np.ndarray 

107) -> tuple[np.ndarray, np.ndarray]: 

108 """ 

109 Convert detector-pixel coordinates to Altitude and Azimuth (degrees) 

110 in a vectorized manner. 

111 

112 Parameters 

113 ---------- 

114 wcs :lsst.afw.geom.SkyWcs 

115 World Coordinate System object with pixelToSkyArray method. 

116 time : astropy.time.Time 

117 Observation time. 

118 xPix : np.ndarray 

119 Array of x pixel coordinates. 

120 yPix : np.ndarray 

121 Array of y pixel coordinates. 

122 

123 Returns 

124 ------- 

125 az : np.ndarray 

126 Array of azimuth values in degrees (same shape as xPix/yPix). 

127 alt : np.ndarray 

128 Array of altitude values in degrees (same shape as xPix/yPix). 

129 """ 

130 

131 # 1) make sure we have numpy arrays, remember their shape 

132 x_arr = np.asarray(xPix) 

133 y_arr = np.asarray(yPix) 

134 shp = x_arr.shape 

135 

136 # 2) flatten for the WCS call 

137 x_flat = x_arr.ravel() 

138 y_flat = y_arr.ravel() 

139 

140 # 3) scalarSKyWcs → ICRS RA/Dec (radians) in bulk: 

141 # (pixelToSkyArray expects floats and returns two 1D arrays) 

142 ra_flat, dec_flat = convertPixelToRadec(wcs, x_flat, y_flat) 

143 

144 # 4) assemble a single SkyCoord and transform once 

145 sc_icrs = SkyCoord( 

146 ra=ra_flat * u.rad, 

147 dec=dec_flat * u.rad, 

148 frame="icrs", 

149 obstime=time, 

150 location=SIMONYI_LOCATION, 

151 ) 

152 # Transform to AltAz 

153 sc_altaz = sc_icrs.transform_to(AltAz(obstime=time, location=SIMONYI_LOCATION)) 

154 

155 # 5) reshape back to original grid 

156 az = sc_altaz.az.deg.reshape(shp) 

157 alt = sc_altaz.alt.deg.reshape(shp) 

158 return az, alt 

159 

160 

161def convertFocalToAltaz( 

162 wcs: SkyWcs, time: Time, detector: Detector, xFocal: np.ndarray, yFocal: np.ndarray 

163) -> tuple[np.ndarray, np.ndarray]: 

164 """ 

165 Map focal-plane coordinates (mm) to Altitude/Azimuth (degrees) by 

166 chaining focal → pixel → AltAz transformations. 

167 

168 Parameters 

169 ---------- 

170 wcs : lsst.afw.geom.SkyWcs 

171 World Coordinate System object with pixelToSkyArray method. 

172 time : astropy.time.Time 

173 Observation time. 

174 detector : lsst.afw.cameraGeom.Detector 

175 Detector object. 

176 xFocal : np.ndarray 

177 Array of x focal-plane coordinates (mm). 

178 yFocal : np.ndarray 

179 Array of y focal-plane coordinates (mm). 

180 

181 Returns 

182 ------- 

183 az : np.ndarray 

184 Array of azimuth values in degrees. 

185 alt : np.ndarray 

186 Array of altitude values in degrees. 

187 """ 

188 x_pix, y_pix = focalToPixel(xFocal, yFocal, detector) 

189 return convertPixelsToAltaz(wcs, time, x_pix, y_pix) 

190 

191 

192# Aaron's code to make transformations from ROI coordinates to Focal Plane and 

193# sky coordinates. 

194def makeRotationTransform(detNquarter: int, direction: int = 1) -> AffineTransform: 

195 """ 

196 Create an AffineTransform representing a rotation by multiples of 90 

197 degrees. 

198 

199 Parameters 

200 ---------- 

201 detNquarter : `int` 

202 Number of 90 degree counter-clockwise rotations (modulo 4). 

203 direction : `int`, optional 

204 1 for forward rotation (default), -1 for reverse rotation (transpose) 

205 

206 Returns 

207 ------- 

208 rotation : `AffineTransform` 

209 The rotation transformation. 

210 

211 Raises 

212 ------ 

213 ValueError 

214 If direction is not 1 or -1. 

215 """ 

216 rot = ROTATION_MATRICES 

217 irot = INVERSE_ROTATION_MATRICES 

218 

219 nq = np.mod(detNquarter, 4) 

220 if direction == 1: 

221 rotation = AffineTransform(rot[nq]) 

222 elif direction == -1: 

223 rotation = AffineTransform(irot[nq]) 

224 else: 

225 raise ValueError(f"direction must be either +/- 1, got {direction}") 

226 return rotation 

227 

228 

229def makeCcdToDvcsTransform( 

230 bboxCcd: Box2I, 

231 detNquarter: int, 

232) -> tuple[AffineTransform, AffineTransform]: 

233 """ 

234 Create forward and backward AffineTransforms for converting between CCD 

235 pixel coordinates and DVCS (focal plane) stamp coordinates. 

236 

237 Parameters 

238 ---------- 

239 bboxCcd : `Box2I` 

240 Bounding box for the ROI in CCD pixel coordinates. 

241 detNquarter : `int` 

242 Number of 90-degree CCW rotations to align CCD view to DVCS view. 

243 

244 Returns 

245 ------- 

246 forwards : `AffineTransform` 

247 Transform from CCD pixel coordinates (CCD view) to stamp pixel 

248 coordinates (DVCS view). 

249 backwards : `AffineTransform` 

250 Transform from stamp pixel coordinates (DVCS view) to CCD pixel 

251 coordinates (CCD view). 

252 

253 Notes 

254 ----- 

255 Useful for mapping between the raw detector coordinates and the orientation 

256 of the focal plane as used in DVCS. 

257 """ 

258 # Use shared integer rotation matrices 

259 nq = np.mod(detNquarter, 4) 

260 rot = ROTATION_MATRICES 

261 irot = INVERSE_ROTATION_MATRICES 

262 # number of 90deg CCW rotations 

263 nq = np.mod(detNquarter, 4) 

264 

265 # get LL,Size of the CCD view BBox 

266 lowerLeftCcd = Extent2D(bboxCcd.getCorners()[0]) 

267 nx, ny = bboxCcd.getDimensions() 

268 

269 # get translations to use for each NQ value 

270 boxtranslation = {} 

271 boxtranslation[0] = Extent2D(0, 0) 

272 boxtranslation[1] = Extent2D(ny - 1, 0) 

273 boxtranslation[2] = Extent2D(nx - 1, ny - 1) 

274 boxtranslation[3] = Extent2D(0, nx - 1) 

275 

276 # build the forware Transform 

277 forwardTranslation = AffineTransform.makeTranslation(-lowerLeftCcd) 

278 forwardRotation = AffineTransform(rot[nq]) 

279 forwardBoxTranslation = AffineTransform.makeTranslation(boxtranslation[nq]) 

280 # ordering is third*second*first 

281 forwards = forwardBoxTranslation * forwardRotation * forwardTranslation 

282 

283 backwardTranslation = AffineTransform.makeTranslation(lowerLeftCcd) 

284 backwardRotation = AffineTransform(irot[nq]) 

285 backwardBoxTranslation = AffineTransform.makeTranslation(-boxtranslation[nq]) 

286 backwards = backwardTranslation * backwardRotation * backwardBoxTranslation 

287 

288 return forwards, backwards 

289 

290 

291def makeRoiBbox(md: PropertyList, camera: LsstCam) -> Box2I: 

292 """ 

293 Construct the bounding box for a Guider ROI stamp in CCD coords. 

294 """ 

295 # ROI metadata 

296 roiColumn: int = int(md["ROICOL"]) 

297 roiRow: int = int(md["ROIROW"]) 

298 roiColumns: int = int(md["ROICOLS"]) 

299 roiRows: int = int(md["ROIROWS"]) 

300 

301 # detector and amp 

302 detector, ampName = getAmpNameFromMetadata(md, camera) 

303 

304 # corner 0 (CCD view) 

305 lct = LsstCameraTransforms(camera, detector.getName()) 

306 corner0CcdXFloat, corner0CcdYFloat = lct.ampPixelToCcdPixel(roiColumn, roiRow, ampName) 

307 corner0CcdX: int = int(corner0CcdXFloat) 

308 corner0CcdY: int = int(corner0CcdYFloat) 

309 

310 # corner 2 (depends on flips) 

311 amp = detector[ampName] 

312 corner2CcdX: int = corner0CcdX - (roiColumns - 1) if amp.getRawFlipX() else corner0CcdX + (roiColumns - 1) 

313 corner2CcdY: int = corner0CcdY - (roiRows - 1) if amp.getRawFlipY() else corner0CcdY + (roiRows - 1) 

314 

315 # bbox with LL/UR ordering 

316 llX: int = min(corner0CcdX, corner2CcdX) 

317 urX: int = max(corner0CcdX, corner2CcdX) 

318 llY: int = min(corner0CcdY, corner2CcdY) 

319 urY: int = max(corner0CcdY, corner2CcdY) 

320 

321 llCcd: Point2D = Point2D(llX, llY) 

322 urCcd: Point2D = Point2D(urX, urY) 

323 ccdViewBbox: Box2I = Box2I(Box2D(llCcd, urCcd)) 

324 return ccdViewBbox 

325 

326 

327def getAmpNameFromMetadata(md: PropertyList, camera: LsstCam) -> tuple[Detector, str]: 

328 """ 

329 Unpack detector and amplifier information from metadata. 

330 

331 Parameters 

332 ---------- 

333 md : lsst.daf.base.PropertyList 

334 Metadata from one Guider CCD. 

335 camera : lsst.obs.lsst.LsstCam 

336 LsstCam object. 

337 

338 Returns 

339 ------- 

340 detector : lsst.afw.cameraGeom.Detector 

341 CCD Detector object. 

342 ampName : str 

343 Amplifier name, e.g., 'C00'. 

344 """ 

345 raftBay: str = md["RAFTBAY"] 

346 ccdSlot: str = md["CCDSLOT"] 

347 

348 segment: str = md["ROISEG"] 

349 ampName: str = "C" + segment[7:] 

350 detName: str = f"{raftBay}_{ccdSlot}" 

351 detector: Detector = camera[detName] 

352 

353 return detector, ampName 

354 

355 

356def roiImageToDvcs( 

357 roi: np.ndarray, 

358 metadata: PropertyList, 

359 detector: Detector, 

360 ampName: str, 

361 camera: LsstCam, 

362 view: str = "dvcs", 

363) -> ImageF: 

364 """ 

365 Convert ROI image from raw amplifier view to CCD or DVCS views, 

366 or embed in full CCD image. 

367 

368 Parameters 

369 ---------- 

370 roi : np.ndarray 

371 ROI array (raw, amplifier coordinates). 

372 metadata : lsst.daf.base.PropertyList 

373 Metadata from one Guider CCD. 

374 detector : lsst.afw.cameraGeom.Detector 

375 CCD detector object. 

376 ampName : str 

377 Amplifier name, e.g., 'C00'. 

378 camera : lsst.obs.lsst.LsstCam 

379 LsstCam object. 

380 view : str, optional 

381 Desired output view: 'dvcs' (default), 'ccd', or 'ccdfull'. 

382 

383 Returns 

384 ------- 

385 imf : lsst.afw.image.ImageF 

386 Image in the requested view. 

387 

388 Notes 

389 ----- 

390 - 'dvcs': Oriented as in the focal plane. 

391 - 'ccd': Oriented as in CCD coordinates. 

392 - 'ccdfull': ROI embedded in a full CCD-sized image. 

393 """ 

394 

395 # convert image to ccd view 

396 roiCcdView = ampToCcdView(roi, detector, ampName) 

397 

398 # convert image to dvcs view (nb. np.rot90 works opposite to 

399 # afwMath.rotateImageBy90) 

400 roiDvcsView = np.rot90(roiCcdView, -detector.getOrientation().getNQuarter()) 

401 

402 ny, nx = roi.shape 

403 imf = ImageF(nx, ny, 0.0) 

404 

405 # output ImageF 

406 if view == "dvcs": 

407 imf.array[:] = roiDvcsView 

408 elif view == "ccd": 

409 imf.array[:] = roiCcdView 

410 elif view == "ccdfull": 

411 # make the full CCD and place the ROI inside it 

412 nx, ny = detector.getBBox().getDimensions() 

413 imf = ImageF(nx, ny, 0.0) 

414 roiSegment = metadata["ROISEG"] 

415 ampName = "C" + roiSegment[7:] 

416 roiColumn = metadata["ROICOL"] 

417 roiRow = metadata["ROIROW"] 

418 imf.array[:] = stampToCcd(roi, imf.array[:], detector, camera, ampName, roiColumn, roiRow) 

419 

420 return imf 

421 

422 

423def focalToPixel(focalX: np.ndarray, focalY: np.ndarray, det: Detector) -> tuple[np.ndarray, np.ndarray]: 

424 """ 

425 Convert focal plane coordinates (mm, DVCS) to detector pixel coordinates. 

426 

427 Parameters 

428 ---------- 

429 focalX : np.ndarray 

430 Focal plane x coordinates (mm). 

431 focalY : np.ndarray 

432 Focal plane y coordinates (mm). 

433 det : lsst.afw.cameraGeom.Detector 

434 Detector of interest. 

435 

436 Returns 

437 ------- 

438 x : np.ndarray 

439 Pixel x coordinates. 

440 y : np.ndarray 

441 Pixel y coordinates. 

442 """ 

443 transform = det.getTransform(cameraGeom.FOCAL_PLANE, cameraGeom.PIXELS) 

444 x, y = transform.getMapping().applyForward(np.vstack((focalX, focalY))) 

445 return x.ravel(), y.ravel() 

446 

447 

448def pixelToFocal(pixelX: np.ndarray, pixelY: np.ndarray, det: Detector) -> tuple[np.ndarray, np.ndarray]: 

449 """ 

450 Convert pixel coordinates to focal plane coordinates (mm, DVCS). 

451 

452 Parameters 

453 ---------- 

454 pixelX : np.ndarray 

455 Pixel x coordinates. 

456 pixelY : np.ndarray 

457 Pixel y coordinates. 

458 det : lsst.afw.cameraGeom.Detector 

459 Detector of interest. 

460 

461 Returns 

462 ------- 

463 focalX : np.ndarray 

464 Focal plane x coordinates (mm). 

465 focalY : np.ndarray 

466 Focal plane y coordinates (mm). 

467 """ 

468 transform = det.getTransform(cameraGeom.PIXELS, cameraGeom.FOCAL_PLANE) 

469 focalX, focalY = transform.getMapping().applyForward(np.vstack((pixelX, pixelY))) 

470 return focalX.ravel(), focalY.ravel() 

471 

472 

473def stampToCcd(stamp, ccdimg, detector, camera, ampName, ampCol, ampRow): 

474 """ 

475 Place ROI (stamp) into an existing full CCD array, all in CCD view. 

476 

477 Parameters 

478 ---------- 

479 stamp : np.ndarray 

480 Raw ROI array in amplifier coordinates. 

481 ccdimg : np.ndarray 

482 CCD image array in CCD coordinates (to be filled). 

483 detector : lsst.afw.cameraGeom.Detector 

484 Detector of interest. 

485 camera : lsst.afw.cameraGeom.Camera 

486 Camera object. 

487 ampName : str 

488 Amplifier name, e.g., 'C00'. 

489 ampCol : int 

490 Starting column number for ROI. 

491 ampRow : int 

492 Starting row number for ROI. 

493 

494 Returns 

495 ------- 

496 ccdimg : np.ndarray 

497 CCD image array with the ROI inserted in the appropriate location. 

498 

499 Notes 

500 ----- 

501 This function may print debugging information if placement fails. 

502 """ 

503 

504 stamp_ccd = ampToCcdView(stamp, detector, ampName) 

505 ampRows, ampCols = stamp_ccd.shape 

506 

507 # get corner0 of the location for the ROI 

508 lct = LsstCameraTransforms(camera, detector.getName()) 

509 corner0_CCDX, corner0_CCDY = lct.ampPixelToCcdPixel(ampCol, ampRow, ampName) 

510 corner0_CCDX = int(corner0_CCDX) 

511 corner0_CCDY = int(corner0_CCDY) 

512 

513 # get opposite corner, corner2, of the ROI 

514 amp = detector[ampName] 

515 if amp.getRawFlipX(): 

516 corner2_CCDX = corner0_CCDX - ampCols 

517 else: 

518 corner2_CCDX = corner0_CCDX + ampCols 

519 

520 if amp.getRawFlipY(): 

521 corner2_CCDY = corner0_CCDY - ampRows 

522 else: 

523 corner2_CCDY = corner0_CCDY + ampRows 

524 

525 # now place the ROI 

526 yStart = min(corner0_CCDY, corner2_CCDY) 

527 yEnd = max(corner0_CCDY, corner2_CCDY) 

528 xStart = min(corner0_CCDX, corner2_CCDX) 

529 xEnd = max(corner0_CCDX, corner2_CCDX) 

530 

531 if stamp_ccd.shape != (yEnd - yStart, xEnd - xStart): 

532 raise ValueError( 

533 f"ROI shape mismatch for detector '{detector.getName()}', amp '{ampName}'.\n" 

534 f"Expected shape: {(yEnd - yStart, xEnd - xStart)}, " 

535 f"Got: {stamp_ccd.shape}\n" 

536 f"ROI box X: {xStart}:{xEnd}, Y: {yStart}:{yEnd}" 

537 ) 

538 

539 ccdimg[yStart:yEnd, xStart:xEnd] = stamp_ccd 

540 

541 return ccdimg 

542 

543 

544def ampToCcdView(stamp: np.ndarray, detector: Detector, ampName: str) -> np.ndarray: 

545 """ 

546 Convert a Guider ROI stamp image from amplifier view to CCD view. 

547 

548 Parameters 

549 ---------- 

550 stamp : np.ndarray 

551 Raw ROI array in amplifier coordinates. 

552 detector : lsst.afw.cameraGeom.Detector 

553 Detector of interest. 

554 ampName : str 

555 Amplifier name, e.g., 'C00'. 

556 

557 Returns 

558 ------- 

559 ccdImage : np.ndarray 

560 ROI image array flipped to be in CCD coordinate orientation. 

561 """ 

562 amplifier = detector[ampName] 

563 ccdImage = stamp.copy() 

564 if amplifier.getRawFlipX(): 

565 ccdImage = np.fliplr(ccdImage) 

566 if amplifier.getRawFlipY(): 

567 ccdImage = np.flipud(ccdImage) 

568 return ccdImage 

569 

570 

571def convertToFocalPlane(xccd: np.ndarray, yccd: np.ndarray, detNum: int) -> tuple[np.ndarray, np.ndarray]: 

572 """ 

573 Convert from CCD pixel coordinates to focal plane coordinates (mm). 

574 

575 Parameters 

576 ---------- 

577 xccd : np.ndarray 

578 Array of x CCD pixel coordinates. 

579 yccd : np.ndarray 

580 Array of y CCD pixel coordinates. 

581 detNum : int 

582 Detector number (e.g., 189 for R22_S11). 

583 

584 Returns 

585 ------- 

586 xfp : np.ndarray 

587 Focal plane x coordinates (mm). 

588 yfp : np.ndarray 

589 Focal plane y coordinates (mm). 

590 """ 

591 if len(xccd) > 0: 

592 detector = LsstCam.getCamera()[detNum] 

593 # Convert the star positions to focal plane coordinates 

594 xfp, yfp = pixelToFocal(xccd, yccd, detector) 

595 else: 

596 xfp, yfp = np.array([]), np.array([]) 

597 return xfp, yfp 

598 

599 

600def convertToAltaz( 

601 xccd: np.ndarray, yccd: np.ndarray, wcs: SkyWcs, obsTime: Time 

602) -> tuple[np.ndarray, np.ndarray]: 

603 """ 

604 Convert CCD pixel coordinates to altitude and azimuth (degrees). 

605 

606 Parameters 

607 ---------- 

608 xccd : np.ndarray 

609 Array of x CCD pixel coordinates. 

610 yccd : np.ndarray 

611 Array of y CCD pixel coordinates. 

612 wcs : lsst.afw.image.Wcs 

613 WCS object for the guider detector. 

614 obsTime : astropy.time.Time 

615 Observation time. 

616 

617 Returns 

618 ------- 

619 alt : np.ndarray 

620 Array of altitude coordinates (degrees). 

621 az : np.ndarray 

622 Array of azimuth coordinates (degrees). 

623 """ 

624 if len(xccd) > 0: 

625 az, alt = convertPixelsToAltaz(wcs, obsTime, xccd, yccd) 

626 else: 

627 alt, az = np.array([]), np.array([]) 

628 

629 return alt, az 

630 

631 

632def convertRoiToCcd( 

633 xroi: np.ndarray, yroi: np.ndarray, guiderData: GuiderData, guiderName: str 

634) -> tuple[np.ndarray, np.ndarray]: 

635 """ 

636 Convert ROI coordinates to CCD pixel coordinates. 

637 

638 Parameters 

639 ---------- 

640 xroi : np.ndarray 

641 Array of x ROI coordinates (pixels within the ROI). 

642 yroi : np.ndarray 

643 Array of y ROI coordinates (pixels within the ROI). 

644 guiderData : GuiderData 

645 GuiderData object containing the guider datasets and view information. 

646 guiderName : str 

647 Name of the guider (e.g., 'R44_SG0'). 

648 

649 Returns 

650 ------- 

651 xccd : np.ndarray 

652 Array of x CCD pixel coordinates. 

653 yccd : np.ndarray 

654 Array of y CCD pixel coordinates. 

655 

656 Raises 

657 ------ 

658 ValueError 

659 If the view in GuiderData is not supported. 

660 """ 

661 view = guiderData.view 

662 stamps = guiderData[guiderName] 

663 

664 xroi = np.asarray(xroi) 

665 yroi = np.asarray(yroi) 

666 

667 box, _, roi2ccd = stamps.getArchiveElements()[0] 

668 if view == "ccd": 

669 # convert roi coords to ccd coords 

670 # by adding the lower left corner of the box 

671 lower_left_corner = box.getMin() 

672 xmin, ymin = lower_left_corner.getX(), lower_left_corner.getY() 

673 xccd, yccd = xroi + xmin, yroi + ymin 

674 

675 elif view == "dvcs": 

676 # convert roi coords to ccd coords using the roi2ccd transform 

677 # roi2ccd is an AffineTransform that converts 

678 # from roi coords to ccd coords 

679 xccd, yccd = roi2ccd(xroi, yroi) 

680 else: 

681 raise ValueError(f"Unsupported view '{view}' in convertRoiToCcd" "must be 'ccd' or 'dvcs'") 

682 

683 return xccd, yccd 

684 

685 

686def convertCcdToRoi( 

687 xccd: np.ndarray, yccd: np.ndarray, guiderData: GuiderData, guiderName: str 

688) -> tuple[np.ndarray, np.ndarray]: 

689 """ 

690 Convert CCD pixel coordinates to ROI pixel coordinates. 

691 

692 Parameters 

693 ---------- 

694 xccd : np.ndarray 

695 Array of x CCD pixel coordinates. 

696 yccd : np.ndarray 

697 Array of y CCD pixel coordinates. 

698 guiderData : GuiderData 

699 GuiderData object containing the guider datasets and view information. 

700 guiderName : str 

701 Name of the guider (e.g., 'R44_SG0'). 

702 

703 Returns 

704 ------- 

705 xroi : np.ndarray 

706 Array of x ROI pixel coordinates. 

707 yroi : np.ndarray 

708 Array of y ROI pixel coordinates. 

709 

710 Raises 

711 ------ 

712 ValueError 

713 If the view in GuiderData is not supported. 

714 """ 

715 view = guiderData.view 

716 stamps = guiderData[guiderName] 

717 

718 if len(xccd) == 0: 

719 return np.array([]), np.array([]) 

720 

721 box, ccd2dvcs, _ = stamps.getArchiveElements()[0] 

722 

723 if view == "ccd": 

724 lower_left_corner = box.getMin() 

725 xmin, ymin = lower_left_corner.getX(), lower_left_corner.getY() 

726 xroi, yroi = xccd - xmin, yccd - ymin 

727 

728 elif view == "dvcs": 

729 xroi, yroi = ccd2dvcs(xccd, yccd) 

730 

731 else: 

732 raise ValueError(f"Unsupported view '{view}' in convertCcdToRoi", "must be 'ccd' or 'dvcs'") 

733 

734 return xroi, yroi 

735 

736 

737def convertCcdToDvcs( 

738 xccd: np.ndarray, yccd: np.ndarray, guiderData: GuiderData, guiderName: str 

739) -> tuple[np.ndarray, np.ndarray]: 

740 """ 

741 Convert CCD pixel coordinates to DVCS (focal plane) 

742 coordinates (pixels or mm). 

743 

744 Parameters 

745 ---------- 

746 xccd : np.ndarray 

747 Array of x CCD pixel coordinates. 

748 yccd : np.ndarray 

749 Array of y CCD pixel coordinates. 

750 guiderData : GuiderData 

751 GuiderData object containing the guider datasets and view information. 

752 guiderName : str 

753 Name of the guider (e.g., 'R44_SG0'). 

754 

755 Returns 

756 ------- 

757 xdvcs : np.ndarray 

758 Array of x DVCS focal plane coordinates. 

759 ydvcs : np.ndarray 

760 Array of y DVCS focal plane coordinates. 

761 """ 

762 stamps = guiderData[guiderName] 

763 

764 if len(xccd) == 0: 

765 return np.array([]), np.array([]) 

766 

767 _, ccd2dvcs, _ = stamps.getArchiveElements()[0] 

768 xdvcs, ydvcs = ccd2dvcs(xccd, yccd) 

769 return xdvcs, ydvcs 

770 

771 

772def getCamRotAngle(visitinfo) -> float: 

773 """ 

774 Compute the camera rotation angle in degrees from visit information. 

775 

776 Parameters 

777 ---------- 

778 visitinfo : lsst.afw.image.VisitInfo 

779 Visit information from an exposure. 

780 

781 Returns 

782 ------- 

783 float 

784 Camera rotation angle in degrees, wrapped near zero. 

785 """ 

786 camRot = ( 

787 visitinfo.getBoresightParAngle().asDegrees() - visitinfo.getBoresightRotAngle().asDegrees() - 90.0 

788 ) 

789 camRotAngle = Angle(camRot, degrees) 

790 camRotAngleW = camRotAngle.wrapNear(Angle(0.0)) 

791 return camRotAngleW.asDegrees() 

792 

793 

794# Putting former guiderwcs here 

795 

796 

797def makeInitGuiderWcs(camera, visitInfo) -> dict[str, SkyWcs]: 

798 """ 

799 Create initial WCS for each guider detector in the camera. 

800 

801 Parameters 

802 ---------- 

803 camera : lsst.afw.cameraGeom.Camera 

804 Camera object containing detectors. 

805 visitInfo : lsst.afw.image.VisitInfo 

806 Visit information from an exposure. 

807 

808 Returns 

809 ------- 

810 dict[str, SkyWcs] 

811 Dictionary of WCS objects keyed by detector name for guider detectors. 

812 """ 

813 orientation = visitInfo.getBoresightRotAngle() 

814 boresight = visitInfo.getBoresightRaDec() 

815 

816 # Get WCS for each CCD in camera 

817 camWcs: dict[str, SkyWcs] = {} 

818 

819 for det in camera: 

820 if det.getType() == cameraGeom.DetectorType.GUIDER: 

821 # get WCS 

822 args = boresight, orientation, det 

823 camWcs[det.getName()] = createInitialSkyWcsFromBoresight(*args, flipX=False) 

824 

825 return camWcs 

826 

827 

828# Data class for drift results 

829@dataclass 

830class DriftResult: 

831 """ 

832 Data class to store results of telescope drift calculations. 

833 

834 Attributes 

835 ---------- 

836 ra_point : float 

837 Telescope pointing right ascension in degrees. 

838 dec_point : float 

839 Telescope pointing declination in degrees. 

840 ra_real : float 

841 Actual right ascension in degrees. 

842 dec_real : float 

843 Actual declination in degrees. 

844 delta_ra_arcsec : float 

845 Pointing error in RA in arcseconds. 

846 delta_dec_arcsec : float 

847 Pointing error in Dec in arcseconds. 

848 az_drift_arcsec : float 

849 Azimuth drift in arcseconds. 

850 el_drift_arcsec : float 

851 Elevation drift in arcseconds. 

852 total_drift_arcsec : float 

853 Total drift in arcseconds. 

854 el_start : float 

855 Starting elevation in degrees. 

856 pixel_offset : tuple of float 

857 Pixel offset from raw WCS origin to calibrated WCS origin. 

858 """ 

859 

860 ra_point: float 

861 dec_point: float 

862 ra_real: float 

863 dec_real: float 

864 delta_ra_arcsec: float 

865 delta_dec_arcsec: float 

866 az_drift_arcsec: float 

867 el_drift_arcsec: float 

868 total_drift_arcsec: float 

869 el_start: float 

870 pixel_offset: tuple[float, float] 

871 

872 def __str__(self, expid=""): 

873 """ 

874 Return a formatted string summarizing the drift results. 

875 

876 Parameters 

877 ---------- 

878 expid : str, optional 

879 Exposure identifier to include in the summary (default is ""). 

880 

881 Returns 

882 ------- 

883 str 

884 Formatted summary string. 

885 """ 

886 s = ( 

887 f"Exposure summary: {expid}\n" 

888 f"Telescope pointing (RA, Dec): ({self.ra_point:.6f}, {self.dec_point:.6f})\n" 

889 f"Actual pointing (RA, Dec) : ({self.ra_real:.6f}, {self.dec_real:.6f})\n" 

890 f"Pointing error (RA, Dec). : ({self.delta_ra_arcsec:.1f}," 

891 f"{self.delta_dec_arcsec:.1f}) arcseconds\n" 

892 f"Starting elevation. : {self.el_start:.2f} degrees\n" 

893 f"Pixel offset from origin : ({self.pixel_offset[0]:.4f}, {self.pixel_offset[1]:.4f}) pixels\n" 

894 f"\n" 

895 f"-------------------- Drift Summary --------------------\n" 

896 f"Azimuth drift : {self.az_drift_arcsec:.2f} arcseconds\n" 

897 f"Elevation drift : {self.el_drift_arcsec:.2f} arcseconds\n" 

898 f"Total drift. : {self.total_drift_arcsec:.2f} arcseconds\n" 

899 ) 

900 return s 

901 

902 def summary(self, expid=""): 

903 """ 

904 Print a summary of the drift results. 

905 

906 Parameters 

907 ---------- 

908 expid : str, optional 

909 Exposure identifier to include in the summary (default is ""). 

910 """ 

911 print(self.__str__(expid)) 

912 

913 

914def getObsAltAz(ra, dec, pressure, hum, temperature, wl, time): 

915 """ 

916 Compute the observed AltAz coordinates for given RA/Dec and conditions. 

917 

918 Parameters 

919 ---------- 

920 ra : float 

921 Right ascension in degrees. 

922 dec : float 

923 Declination in degrees. 

924 pressure : astropy.units.Quantity 

925 Atmospheric pressure with units. 

926 hum : float 

927 Relative humidity as a fraction (0-1). 

928 temperature : astropy.units.Quantity 

929 Temperature with units. 

930 wl : astropy.units.Quantity 

931 Wavelength of observation with units. 

932 time : astropy.time.Time 

933 Observation time. 

934 

935 Returns 

936 ------- 

937 astropy.coordinates.AltAz 

938 Altitude and azimuth coordinates of the object. 

939 """ 

940 skyLocation = SkyCoord(ra * u.deg, dec * u.deg) 

941 altAz1 = AltAz( 

942 obstime=time, 

943 location=SIMONYI_LOCATION, 

944 pressure=pressure, 

945 temperature=temperature, 

946 relative_humidity=hum, 

947 obswl=wl,