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

1# This file is part of summit_utils.

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

21from __future__ import annotations

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]

46from collections.abc import Mapping

47from dataclasses import dataclass

48from typing import TYPE_CHECKING, Any

50import astropy.units as u

51import numpy as np

52from astropy.coordinates import AltAz, SkyCoord

53from astropy.time import Time

55from lsst.afw import cameraGeom

56from lsst.afw.cameraGeom import Detector

57from lsst.afw.image import ImageF, VisitInfo

59# at top-level imports

60from lsst.daf.base import PropertyList

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

62from lsst.obs.base import createInitialSkyWcsFromBoresight

63from lsst.obs.lsst import LsstCam

64from lsst.obs.lsst.cameraTransforms import LsstCameraTransforms

65from lsst.obs.lsst.translators.lsst import SIMONYI_LOCATION

67if TYPE_CHECKING:

68 from lsst.geom import SkyWcs

70 from .reading import GuiderData

72# build integer rotation matrices

73ROTATION_MATRICES = {

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

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

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

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

78}

79# build inverse rotation matrices

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

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

84 """

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

87 Parameters

88 ----------

89 wcs : `SkyWcs`

90 World Coordinate System object with pixelToSkyArray method.

91 x_flat : `np.ndarray`

92 Flattened array of x pixel coordinates.

93 y_flat : `np.ndarray`

94 Flattened array of y pixel coordinates.

96 Returns

97 -------

98 raFlat : `np.ndarray`

99 Array of right ascension values in radians.

100 decFlat : `np.ndarray`

101 Array of declination values in radians.

102 """

103 return wcs.pixelToSkyArray(x_flat, y_flat)

104

105

