Coverage for tests / test_dcrModel.py: 13%
290 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-01 08:37 +0000
1# This file is part of ip_diffim.
2#
3# LSST Data Management System
4# This product includes software developed by the
5# LSST Project (http://www.lsst.org/).
6# See COPYRIGHT file at the top of the source tree.
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <https://www.lsstcorp.org/LegalNotices/>.
22from astropy import units as u
23from astropy.coordinates import SkyCoord, EarthLocation, Angle
24from astropy.time import Time
25import numpy as np
26from scipy import ndimage
27import unittest
29from astro_metadata_translator import makeObservationInfo
30from lsst.afw.coord import differentialRefraction
31import lsst.afw.geom as afwGeom
32import lsst.afw.image as afwImage
33import lsst.afw.math as afwMath
34import lsst.geom as geom
35from lsst.geom import arcseconds, degrees, radians, arcminutes
36from lsst.ip.diffim.dcrModel import (DcrModel, calculateDcr, calculateImageParallacticAngle,
37 applyDcr, wavelengthGenerator)
38from lsst.obs.base import MakeRawVisitInfoViaObsInfo
39from lsst.meas.algorithms.testUtils import plantSources
40import lsst.utils.tests
43# Our calculation of hour angle and parallactic angle ignore precession
44# and nutation, so calculations depending on these are not precise. DM-20133
45coordinateTolerance = 1.*arcminutes
48class DcrModelTestTask(lsst.utils.tests.TestCase):
49 """A test case for the DCR-aware image coaddition algorithm.
51 Attributes
52 ----------
53 bbox : `lsst.afw.geom.Box2I`
54 Bounding box of the test model.
55 bufferSize : `int`
56 Distance from the inner edge of the bounding box
57 to avoid placing test sources in the model images.
58 dcrNumSubfilters : int
59 Number of sub-filters used to model chromatic effects within a band.
60 effectiveWavelength : `float`
61 Effective wavelength of the full band.
62 lambdaMax : `float`
63 Maximum wavelength where the relative throughput
64 of the band is greater than 1%.
65 lambdaMin : `float`
66 Minimum wavelength where the relative throughput
67 of the band is greater than 1%.
68 mask : `lsst.afw.image.Mask`
69 Reference mask of the unshifted model.
70 """
72 def setUp(self):
73 """Define the filter, DCR parameters, and the bounding box for the tests.
74 """
75 self.rng = np.random.RandomState(5)
76 self.nRandIter = 10 # Number of iterations to repeat each test with random numbers.
77 self.dcrNumSubfilters = 3
78 self.effectiveWavelength = 476.31 # Use LSST g band values for the test.
79 self.bandwidth = 552. - 405.
80 self.bufferSize = 5
81 xSize = 40
82 ySize = 42
83 x0 = 12345
84 y0 = 67890
85 self.bbox = geom.Box2I(geom.Point2I(x0, y0), geom.Extent2I(xSize, ySize))
87 def makeTestImages(self, seed=5, nSrc=5, psfSize=2., noiseLevel=5.,
88 detectionSigma=5., sourceSigma=20., fluxRange=2.):
89 """Make reproduceable PSF-convolved masked images for testing.
91 Parameters
92 ----------
93 seed : `int`, optional
94 Seed value to initialize the random number generator.
95 nSrc : `int`, optional
96 Number of sources to simulate.
97 psfSize : `float`, optional
98 Width of the PSF of the simulated sources, in pixels.
99 noiseLevel : `float`, optional
100 Standard deviation of the noise to add to each pixel.
101 detectionSigma : `float`, optional
102 Threshold amplitude of the image to set the "DETECTED" mask.
103 sourceSigma : `float`, optional
104 Average amplitude of the simulated sources,
105 relative to ``noiseLevel``
106 fluxRange : `float`, optional
107 Range in flux amplitude of the simulated sources.
109 Returns
110 -------
111 modelImages : `list` of `lsst.afw.image.Image`
112 A list of images, each containing the model for one subfilter
113 """
114 rng = np.random.RandomState(seed)
115 x0, y0 = self.bbox.getBegin()
116 xSize, ySize = self.bbox.getDimensions()
117 xLoc = rng.rand(nSrc)*(xSize - 2*self.bufferSize) + self.bufferSize + x0
118 yLoc = rng.rand(nSrc)*(ySize - 2*self.bufferSize) + self.bufferSize + y0
119 modelImages = []
121 imageSum = np.zeros((ySize, xSize))
122 for subfilter in range(self.dcrNumSubfilters):
123 flux = (rng.rand(nSrc)*(fluxRange - 1.) + 1.)*sourceSigma*noiseLevel
124 sigmas = [psfSize for src in range(nSrc)]
125 coordList = list(zip(xLoc, yLoc, flux, sigmas))
126 model = plantSources(self.bbox, 10, 0, coordList, addPoissonNoise=False)
127 model.image.array += rng.rand(ySize, xSize)*noiseLevel
128 imageSum += model.image.array
129 model.mask.addMaskPlane("CLIPPED")
130 modelImages.append(model.image)
131 maskVals = np.zeros_like(imageSum)
132 maskVals[imageSum > detectionSigma*noiseLevel] = afwImage.Mask.getPlaneBitMask('DETECTED')
133 model.mask.array[:] = maskVals
134 self.mask = model.mask
135 return modelImages
137 def prepareStats(self):
138 """Make a simple statistics object for testing.
140 Returns
141 -------
142 statsCtrl : `lsst.afw.math.StatisticsControl`
143 Statistics control object for coaddition.
144 """
145 statsCtrl = afwMath.StatisticsControl()
146 statsCtrl.setNumSigmaClip(5)
147 statsCtrl.setNumIter(3)
148 statsCtrl.setNanSafe(True)
149 statsCtrl.setWeighted(True)
150 statsCtrl.setCalcErrorFromInputVariance(False)
151 return statsCtrl
153 def makeDummyWcs(self, rotAngle, pixelScale, crval, flipX=True):
154 """Make a World Coordinate System object for testing.
156 Parameters
157 ----------
158 rotAngle : `lsst.geom.Angle`
159 rotation of the CD matrix, East from North
160 pixelScale : `lsst.geom.Angle`
161 Pixel scale of the projection.
162 crval : `lsst.afw.geom.SpherePoint`
163 Coordinates of the reference pixel of the wcs.
164 flipX : `bool`, optional
165 Flip the direction of increasing Right Ascension.
167 Returns
168 -------
169 `lsst.afw.geom.skyWcs.SkyWcs`
170 A wcs that matches the inputs.
171 """
172 crpix = geom.Box2D(self.bbox).getCenter()
173 cdMatrix = afwGeom.makeCdMatrix(scale=pixelScale, orientation=rotAngle, flipX=flipX)
174 wcs = afwGeom.makeSkyWcs(crpix=crpix, crval=crval, cdMatrix=cdMatrix)
175 return wcs
177 def makeDummyVisitInfo(self, azimuth, elevation, exposureId=12345, randomizeTime=False, rotAngle=0.):
178 """Make a self-consistent visitInfo object for testing.
180 Parameters
181 ----------
182 azimuth : `lsst.geom.Angle`
183 Azimuth angle of the simulated observation.
184 elevation : `lsst.geom.Angle`
185 Elevation angle of the simulated observation.
186 exposureId : `int`, optional
187 Unique integer identifier for this observation (detector+visit id).
188 randomizeTime : `bool`, optional
189 Add a random offset to the observation time.
191 Returns
192 -------
193 `lsst.afw.image.VisitInfo`
194 VisitInfo for the exposure.
195 """
196 lsstLat = -30.244639*u.degree
197 lsstLon = -70.749417*u.degree
198 lsstAlt = 2663.*u.m
199 lsstTemperature = 20.*u.Celsius
200 lsstHumidity = 40. # in percent
201 lsstPressure = 73892.*u.pascal
202 loc = EarthLocation(lat=lsstLat,
203 lon=lsstLon,
204 height=lsstAlt)
205 airmass = 1.0/np.sin(elevation.asDegrees())
207 time = Time(2000.0, format="jyear", scale="tt")
208 if randomizeTime:
209 # Pick a random date and time within a 20-year span
210 time += 20*u.year*self.rng.rand()
211 altaz = SkyCoord(alt=elevation.asDegrees(), az=azimuth.asDegrees(),
212 unit='deg', obstime=time, frame='altaz', location=loc)
213 obsInfo = makeObservationInfo(location=loc,
214 detector_exposure_id=exposureId,
215 datetime_begin=time,
216 datetime_end=time,
217 boresight_airmass=airmass,
218 boresight_rotation_angle=Angle(rotAngle*u.degree),
219 boresight_rotation_coord='sky',
220 temperature=lsstTemperature,
221 pressure=lsstPressure,
222 relative_humidity=lsstHumidity,
223 tracking_radec=altaz.icrs,
224 altaz_begin=altaz,
225 observation_type='science',
226 )
227 visitInfo = MakeRawVisitInfoViaObsInfo.observationInfo2visitInfo(obsInfo)
228 return visitInfo
230 def testDummyVisitInfo(self):
231 """Verify the implementation of the visitInfo used for tests.
232 """
233 azimuth = 0*degrees
234 for testIter in range(self.nRandIter):
235 # Restrict to 45 < elevation < 85 degrees
236 elevation = (45. + self.rng.rand()*40.)*degrees
237 visitInfo = self.makeDummyVisitInfo(azimuth, elevation)
238 dec = visitInfo.getBoresightRaDec().getLatitude()
239 lat = visitInfo.getObservatory().getLatitude()
240 # An observation made with azimuth=0 should be pointed to the North
241 # So the RA should be North of the telescope's latitude
242 self.assertGreater(dec.asDegrees(), lat.asDegrees())
244 # The hour angle should be zero for azimuth=0
245 HA = visitInfo.getBoresightHourAngle()
246 refHA = 0.*degrees
247 self.assertAnglesAlmostEqual(HA, refHA, maxDiff=coordinateTolerance)
248 # If the observation is North of the telescope's latitude, the
249 # direction to zenith should be along the -y axis
250 # with a parallactic angle of 180 degrees
251 parAngle = visitInfo.getBoresightParAngle()
252 refParAngle = 180.*degrees
253 self.assertAnglesAlmostEqual(parAngle, refParAngle, maxDiff=coordinateTolerance)
255 def testDcrCalculation(self):
256 """Test that the shift in pixels due to DCR is consistently computed.
258 The shift is compared to pre-computed values.
259 """
260 rotAngle = 0.*degrees
261 azimuth = 30.*degrees
262 elevation = 65.*degrees
263 pixelScale = 0.2*arcseconds
264 visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees())
265 wcs = self.makeDummyWcs(rotAngle, pixelScale, crval=visitInfo.getBoresightRaDec())
266 dcrShift = calculateDcr(visitInfo, wcs, self.effectiveWavelength,
267 self.bandwidth, self.dcrNumSubfilters)
268 # Compare to precomputed values.
269 refShift = [(-0.587550080368824, -0.28495575189083666),
270 (-0.019158809951774006, -0.009291825969480522),
271 (0.38855779160585996, 0.18844653649028056)]
272 for shiftOld, shiftNew in zip(refShift, dcrShift):
273 self.assertFloatsAlmostEqual(shiftOld[1], shiftNew[1], rtol=1e-6, atol=1e-8)
274 self.assertFloatsAlmostEqual(shiftOld[0], shiftNew[0], rtol=1e-6, atol=1e-8)
276 def testCoordinateTransformDcrCalculation(self):
277 """Check the DCR calculation using astropy coordinate transformations.
279 Astmospheric refraction causes sources to appear closer to zenith than
280 they really are. An alternate calculation of the shift due to DCR is to
281 transform the pixel coordinates to altitude and azimuth, add the DCR
282 amplitude to the altitude, and transform back to pixel coordinates.
283 """
284 pixelScale = 0.2*arcseconds
286 for testIter in range(self.nRandIter):
287 rotAngle = 360.*self.rng.rand()*degrees
288 azimuth = 360.*self.rng.rand()*degrees
289 elevation = (45. + self.rng.rand()*40.)*degrees # Restrict to 45 < elevation < 85 degrees
290 visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees())
291 for flip in [True, False]:
292 wcs = self.makeDummyWcs(rotAngle, pixelScale, crval=visitInfo.getBoresightRaDec(), flipX=flip)
293 dcrShifts = calculateDcr(visitInfo, wcs,
294 self.effectiveWavelength, self.bandwidth,
295 self.dcrNumSubfilters)
296 refShifts = calculateAstropyDcr(visitInfo, wcs,
297 self.effectiveWavelength, self.bandwidth,
298 self.dcrNumSubfilters)
299 for refShift, dcrShift in zip(refShifts, dcrShifts):
300 # Use a fairly loose tolerance, since 1% of a pixel is good enough agreement.
301 if wcs.isFlipped:
302 self.assertFloatsAlmostEqual(refShift[1], dcrShift[1], rtol=1e-2, atol=1e-2)
303 else:
304 self.assertFloatsAlmostEqual(-refShift[1], dcrShift[1], rtol=1e-2, atol=1e-2)
305 self.assertFloatsAlmostEqual(refShift[0], dcrShift[0], rtol=1e-2, atol=1e-2)
307 def testDcrSubfilterOrder(self):
308 """Test that the bluest subfilter always has the largest DCR amplitude.
309 """
310 pixelScale = 0.2*arcseconds
311 for testIter in range(self.nRandIter):
312 rotAngle = 360.*self.rng.rand()*degrees
313 azimuth = 360.*self.rng.rand()*degrees
314 elevation = (45. + self.rng.rand()*40.)*degrees # Restrict to 45 < elevation < 85 degrees
315 visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees())
316 wcs = self.makeDummyWcs(rotAngle, pixelScale, crval=visitInfo.getBoresightRaDec())
317 dcrShift = calculateDcr(visitInfo, wcs, self.effectiveWavelength, self.bandwidth,
318 self.dcrNumSubfilters)
319 # First check that the blue subfilter amplitude is greater than the red subfilter
320 rotation = calculateImageParallacticAngle(visitInfo).asRadians()
321 ampShift = [dcr[1]*np.sin(rotation) + dcr[0]*np.cos(rotation) for dcr in dcrShift]
322 self.assertGreater(ampShift[0], 0.) # The blue subfilter should be shifted towards zenith
323 self.assertLess(ampShift[2], 0.) # The red subfilter should be shifted away from zenith
324 # The absolute amplitude of the blue subfilter should also
325 # be greater than that of the red subfilter
326 self.assertGreater(np.abs(ampShift[0]), np.abs(ampShift[2]))
328 def testApplyDcr(self):
329 """Test that the image transformation reduces to a simple shift.
330 """
331 dxVals = [-2, 1, 0, 1, 2]
332 dyVals = [-2, 1, 0, 1, 2]
333 x0 = 13
334 y0 = 27
335 inputImage = afwImage.MaskedImageF(self.bbox)
336 image = inputImage.image.array
337 image[y0, x0] = 1.
338 for dx in dxVals:
339 for dy in dyVals:
340 shift = (dy, dx)
341 shiftedImage = applyDcr(image, shift, useInverse=False)
342 # Create a blank reference image, and add the fake point source at the shifted location.
343 refImage = afwImage.MaskedImageF(self.bbox)
344 refImage.image.array[y0 + dy, x0 + dx] = 1.
345 self.assertFloatsAlmostEqual(shiftedImage, refImage.image.array, rtol=1e-12, atol=1e-12)
347 def testRotationAngle(self):
348 """Test that the sky rotation angle is consistently computed.
350 The rotation is compared to pre-computed values.
351 """
352 rotAngle = 0.*degrees
353 azimuth = 130.*degrees
354 elevation = 70.*degrees
355 visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees())
356 rotAngle = calculateImageParallacticAngle(visitInfo)
357 refAngle = -1.0848032636337364*radians
358 self.assertAnglesAlmostEqual(refAngle, rotAngle)
360 def testRotationSouthZero(self):
361 """Test that an observation pointed due South has zero rotation angle.
363 An observation pointed South and on the meridian should have zenith
364 directly to the North, and a parallactic angle of zero.
365 """
366 azimuth = 180.*degrees # Telescope is pointed South
367 for testIter in range(self.nRandIter):
368 # Any additional arbitrary rotation should fall out of the calculation
369 rotAngle = 360*self.rng.rand()*degrees
370 elevation = (45. + self.rng.rand()*40.)*degrees # Restrict to 45 < elevation < 85 degrees
371 visitInfo = self.makeDummyVisitInfo(azimuth, elevation, rotAngle=rotAngle.asDegrees())
372 parAngle = calculateImageParallacticAngle(visitInfo)
373 self.assertAnglesAlmostEqual(visitInfo.boresightParAngle, geom.Angle(0.),
374 maxDiff=coordinateTolerance)
375 self.assertAnglesAlmostEqual(rotAngle, -parAngle, maxDiff=coordinateTolerance)
377 def testConditionDcrModelNoChange(self):
378 """Conditioning should not change the model if it equals the reference.
379 """
380 modelImages = self.makeTestImages()
381 dcrModels = DcrModel(modelImages=modelImages, mask=self.mask,
382 effectiveWavelength=self.effectiveWavelength, bandwidth=self.bandwidth)
383 newModels = [model.clone() for model in dcrModels]
384 dcrModels.conditionDcrModel(newModels, self.bbox, gain=1.)
385 for refModel, newModel in zip(dcrModels, newModels):
386 self.assertFloatsAlmostEqual(refModel.array, newModel.array)
388 def testConditionDcrModelNoChangeHighGain(self):
389 """Conditioning should not change the model if it equals the reference.
390 """
391 modelImages = self.makeTestImages()
392 dcrModels = DcrModel(modelImages=modelImages, mask=self.mask,
393 effectiveWavelength=self.effectiveWavelength, bandwidth=self.bandwidth)
394 newModels = [model.clone() for model in dcrModels]
395 dcrModels.conditionDcrModel(newModels, self.bbox, gain=3.)
396 for refModel, newModel in zip(dcrModels, newModels):
397 self.assertFloatsAlmostEqual(refModel.array, newModel.array)
399 def testConditionDcrModelWithChange(self):
400 """Verify conditioning when the model changes by a known amount.
401 """
402 modelImages = self.makeTestImages()
403 dcrModels = DcrModel(modelImages=modelImages, mask=self.mask,
404 effectiveWavelength=self.effectiveWavelength, bandwidth=self.bandwidth)
405 newModels = [model.clone() for model in dcrModels]
406 for model in newModels:
407 model.array[:] *= 3.
408 dcrModels.conditionDcrModel(newModels, self.bbox, gain=1.)
409 for refModel, newModel in zip(dcrModels, newModels):
410 refModel.array[:] *= 2.
411 self.assertFloatsAlmostEqual(refModel.array, newModel.array)
413 def testRegularizationSmallClamp(self):
414 """Test that large variations between model planes are reduced.
416 This also tests that noise-like pixels are not regularized.
417 """
418 clampFrequency = 2
419 regularizationWidth = 2
420 fluxRange = 10.
421 modelImages = self.makeTestImages(fluxRange=fluxRange)
422 dcrModels = DcrModel(modelImages=modelImages, mask=self.mask,
423 effectiveWavelength=self.effectiveWavelength, bandwidth=self.bandwidth)
424 newModels = [model.clone() for model in dcrModels]
425 templateImage = dcrModels.getReferenceImage(self.bbox)
427 statsCtrl = self.prepareStats()
428 dcrModels.regularizeModelFreq(newModels, self.bbox, statsCtrl, clampFrequency, regularizationWidth)
429 for model, refModel in zip(newModels, dcrModels):
430 # Make sure the test parameters do reduce the outliers
431 self.assertGreater(np.max(refModel.array - templateImage),
432 np.max(model.array - templateImage))
433 highThreshold = templateImage*clampFrequency
434 highPix = model.array > highThreshold
435 highPix = ndimage.binary_opening(highPix, iterations=regularizationWidth)
436 self.assertFalse(np.all(highPix))
437 lowThreshold = templateImage/clampFrequency
438 lowPix = model.array < lowThreshold
439 lowPix = ndimage.binary_opening(lowPix, iterations=regularizationWidth)
440 self.assertFalse(np.all(lowPix))
442 def testRegularizationSidelobes(self):
443 """Test that artificial chromatic sidelobes are suppressed.
444 """
445 clampFrequency = 2.
446 regularizationWidth = 2
447 noiseLevel = 0.1
448 sourceAmplitude = 100.
449 modelImages = self.makeTestImages(seed=5, nSrc=5, psfSize=3., noiseLevel=noiseLevel,
450 detectionSigma=5., sourceSigma=sourceAmplitude, fluxRange=2.)
451 templateImage = np.mean([model.array for model in modelImages], axis=0)
452 sidelobeImages = self.makeTestImages(seed=5, nSrc=5, psfSize=1.5, noiseLevel=noiseLevel/10.,
453 detectionSigma=5., sourceSigma=sourceAmplitude*5., fluxRange=2.)
454 statsCtrl = self.prepareStats()
455 signList = [-1., 0., 1.]
456 sidelobeShift = (0., 4.)
457 for model, sidelobe, sign in zip(modelImages, sidelobeImages, signList):
458 sidelobe.array *= sign
459 model.array += applyDcr(sidelobe.array, sidelobeShift, useInverse=False)
460 model.array += applyDcr(sidelobe.array, sidelobeShift, useInverse=True)
462 dcrModels = DcrModel(modelImages=modelImages, mask=self.mask,
463 effectiveWavelength=self.effectiveWavelength, bandwidth=self.bandwidth)
464 refModels = [dcrModels[subfilter].clone() for subfilter in range(self.dcrNumSubfilters)]
466 dcrModels.regularizeModelFreq(modelImages, self.bbox, statsCtrl, clampFrequency,
467 regularizationWidth=regularizationWidth)
468 for model, refModel, sign in zip(modelImages, refModels, signList):
469 # Make sure the test parameters do reduce the outliers
470 self.assertGreater(np.sum(np.abs(refModel.array - templateImage)),
471 np.sum(np.abs(model.array - templateImage)))
473 def testRegularizeModelIter(self):
474 """Test that large amplitude changes between iterations are restricted.
476 This also tests that noise-like pixels are not regularized.
477 """
478 modelClampFactor = 2.
479 regularizationWidth = 2
480 subfilter = 0
481 dcrModels = DcrModel(modelImages=self.makeTestImages(),
482 effectiveWavelength=self.effectiveWavelength, bandwidth=self.bandwidth)
483 oldModel = dcrModels[0]
484 xSize, ySize = self.bbox.getDimensions()
485 newModel = oldModel.clone()
486 newModel.array[:] += self.rng.rand(ySize, xSize)*np.max(oldModel.array)
487 newModelRef = newModel.clone()
489 dcrModels.regularizeModelIter(subfilter, newModel, self.bbox, modelClampFactor, regularizationWidth)
491 # Make sure the test parameters do reduce the outliers
492 self.assertGreater(np.max(newModelRef.array),
493 np.max(newModel.array - oldModel.array))
494 # Check that all of the outliers are clipped
495 highThreshold = oldModel.array*modelClampFactor
496 highPix = newModel.array > highThreshold
497 highPix = ndimage.binary_opening(highPix, iterations=regularizationWidth)
498 self.assertFalse(np.all(highPix))
499 lowThreshold = oldModel.array/modelClampFactor
500 lowPix = newModel.array < lowThreshold
501 lowPix = ndimage.binary_opening(lowPix, iterations=regularizationWidth)
502 self.assertFalse(np.all(lowPix))
504 def testIterateModel(self):
505 """Test that the DcrModel is iterable, and has the right values.
506 """
507 testModels = self.makeTestImages()
508 refVals = [np.sum(model.array) for model in testModels]
509 dcrModels = DcrModel(modelImages=testModels,
510 effectiveWavelength=self.effectiveWavelength, bandwidth=self.bandwidth)
511 for refVal, model in zip(refVals, dcrModels):
512 self.assertFloatsEqual(refVal, np.sum(model.array))
513 # Negative indices are allowed, so check that those return models from the end.
514 self.assertFloatsEqual(refVals[-1], np.sum(dcrModels[-1].array))
517def calculateAstropyDcr(visitInfo, wcs, effectiveWavelength, bandwidth, dcrNumSubfilters):
518 """Calculate the DCR shift using astropy coordinate transformations.
520 Parameters
521 ----------
522 visitInfo : `lsst.afw.image.VisitInfo`
523 VisitInfo for the exposure.
524 wcs : `lsst.afw.geom.skyWcs.SkyWcs`
525 A wcs that matches the inputs.
526 dcrNumSubfilters : `int`
527 Number of sub-filters used to model chromatic effects within a band.
529 Returns
530 -------
531 dcrShift : `tuple` of two `float`
532 The 2D shift due to DCR, in pixels.
533 Uses numpy axes ordering (Y, X).
534 """
535 elevation = visitInfo.getBoresightAzAlt().getLatitude()
536 azimuth = visitInfo.getBoresightAzAlt().getLongitude()
537 loc = EarthLocation(lat=visitInfo.getObservatory().getLatitude().asDegrees()*u.degree,
538 lon=visitInfo.getObservatory().getLongitude().asDegrees()*u.degree,
539 height=visitInfo.getObservatory().getElevation()*u.m)
540 date = visitInfo.getDate()
541 time = Time(date.get(date.MJD, date.TAI), format='mjd', location=loc, scale='tai')
542 altaz = SkyCoord(alt=elevation.asDegrees(), az=azimuth.asDegrees(),
543 unit='deg', obstime=time, frame='altaz', location=loc)
544 # The DCR calculations are performed at the boresight
545 ra0 = altaz.icrs.ra.degree*degrees
546 dec0 = altaz.icrs.dec.degree*degrees
547 x0, y0 = wcs.skyToPixel(geom.SpherePoint(ra0, dec0))
548 dcrShift = []
549 # We divide the filter into "subfilters" with the full wavelength range
550 # divided into equal sub-ranges.
551 for wl0, wl1 in wavelengthGenerator(effectiveWavelength, bandwidth, dcrNumSubfilters):
552 # Note that diffRefractAmp can be negative,
553 # since it is relative to the midpoint of the full band
554 diffRefractAmp0 = differentialRefraction(wavelength=wl0, wavelengthRef=effectiveWavelength,
555 elevation=elevation,
556 observatory=visitInfo.getObservatory(),
557 weather=visitInfo.getWeather())
558 diffRefractAmp1 = differentialRefraction(wavelength=wl1, wavelengthRef=effectiveWavelength,
559 elevation=elevation,
560 observatory=visitInfo.getObservatory(),
561 weather=visitInfo.getWeather())
562 diffRefractAmp = (diffRefractAmp0 + diffRefractAmp1)/2.
564 elevation1 = elevation + diffRefractAmp
565 altaz = SkyCoord(alt=elevation1.asDegrees(), az=azimuth.asDegrees(),
566 unit='deg', obstime=time, frame='altaz', location=loc)
567 ra1 = altaz.icrs.ra.degree*degrees
568 dec1 = altaz.icrs.dec.degree*degrees
569 x1, y1 = wcs.skyToPixel(geom.SpherePoint(ra1, dec1))
570 dcrShift.append((y1-y0, x1-x0))
571 return dcrShift
574class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
575 pass
578def setup_module(module):
579 lsst.utils.tests.init()
582if __name__ == "__main__": 582 ↛ 583line 582 didn't jump to line 583 because the condition on line 582 was never true
583 lsst.utils.tests.init()
584 unittest.main()