From 0bb107d2c4da13192d19779e77c2e7381168f694 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 14 Nov 2025 15:06:18 -0800 Subject: [PATCH] Update the parallactic angle calculation for DCR --- python/lsst/ip/diffim/dcrModel.py | 16 +++----- tests/test_dcrModel.py | 63 ++++++++++--------------------- 2 files changed, 25 insertions(+), 54 deletions(-) diff --git a/python/lsst/ip/diffim/dcrModel.py b/python/lsst/ip/diffim/dcrModel.py index ea47a9a96..f55e2ff1b 100644 --- a/python/lsst/ip/diffim/dcrModel.py +++ b/python/lsst/ip/diffim/dcrModel.py @@ -804,7 +804,7 @@ def calculateDcr(visitInfo, wcs, effectiveWavelength, bandwidth, dcrNumSubfilter The 2D shift due to DCR, in pixels. Uses numpy axes ordering (Y, X). """ - rotation = calculateImageParallacticAngle(visitInfo, wcs) + rotation = calculateImageParallacticAngle(visitInfo) dcrShift = [] weight = [0.75, 0.25] for wl0, wl1 in wavelengthGenerator(effectiveWavelength, bandwidth, dcrNumSubfilters): @@ -840,15 +840,16 @@ def calculateDcr(visitInfo, wcs, effectiveWavelength, bandwidth, dcrNumSubfilter return dcrShift -def calculateImageParallacticAngle(visitInfo, wcs): +def calculateImageParallacticAngle(visitInfo, wcs=None): """Calculate the total sky rotation angle of an exposure. Parameters ---------- visitInfo : `lsst.afw.image.VisitInfo` Metadata for the exposure. - wcs : `lsst.afw.geom.SkyWcs` + wcs : `lsst.afw.geom.SkyWcs`, optional Coordinate system definition (wcs) for the exposure. + Deprecated, will be removed following v30. Returns ------- @@ -861,14 +862,7 @@ def calculateImageParallacticAngle(visitInfo, wcs): A rotation angle of 90 degrees is defined with North along the +x axis and East along the -y axis. """ - parAngle = visitInfo.getBoresightParAngle().asRadians() - cd = wcs.getCdMatrix() - if wcs.isFlipped: - cdAngle = (np.arctan2(-cd[0, 1], cd[0, 0]) + np.arctan2(cd[1, 0], cd[1, 1]))/2. - rotAngle = (cdAngle + parAngle)*geom.radians - else: - cdAngle = (np.arctan2(cd[0, 1], -cd[0, 0]) + np.arctan2(cd[1, 0], cd[1, 1]))/2. - rotAngle = (cdAngle - parAngle)*geom.radians + rotAngle = visitInfo.boresightParAngle - visitInfo.boresightRotAngle return rotAngle diff --git a/tests/test_dcrModel.py b/tests/test_dcrModel.py index 030d2d44c..2fd4e76f8 100644 --- a/tests/test_dcrModel.py +++ b/tests/test_dcrModel.py @@ -174,7 +174,7 @@ def makeDummyWcs(self, rotAngle, pixelScale, crval, flipX=True): wcs = afwGeom.makeSkyWcs(crpix=crpix, crval=crval, cdMatrix=cdMatrix) return wcs - def makeDummyVisitInfo(self, azimuth, elevation, exposureId=12345, randomizeTime=False): + def makeDummyVisitInfo(self, azimuth, elevation, exposureId=12345, randomizeTime=False, rotAngle=0.): """Make a self-consistent visitInfo object for testing. Parameters @@ -215,7 +215,7 @@ def makeDummyVisitInfo(self, azimuth, elevation, exposureId=12345, randomizeTime datetime_begin=time, datetime_end=time, boresight_airmass=airmass, - boresight_rotation_angle=Angle(0.*u.degree), + boresight_rotation_angle=Angle(rotAngle*u.degree), boresight_rotation_coord='sky', temperature=lsstTemperature, pressure=lsstPressure, @@ -261,7 +261,7 @@ def testDcrCalculation(self): azimuth = 30.*degrees elevation = 65.*degrees pixelScale = 0.2*arcseconds - visitInfo = self.makeDummyVisitInfo(azimuth, elevation) + visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees()) wcs = self.makeDummyWcs(rotAngle, pixelScale, crval=visitInfo.getBoresightRaDec()) dcrShift = calculateDcr(visitInfo, wcs, self.effectiveWavelength, self.bandwidth, self.dcrNumSubfilters) @@ -282,15 +282,13 @@ def testCoordinateTransformDcrCalculation(self): amplitude to the altitude, and transform back to pixel coordinates. """ pixelScale = 0.2*arcseconds - doFlip = [False, True] for testIter in range(self.nRandIter): rotAngle = 360.*self.rng.rand()*degrees azimuth = 360.*self.rng.rand()*degrees elevation = (45. + self.rng.rand()*40.)*degrees # Restrict to 45 < elevation < 85 degrees - visitInfo = self.makeDummyVisitInfo(azimuth, elevation) - for flip in doFlip: - # Repeat the calculation for both WCS orientations + visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees()) + for flip in [True, False]: wcs = self.makeDummyWcs(rotAngle, pixelScale, crval=visitInfo.getBoresightRaDec(), flipX=flip) dcrShifts = calculateDcr(visitInfo, wcs, self.effectiveWavelength, self.bandwidth, @@ -300,7 +298,10 @@ def testCoordinateTransformDcrCalculation(self): self.dcrNumSubfilters) for refShift, dcrShift in zip(refShifts, dcrShifts): # Use a fairly loose tolerance, since 1% of a pixel is good enough agreement. - self.assertFloatsAlmostEqual(refShift[1], dcrShift[1], rtol=1e-2, atol=1e-2) + if wcs.isFlipped: + self.assertFloatsAlmostEqual(refShift[1], dcrShift[1], rtol=1e-2, atol=1e-2) + else: + self.assertFloatsAlmostEqual(-refShift[1], dcrShift[1], rtol=1e-2, atol=1e-2) self.assertFloatsAlmostEqual(refShift[0], dcrShift[0], rtol=1e-2, atol=1e-2) def testDcrSubfilterOrder(self): @@ -311,12 +312,12 @@ def testDcrSubfilterOrder(self): rotAngle = 360.*self.rng.rand()*degrees azimuth = 360.*self.rng.rand()*degrees elevation = (45. + self.rng.rand()*40.)*degrees # Restrict to 45 < elevation < 85 degrees - visitInfo = self.makeDummyVisitInfo(azimuth, elevation) + visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees()) wcs = self.makeDummyWcs(rotAngle, pixelScale, crval=visitInfo.getBoresightRaDec()) dcrShift = calculateDcr(visitInfo, wcs, self.effectiveWavelength, self.bandwidth, self.dcrNumSubfilters) # First check that the blue subfilter amplitude is greater than the red subfilter - rotation = calculateImageParallacticAngle(visitInfo, wcs).asRadians() + rotation = calculateImageParallacticAngle(visitInfo).asRadians() ampShift = [dcr[1]*np.sin(rotation) + dcr[0]*np.cos(rotation) for dcr in dcrShift] self.assertGreater(ampShift[0], 0.) # The blue subfilter should be shifted towards zenith self.assertLess(ampShift[2], 0.) # The red subfilter should be shifted away from zenith @@ -348,13 +349,11 @@ def testRotationAngle(self): The rotation is compared to pre-computed values. """ - cdRotAngle = 0.*degrees + rotAngle = 0.*degrees azimuth = 130.*degrees elevation = 70.*degrees - pixelScale = 0.2*arcseconds - visitInfo = self.makeDummyVisitInfo(azimuth, elevation) - wcs = self.makeDummyWcs(cdRotAngle, pixelScale, crval=visitInfo.getBoresightRaDec()) - rotAngle = calculateImageParallacticAngle(visitInfo, wcs) + visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees()) + rotAngle = calculateImageParallacticAngle(visitInfo) refAngle = -1.0848032636337364*radians self.assertAnglesAlmostEqual(refAngle, rotAngle) @@ -364,38 +363,16 @@ def testRotationSouthZero(self): An observation pointed South and on the meridian should have zenith directly to the North, and a parallactic angle of zero. """ - refAngle = 0.*degrees azimuth = 180.*degrees # Telescope is pointed South - pixelScale = 0.2*arcseconds - for testIter in range(self.nRandIter): - # Any additional arbitrary rotation should fall out of the calculation - cdRotAngle = 360*self.rng.rand()*degrees - elevation = (45. + self.rng.rand()*40.)*degrees # Restrict to 45 < elevation < 85 degrees - visitInfo = self.makeDummyVisitInfo(azimuth, elevation) - wcs = self.makeDummyWcs(cdRotAngle, pixelScale, crval=visitInfo.getBoresightRaDec(), flipX=True) - rotAngle = calculateImageParallacticAngle(visitInfo, wcs) - self.assertAnglesAlmostEqual(refAngle - cdRotAngle, rotAngle, maxDiff=coordinateTolerance) - - def testRotationFlipped(self): - """Check the interpretation of rotations in the WCS. - """ - doFlip = [False, True] for testIter in range(self.nRandIter): # Any additional arbitrary rotation should fall out of the calculation - cdRotAngle = 360*self.rng.rand()*degrees - # Make the telescope be pointed South, so that the parallactic angle is zero. - azimuth = 180.*degrees + rotAngle = 360*self.rng.rand()*degrees elevation = (45. + self.rng.rand()*40.)*degrees # Restrict to 45 < elevation < 85 degrees - pixelScale = 0.2*arcseconds - visitInfo = self.makeDummyVisitInfo(azimuth, elevation) - for flip in doFlip: - wcs = self.makeDummyWcs(cdRotAngle, pixelScale, - crval=visitInfo.getBoresightRaDec(), - flipX=flip) - rotAngle = calculateImageParallacticAngle(visitInfo, wcs) - if flip: - rotAngle *= -1 - self.assertAnglesAlmostEqual(cdRotAngle, rotAngle, maxDiff=coordinateTolerance) + visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees()) + parAngle = calculateImageParallacticAngle(visitInfo) + self.assertAnglesAlmostEqual(visitInfo.boresightParAngle, geom.Angle(0.), + maxDiff=coordinateTolerance) + self.assertAnglesAlmostEqual(rotAngle, -parAngle, maxDiff=coordinateTolerance) def testConditionDcrModelNoChange(self): """Conditioning should not change the model if it equals the reference.