106def convertPixelsToAltaz(

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

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

109 """

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

111 in a vectorized manner.

112

113 Parameters

114 ----------

115 wcs :lsst.afw.geom.SkyWcs

116 World Coordinate System object with pixelToSkyArray method.

117 time : astropy.time.Time

118 Observation time.

119 xPix : np.ndarray

120 Array of x pixel coordinates.

121 yPix : np.ndarray

122 Array of y pixel coordinates.

123

124 Returns

125 -------

126 az : np.ndarray

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

128 alt : np.ndarray

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

130 """

131

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

133 x_arr = np.asarray(xPix)

134 y_arr = np.asarray(yPix)

135 shp = x_arr.shape

136

137 # 2) flatten for the WCS call

138 x_flat = x_arr.ravel()

139 y_flat = y_arr.ravel()

140

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

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

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

144

145 # 4) assemble a single SkyCoord and transform once

146 sc_icrs = SkyCoord(

147 ra=ra_flat * u.rad,

148 dec=dec_flat * u.rad,

149 frame="icrs",

150 obstime=time,

151 location=SIMONYI_LOCATION,

152 )

153 # Transform to AltAz

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

155

156 # 5) reshape back to original grid

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

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

159 return az, alt

160

161

162def convertFocalToAltaz(

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

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

165 """

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

167 chaining focal → pixel → AltAz transformations.

168

169 Parameters

170 ----------

171 wcs : lsst.afw.geom.SkyWcs

172 World Coordinate System object with pixelToSkyArray method.

173 time : astropy.time.Time

174 Observation time.

175 detector : lsst.afw.cameraGeom.Detector

176 Detector object.

177 xFocal : np.ndarray

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

179 yFocal : np.ndarray

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

181

182 Returns

183 -------

184 az : np.ndarray

185 Array of azimuth values in degrees.

186 alt : np.ndarray

187 Array of altitude values in degrees.

188 """

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

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

191

192

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

194# sky coordinates.

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

196 """

197 Create an AffineTransform representing a rotation by multiples of 90

198 degrees.

199

200 Parameters

201 ----------

202 detNquarter : `int`

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

204 direction : `int`, optional

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

206

207 Returns

208 -------

209 rotation : `AffineTransform`

210 The rotation transformation.

211

212 Raises

213 ------

214 ValueError

215 If direction is not 1 or -1.

216 """

217 rot = ROTATION_MATRICES

218 irot = INVERSE_ROTATION_MATRICES

219

220 nq = np.mod(detNquarter, 4)

221 if direction == 1:

222 rotation = AffineTransform(rot[nq])

223 elif direction == -1:

224 rotation = AffineTransform(irot[nq])

225 else:

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

227 return rotation

228

229

230def makeCcdToDvcsTransform(

231 bboxCcd: Box2I,

232 detNquarter: int,

233) -> tuple[AffineTransform, AffineTransform]:

234 """

235 Create forward and backward AffineTransforms for converting between CCD

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

237

238 Parameters

239 ----------

240 bboxCcd : `Box2I`

241 Bounding box for the ROI in CCD pixel coordinates.

242 detNquarter : `int`

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

244

245 Returns

246 -------

247 forwards : `AffineTransform`

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

249 coordinates (DVCS view).

250 backwards : `AffineTransform`

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

252 coordinates (CCD view).

253

254 Notes

255 -----

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

257 of the focal plane as used in DVCS.

258 """

259 # Use shared integer rotation matrices

260 nq = np.mod(detNquarter, 4)

261 rot = ROTATION_MATRICES

262 irot = INVERSE_ROTATION_MATRICES

263 # number of 90deg CCW rotations

264 nq = np.mod(detNquarter, 4)

265

266 # get LL,Size of the CCD view BBox

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

268 nx, ny = bboxCcd.getDimensions()

269

270 # get translations to use for each NQ value

271 boxtranslation = {}

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

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

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

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

276

277 # build the forware Transform

278 forwardTranslation = AffineTransform.makeTranslation(-lowerLeftCcd)

279 forwardRotation = AffineTransform(rot[nq])

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

281 # ordering is third*second*first

282 forwards = forwardBoxTranslation * forwardRotation * forwardTranslation

283

284 backwardTranslation = AffineTransform.makeTranslation(lowerLeftCcd)

285 backwardRotation = AffineTransform(irot[nq])

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

287 backwards = backwardTranslation * backwardRotation * backwardBoxTranslation

288

289 return forwards, backwards

290

291

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

293 """

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

295 """

296 # ROI metadata

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

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

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

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

301

302 # detector and amp

303 detector, ampName = getAmpNameFromMetadata(md, camera)

304

305 # corner 0 (CCD view)

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

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

308 corner0CcdX: int = int(corner0CcdXFloat)

309 corner0CcdY: int = int(corner0CcdYFloat)

310

311 # corner 2 (depends on flips)

312 amp = detector[ampName]

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

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

315

316 # bbox with LL/UR ordering

317 llX: int = min(corner0CcdX, corner2CcdX)

318 urX: int = max(corner0CcdX, corner2CcdX)

319 llY: int = min(corner0CcdY, corner2CcdY)

320 urY: int = max(corner0CcdY, corner2CcdY)

321

322 llCcd: Point2D = Point2D(llX, llY)

323 urCcd: Point2D = Point2D(urX, urY)

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

325 return ccdViewBbox

326

327

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

329 """

330 Unpack detector and amplifier information from metadata.

331

332 Parameters

333 ----------

334 md : lsst.daf.base.PropertyList

335 Metadata from one Guider CCD.

336 camera : lsst.obs.lsst.LsstCam

337 LsstCam object.

338

339 Returns

340 -------

341 detector : lsst.afw.cameraGeom.Detector

342 CCD Detector object.

343 ampName : str

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

345 """

346 raftBay: str = md["RAFTBAY"]

347 ccdSlot: str = md["CCDSLOT"]

348

349 segment: str = md["ROISEG"]

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

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

352 detector: Detector = camera[detName]

353

354 return detector, ampName

355

356

357def roiImageToDvcs(

358 roi: np.ndarray,

359 metadata: PropertyList,

360 detector: Detector,

361 ampName: str,

362 camera: LsstCam,

363 view: str = "dvcs",

364) -> ImageF:

365 """

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

367 or embed in full CCD image.

368

369 Parameters

370 ----------

371 roi : np.ndarray

372 ROI array (raw, amplifier coordinates).

373 metadata : lsst.daf.base.PropertyList

374 Metadata from one Guider CCD.

375 detector : lsst.afw.cameraGeom.Detector

376 CCD detector object.

377 ampName : str

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

379 camera : lsst.obs.lsst.LsstCam

380 LsstCam object.

381 view : str, optional

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

383

384 Returns

385 -------

386 imf : lsst.afw.image.ImageF

387 Image in the requested view.

388

389 Notes

390 -----

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

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

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

394 """

395

396 # convert image to ccd view

397 roiCcdView = ampToCcdView(roi, detector, ampName)

398

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

400 # afwMath.rotateImageBy90)

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

402

403 ny, nx = roi.shape

404 imf = ImageF(nx, ny, 0.0)

405

406 # output ImageF

407 if view == "dvcs":

408 imf.array[:] = roiDvcsView

409 elif view == "ccd":

410 imf.array[:] = roiCcdView

411 elif view == "ccdfull":

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

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

414 imf = ImageF(nx, ny, 0.0)

415 roiSegment = metadata["ROISEG"]

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

417 roiColumn = metadata["ROICOL"]

418 roiRow = metadata["ROIROW"]

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

420

421 return imf

422

423

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

425 """

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

427

428 Parameters

429 ----------

430 focalX : np.ndarray

431 Focal plane x coordinates (mm).

432 focalY : np.ndarray

433 Focal plane y coordinates (mm).

434 det : lsst.afw.cameraGeom.Detector

435 Detector of interest.

436

437 Returns

438 -------

439 x : np.ndarray

440 Pixel x coordinates.

441 y : np.ndarray

442 Pixel y coordinates.

443 """

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

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

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

447

448

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

450 """

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

452

453 Parameters

454 ----------

455 pixelX : np.ndarray

456 Pixel x coordinates.

457 pixelY : np.ndarray

458 Pixel y coordinates.

459 det : lsst.afw.cameraGeom.Detector

460 Detector of interest.

461

462 Returns

463 -------

464 focalX : np.ndarray

465 Focal plane x coordinates (mm).

466 focalY : np.ndarray

467 Focal plane y coordinates (mm).

468 """

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

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

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

472

473

474def stampToCcd(

475 stamp: np.ndarray,

476 ccdimg: np.ndarray,

477 detector: Detector,

478 camera: cameraGeom.Camera,

479 ampName: str,

480 ampCol: int,

481 ampRow: int,

482) -> np.ndarray:

483 """

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

485

486 Parameters

487 ----------

488 stamp : np.ndarray

489 Raw ROI array in amplifier coordinates.

490 ccdimg : np.ndarray

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

492 detector : lsst.afw.cameraGeom.Detector

493 Detector of interest.

494 camera : lsst.afw.cameraGeom.Camera

495 Camera object.

496 ampName : str

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

498 ampCol : int

499 Starting column number for ROI.

500 ampRow : int

501 Starting row number for ROI.

502

503 Returns

504 -------

505 ccdimg : np.ndarray

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

507

508 Notes

509 -----

510 This function may print debugging information if placement fails.

511 """

512

513 stamp_ccd = ampToCcdView(stamp, detector, ampName)

514 ampRows, ampCols = stamp_ccd.shape

515

516 # get corner0 of the location for the ROI

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

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

519 corner0_CCDX = int(corner0_CCDX)

520 corner0_CCDY = int(corner0_CCDY)

521

522 # get opposite corner, corner2, of the ROI

523 amp = detector[ampName]

524 if amp.getRawFlipX():

525 corner2_CCDX = corner0_CCDX - ampCols

526 else:

527 corner2_CCDX = corner0_CCDX + ampCols

528

529 if amp.getRawFlipY():

530 corner2_CCDY = corner0_CCDY - ampRows

531 else:

532 corner2_CCDY = corner0_CCDY + ampRows

533

534 # now place the ROI

535 yStart = min(corner0_CCDY, corner2_CCDY)

536 yEnd = max(corner0_CCDY, corner2_CCDY)

537 xStart = min(corner0_CCDX, corner2_CCDX)

538 xEnd = max(corner0_CCDX, corner2_CCDX)

539

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

541 raise ValueError(

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

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

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

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

546 )

547

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

549

550 return ccdimg

551

552

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

554 """

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

556

557 Parameters

558 ----------

559 stamp : np.ndarray

560 Raw ROI array in amplifier coordinates.

561 detector : lsst.afw.cameraGeom.Detector

562 Detector of interest.

563 ampName : str

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

565

566 Returns

567 -------

568 ccdImage : np.ndarray

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

570 """

571 amplifier = detector[ampName]

572 ccdImage = stamp.copy()

573 if amplifier.getRawFlipX():

574 ccdImage = np.fliplr(ccdImage)

575 if amplifier.getRawFlipY():

576 ccdImage = np.flipud(ccdImage)

577 return ccdImage

578

579

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

581 """

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

583

584 Parameters

585 ----------

586 xccd : np.ndarray

587 Array of x CCD pixel coordinates.

588 yccd : np.ndarray

589 Array of y CCD pixel coordinates.

590 detNum : int

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

592

593 Returns

594 -------

595 xfp : np.ndarray

596 Focal plane x coordinates (mm).

597 yfp : np.ndarray

598 Focal plane y coordinates (mm).

599 """

600 if len(xccd) > 0:

601 detector = LsstCam.getCamera()[detNum]

602 # Convert the star positions to focal plane coordinates

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

604 else:

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

606 return xfp, yfp

607

608

609def convertToAltaz(

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

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

612 """

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

614

615 Parameters

616 ----------

617 xccd : np.ndarray

618 Array of x CCD pixel coordinates.

619 yccd : np.ndarray

620 Array of y CCD pixel coordinates.

621 wcs : lsst.afw.image.Wcs

622 WCS object for the guider detector.

623 obsTime : astropy.time.Time

624 Observation time.

625

626 Returns

627 -------

628 alt : np.ndarray

629 Array of altitude coordinates (degrees).

630 az : np.ndarray

631 Array of azimuth coordinates (degrees).

632 """

633 if len(xccd) > 0:

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

635 else:

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

637

638 return alt, az

639

640

641def convertRoiToCcd(

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

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

644 """

645 Convert ROI coordinates to CCD pixel coordinates.

646

647 Parameters

648 ----------

649 xroi : np.ndarray

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

651 yroi : np.ndarray

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

653 guiderData : GuiderData

654 GuiderData object containing the guider datasets and view information.

655 guiderName : str

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

657

658 Returns

659 -------

660 xccd : np.ndarray

661 Array of x CCD pixel coordinates.

662 yccd : np.ndarray

663 Array of y CCD pixel coordinates.

664

665 Raises

666 ------

667 ValueError

668 If the view in GuiderData is not supported.

669 """

670 view = guiderData.view

671 stamps = guiderData[guiderName]

672

673 xroi = np.asarray(xroi)

674 yroi = np.asarray(yroi)

675

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

677 if view == "ccd":

678 # convert roi coords to ccd coords

679 # by adding the lower left corner of the box

680 lower_left_corner = box.getMin()

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

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

683

684 elif view == "dvcs":

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

686 # roi2ccd is an AffineTransform that converts

687 # from roi coords to ccd coords

688 xccd, yccd = roi2ccd(xroi, yroi)

689 else:

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

691

692 return xccd, yccd

693

694

695def convertCcdToRoi(

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

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

698 """

699 Convert CCD pixel coordinates to ROI pixel coordinates.

700

701 Parameters

702 ----------

703 xccd : np.ndarray

704 Array of x CCD pixel coordinates.

705 yccd : np.ndarray

706 Array of y CCD pixel coordinates.

707 guiderData : GuiderData

708 GuiderData object containing the guider datasets and view information.

709 guiderName : str

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

711

712 Returns

713 -------

714 xroi : np.ndarray

715 Array of x ROI pixel coordinates.

716 yroi : np.ndarray

717 Array of y ROI pixel coordinates.

718

719 Raises

720 ------

721 ValueError

722 If the view in GuiderData is not supported.

723 """

724 view = guiderData.view

725 stamps = guiderData[guiderName]

726

727 if len(xccd) == 0:

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

729

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

731

732 if view == "ccd":

733 lower_left_corner = box.getMin()

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

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

736

737 elif view == "dvcs":

738 xroi, yroi = ccd2dvcs(xccd, yccd)

739

740 else:

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

742

743 return xroi, yroi

744

745

746def convertCcdToDvcs(

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

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

749 """

750 Convert CCD pixel coordinates to DVCS (focal plane)

751 coordinates (pixels or mm).

752

753 Parameters

754 ----------

755 xccd : np.ndarray

756 Array of x CCD pixel coordinates.

757 yccd : np.ndarray

758 Array of y CCD pixel coordinates.

759 guiderData : GuiderData

760 GuiderData object containing the guider datasets and view information.

761 guiderName : str

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

763

764 Returns

765 -------

766 xdvcs : np.ndarray

767 Array of x DVCS focal plane coordinates.

768 ydvcs : np.ndarray

769 Array of y DVCS focal plane coordinates.

770 """

771 stamps = guiderData[guiderName]

772

773 if len(xccd) == 0:

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

775

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

777 xdvcs, ydvcs = ccd2dvcs(xccd, yccd)

778 return xdvcs, ydvcs

779

780

781def getCamRotAngle(visitinfo: VisitInfo) -> float:

782 """

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

784

785 Parameters

786 ----------

787 visitinfo : lsst.afw.image.VisitInfo

788 Visit information from an exposure.

789

790 Returns

791 -------

792 float

793 Camera rotation angle in degrees, wrapped near zero.

794 """

795 camRot = (

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

797 )

798 camRotAngle = Angle(camRot, degrees)

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

800 return camRotAngleW.asDegrees()

801

802

803# Putting former guiderwcs here

804

805

806def makeInitGuiderWcs(camera: cameraGeom.Camera, visitInfo: VisitInfo) -> dict[str, SkyWcs]:

807 """

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

809

810 Parameters

811 ----------

812 camera : lsst.afw.cameraGeom.Camera

813 Camera object containing detectors.

814 visitInfo : lsst.afw.image.VisitInfo

815 Visit information from an exposure.

816

817 Returns

818 -------

819 dict[str, SkyWcs]

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

821 """

822 orientation = visitInfo.getBoresightRotAngle()

823 boresight = visitInfo.getBoresightRaDec()

824

825 # Get WCS for each CCD in camera

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

827

828 for det in camera:

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

830 # get WCS

831 args = boresight, orientation, det

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

833

834 return camWcs

835

836

837# Data class for drift results

838@dataclass

839class DriftResult:

840 """

841 Data class to store results of telescope drift calculations.

842

843 Attributes

844 ----------

845 ra_point : float

846 Telescope pointing right ascension in degrees.

847 dec_point : float

848 Telescope pointing declination in degrees.

849 ra_real : float

850 Actual right ascension in degrees.

851 dec_real : float

852 Actual declination in degrees.

853 delta_ra_arcsec : float

854 Pointing error in RA in arcseconds.

855 delta_dec_arcsec : float

856 Pointing error in Dec in arcseconds.

857 az_drift_arcsec : float

858 Azimuth drift in arcseconds.

859 el_drift_arcsec : float

860 Elevation drift in arcseconds.

861 total_drift_arcsec : float

862 Total drift in arcseconds.

863 el_start : float

864 Starting elevation in degrees.

865 pixel_offset : tuple of float

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

867 """

868

869 ra_point: float

870 dec_point: float

871 ra_real: float

872 dec_real: float

873 delta_ra_arcsec: float

874 delta_dec_arcsec: float

875 az_drift_arcsec: float

876 el_drift_arcsec: float

877 total_drift_arcsec: float

878 el_start: float

879 pixel_offset: tuple[float, float]

880

881 def __str__(self) -> str:

882 """

883 Return a formatted string summarizing the drift results.

884

885 The exposure id is not included: a DriftResult does not carry one.

886 Use `summary` to print the results under an exposure id header.

887

888 Returns

889 -------

890 str

891 Formatted summary string.

892 """

893 s = (

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

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

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

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

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

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

900 f"\n"

901 f"-------------------- Drift Summary --------------------\n"

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

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

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

905 )

906 return s

907

908 def summary(self, expid: int) -> None:

909 """

910 Print a summary of the drift results under an exposure id header.

911

912 Parameters

913 ----------

914 expid : int

915 Exposure identifier to head the summary with.

916 """

917 print(f"Exposure summary: {expid}")

918 print(self)

919

920

921def getObsAltAz(

922 ra: float,

923 dec: float,

924 pressure: u.Quantity,

925 hum: float,

926 temperature: u.Quantity,

927 wl: u.Quantity,

928 time: Time,

929) -> SkyCoord:

930 """

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

932

933 Parameters

934 ----------

935 ra : float

936 Right ascension in degrees.

937 dec : float

938 Declination in degrees.

939 pressure : astropy.units.Quantity

940 Atmospheric pressure with units.

941 hum : float

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

943 temperature : astropy.units.Quantity

944 Temperature with units.

945 wl : astropy.units.Quantity

946 Wavelength of observation with units.

947 time : astropy.time.Time

948 Observation time.

949

950 Returns

951 -------

952 astropy.coordinates.SkyCoord

953 Altitude and azimuth coordinates of the object.

954 """

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

956 altAz1 = AltAz(

957 obstime=time,

958 location=SIMONYI_LOCATION,

959 pressure=pressure,

960 temperature=temperature,

961 relative_humidity=hum,

962 obswl=wl,

963 )

964 obsAltAz1 = skyLocation.transform_to(altAz1)

965 return obsAltAz1

966

967

968def DeltaAltAz(

969 ra: float,

970 dec: float,

971 pressure: u.Quantity,

972 hum: float,

973 temperature: u.Quantity,

974 wl: u.Quantity,

975 time1: Time,

976 time2: Time,

977) -> list[float]:

978 """

979 Calculate the change in AltAz coordinates between two times.

980

981 Parameters

982 ----------

983 ra : float

984 Right ascension in degrees.

985 dec : float

986 Declination in degrees.

987 pressure : astropy.units.Quantity

988 Atmospheric pressure with units.

989 hum : float

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

991 temperature : astropy.units.Quantity

992 Temperature with units.

993 wl : astropy.units.Quantity

994 Wavelength of observation with units.

995 time1 : astropy.time.Time

996 Start time of observation.

997 time2 : astropy.time.Time

998 End time of observation.

999

1000 Returns

1001 -------

1002 list of float

1003 List containing azimuth change and elevation change in arcseconds.

1004 """

1005 obsAltAz1 = getObsAltAz(ra, dec, pressure, hum, temperature, wl, time1)

1006 obsAltAz2 = getObsAltAz(ra, dec, pressure, hum, temperature, wl, time2)

1007 # 1 is at the beginning of the exposure, 2 is at the end

1008 # el, az are the actual values, prime values reflect the pointing model

1009 # These are all in degrees

1010 el1 = obsAltAz1.alt.deg

1011 az1 = obsAltAz1.az.deg

1012 el2 = obsAltAz2.alt.deg

1013 az2 = obsAltAz2.az.deg

1014 # Change values are the change from the beginning

1015 # to the end of the exposure,

1016 # in arcseconds

1017 azChange = (az2 - az1) * 3600.0

1018 elChange = (el2 - el1) * 3600.0

1019 return [azChange, elChange]

1020

1021

1022def computePointModelDrift(metadata: Mapping[str, Any], cWcs: SkyWcs, rWcs: SkyWcs) -> DriftResult:

1023 """

1024 Calculate Point Model drift from exposure metadata and

1025 the measured `calexp` WCS.

1026

1027 Parameters

1028 ----------

1029 metadata : dict-like

1030 Must contain at least: 'FILTBAND', 'PRESSURE', 'AIRTEMP', 'HUMIDITY',

1031 'MJD-BEG', 'MJD-END', 'RASTART', 'DECSTART', 'ELSTART'.

1032 cWcs : lsst.afw.geom.SkyWcs

1033 Calibrated WCS object.

1034 rWcs : lsst.afw.geom.SkyWcs

1035 Raw WCS object.

1036

1037 Returns

1038 -------

1039 DriftResult

1040 Object containing drift calculation results.

1041

1042 Example

1043 -------

1044 expId = 2025072300537

1045 rawExp = butler.get('raw', detector=94,

1046 exposure=expId, instrument=instrument)

1047 calExp = butler.get('preliminary_visit_image',

1048 detector=94, visit=expId, instrument=instrument)

1049

1050 md = rawExp.getMetadata()

1051 cWcs = calExp.getWcs()

1052 rWcs = rawExp.getWcs()

1053

1054 drift547 = calculate_drift(md, cWcs, rWcs)

1055 drift547.summary(2025072300537)

1056 """

1057 wavelengths = {"u": 3671, "g": 4827, "r": 6223, "i": 7546, "z": 8691, "y": 9712}

1058 filter = metadata["FILTBAND"]

1059 wl = wavelengths[filter] * u.angstrom

1060 pressure = metadata["PRESSURE"] * u.pascal

1061 temperature = metadata["AIRTEMP"] * u.Celsius

1062 hum = metadata["HUMIDITY"]

1063 time1 = Time(metadata["MJD-BEG"], format="mjd", scale="tai")

1064 time2 = Time(metadata["MJD-END"], format="mjd", scale="tai")

1065 ra_point = metadata["RASTART"]

1066 dec_point = metadata["DECSTART"]

1067 el_start = metadata["ELSTART"]

1068

1069 azChangePoint, elChangePoint = DeltaAltAz(

1070 ra_point, dec_point, pressure, hum, temperature, wl, time1, time2

1071 )

1072

1073 # Use WCS to get actual (ra, dec) at image center

1074 calExpSkyCenter = cWcs.pixelToSky(rWcs.getPixelOrigin())

1075 ra_real = calExpSkyCenter.getRa().asDegrees()

1076 dec_real = calExpSkyCenter.getDec().asDegrees()

1077 delta_ra = (ra_real - ra_point) * 3600.0

1078 delta_dec = (dec_real - dec_point) * 3600.0

1079

1080 # Calculate pixel offset from raw WCS origin to calibrated WCS origin

1081 pixel_offset = tuple(cWcs.getPixelOrigin() - rWcs.getPixelOrigin())

1082

1083 azChangeReal, elChangeReal = DeltaAltAz(ra_real, dec_real, pressure, hum, temperature, wl, time1, time2)

1084

1085 az_drift = azChangeReal - azChangePoint

1086 el_drift = elChangeReal - elChangePoint

1087 total_drift = np.sqrt(el_drift**2 + (az_drift * np.cos(np.radians(el_start))) ** 2)

1088

1089 return DriftResult(

1090 ra_point=ra_point,

1091 dec_point=dec_point,

1092 ra_real=ra_real,

1093 dec_real=dec_real,

1094 delta_ra_arcsec=delta_ra,

1095 delta_dec_arcsec=delta_dec,

1096 az_drift_arcsec=az_drift,

1097 el_drift_arcsec=el_drift,

1098 total_drift_arcsec=total_drift,

1099 el_start=el_start,

1100 pixel_offset=pixel_offset,

1101 )