LSST Applications g0f08755f38+82efc23009,g12f32b3c4e+e7bdf1200e,g1653933729+a8ce1bb630,g1a0ca8cf93+50eff2b06f,g28da252d5a+52db39f6a5,g2bbee38e9b+37c5a29d61,g2bc492864f+37c5a29d61,g2cdde0e794+c05ff076ad,g3156d2b45e+41e33cbcdc,g347aa1857d+37c5a29d61,g35bb328faa+a8ce1bb630,g3a166c0a6a+37c5a29d61,g3e281a1b8c+fb992f5633,g414038480c+7f03dfc1b0,g41af890bb2+11b950c980,g5fbc88fb19+17cd334064,g6b1c1869cb+12dd639c9a,g781aacb6e4+a8ce1bb630,g80478fca09+72e9651da0,g82479be7b0+04c31367b4,g858d7b2824+82efc23009,g9125e01d80+a8ce1bb630,g9726552aa6+8047e3811d,ga5288a1d22+e532dc0a0b,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gc28159a63d+37c5a29d61,gcf0d15dbbd+2acd6d4d48,gd7358e8bfb+778a810b6e,gda3e153d99+82efc23009,gda6a2b7d83+2acd6d4d48,gdaeeff99f8+1711a396fd,ge2409df99d+6b12de1076,ge79ae78c31+37c5a29d61,gf0baf85859+d0a5978c5a,gf3967379c6+4954f8c433,gfb92a5be7c+82efc23009,gfec2e1e490+2aaed99252,w.2024.46
LSST Data Management Base Package
Loading...
Searching...
No Matches
isrMock.py
Go to the documentation of this file.
1# This file is part of ip_isr.
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/>.
21
22__all__ = ["IsrMockConfig", "IsrMock", "RawMock", "TrimmedRawMock", "RawDictMock",
23 "CalibratedRawMock", "MasterMock",
24 "BiasMock", "DarkMock", "FlatMock", "FringeMock", "UntrimmedFringeMock",
25 "BfKernelMock", "DefectMock", "CrosstalkCoeffMock", "TransmissionMock",
26 "MockDataContainer", "MockFringeContainer"]
27
28import copy
29import numpy as np
30import tempfile
31
32import lsst.geom
33import lsst.afw.geom as afwGeom
34import lsst.afw.image as afwImage
35from lsstDebug import getDebugFrame
36
37import lsst.afw.cameraGeom.utils as afwUtils
38import lsst.afw.cameraGeom.testUtils as afwTestUtils
39from lsst.afw.cameraGeom import ReadoutCorner
40import lsst.pex.config as pexConfig
41import lsst.pipe.base as pipeBase
42from .crosstalk import CrosstalkCalib
43from .defects import Defects
44
45
46class IsrMockConfig(pexConfig.Config):
47 """Configuration parameters for isrMock.
48
49 These parameters produce generic fixed position signals from
50 various sources, and combine them in a way that matches how those
51 signals are combined to create real data. The camera used is the
52 test camera defined by the afwUtils code.
53 """
54 # Detector parameters. "Exposure" parameters.
55 isLsstLike = pexConfig.Field(
56 dtype=bool,
57 default=False,
58 doc="If True, products have one raw image per amplifier, otherwise, one raw image per detector.",
59 )
60 plateScale = pexConfig.Field(
61 dtype=float,
62 default=20.0,
63 doc="Plate scale used in constructing mock camera.",
64 )
65 radialDistortion = pexConfig.Field(
66 dtype=float,
67 default=0.925,
68 doc="Radial distortion term used in constructing mock camera.",
69 )
70 isTrimmed = pexConfig.Field(
71 dtype=bool,
72 default=True,
73 doc="If True, amplifiers have been trimmed and mosaicked to remove regions outside the data BBox.",
74 )
75 detectorIndex = pexConfig.Field(
76 dtype=int,
77 default=20,
78 doc="Index for the detector to use. The default value uses a standard 2x4 array of amps.",
79 )
80 rngSeed = pexConfig.Field(
81 dtype=int,
82 default=20000913,
83 doc="Seed for random number generator used to add noise.",
84 )
85 # TODO: DM-18345 Check that mocks scale correctly when gain != 1.0
86 gain = pexConfig.Field(
87 dtype=float,
88 default=1.0,
89 doc="Gain for simulated data in electron/adu.",
90 )
91 readNoise = pexConfig.Field(
92 dtype=float,
93 default=5.0,
94 doc="Read noise of the detector in electron.",
95 )
96 expTime = pexConfig.Field(
97 dtype=float,
98 default=5.0,
99 doc="Exposure time for simulated data.",
100 )
101
102 # Signal parameters
103 skyLevel = pexConfig.Field(
104 dtype=float,
105 default=1000.0,
106 doc="Background contribution to be generated from 'the sky' in "
107 "adu (IsrTask) or electron (IsrTaskLSST).",
108 )
109 sourceFlux = pexConfig.ListField(
110 dtype=float,
111 default=[45000.0],
112 doc="Peak flux level of simulated 'astronomical sources' in "
113 "adu (IsrTask) or electron (IsrTaskLSST).",
114 )
115 sourceAmp = pexConfig.ListField(
116 dtype=int,
117 default=[0],
118 doc="Amplifier to place simulated 'astronomical sources'.",
119 )
120 sourceX = pexConfig.ListField(
121 dtype=float,
122 default=[50.0],
123 doc="Peak position (in amplifier coordinates) of simulated 'astronomical sources'.",
124 )
125 sourceY = pexConfig.ListField(
126 dtype=float,
127 default=[25.0],
128 doc="Peak position (in amplifier coordinates) of simulated 'astronomical sources'.",
129 )
130 overscanScale = pexConfig.Field(
131 dtype=float,
132 default=100.0,
133 doc="Amplitude of the ramp function to add to overscan data in "
134 "adu (IsrTask) or electron (IsrTaskLSST)",
135 )
136 biasLevel = pexConfig.Field(
137 dtype=float,
138 default=8000.0,
139 doc="Background contribution to be generated from the bias offset in adu.",
140 )
141 darkRate = pexConfig.Field(
142 dtype=float,
143 default=5.0,
144 doc="Background level contribution (in electron/s) to be generated from dark current.",
145 )
146 darkTime = pexConfig.Field(
147 dtype=float,
148 default=5.0,
149 doc="Exposure time for the dark current contribution.",
150 )
151 flatDrop = pexConfig.Field(
152 dtype=float,
153 default=0.1,
154 doc="Fractional flux drop due to flat from center to edge of detector along x-axis.",
155 )
156 fringeScale = pexConfig.ListField(
157 dtype=float,
158 default=[200.0],
159 doc="Peak fluxes for the components of the fringe ripple in "
160 "adu (IsrTask) or electron (IsrTaskLSST).",
161 )
162 fringeX0 = pexConfig.ListField(
163 dtype=float,
164 default=[-100],
165 doc="Center position for the fringe ripples.",
166 )
167 fringeY0 = pexConfig.ListField(
168 dtype=float,
169 default=[-0],
170 doc="Center position for the fringe ripples.",
171 )
172
173 # Inclusion parameters
174 doAddSky = pexConfig.Field(
175 dtype=bool,
176 default=True,
177 doc="Apply 'sky' signal to output image.",
178 )
179 doAddSource = pexConfig.Field(
180 dtype=bool,
181 default=True,
182 doc="Add simulated source to output image.",
183 )
184 doAddCrosstalk = pexConfig.Field(
185 dtype=bool,
186 default=False,
187 doc="Apply simulated crosstalk to output image. This cannot be corrected by ISR, "
188 "as detector.hasCrosstalk()==False.",
189 )
190 doAddOverscan = pexConfig.Field(
191 dtype=bool,
192 default=True,
193 doc="If untrimmed, add overscan ramp to overscan and data regions.",
194 )
195 doAddBias = pexConfig.Field(
196 dtype=bool,
197 default=True,
198 doc="Add bias signal to data.",
199 )
200 doAddDark = pexConfig.Field(
201 dtype=bool,
202 default=True,
203 doc="Add dark signal to data.",
204 )
205 doAddFlat = pexConfig.Field(
206 dtype=bool,
207 default=True,
208 doc="Add flat signal to data.",
209 )
210 doAddFringe = pexConfig.Field(
211 dtype=bool,
212 default=True,
213 doc="Add fringe signal to data.",
214 )
215
216 # Datasets to create and return instead of generating an image.
217 doTransmissionCurve = pexConfig.Field(
218 dtype=bool,
219 default=False,
220 doc="Return a simulated transmission curve.",
221 )
222 doDefects = pexConfig.Field(
223 dtype=bool,
224 default=False,
225 doc="Return a simulated defect list.",
226 )
227 doBrighterFatter = pexConfig.Field(
228 dtype=bool,
229 default=False,
230 doc="Return a simulated brighter-fatter kernel.",
231 )
232 doDeferredCharge = pexConfig.Field(
233 dtype=bool,
234 default=False,
235 doc="Return a simulated deferred charge calibration.",
236 )
237 doCrosstalkCoeffs = pexConfig.Field(
238 dtype=bool,
239 default=False,
240 doc="Return the matrix of crosstalk coefficients.",
241 )
242 doLinearizer = pexConfig.Field(
243 dtype=bool,
244 default=False,
245 doc="Return linearizer dataset.",
246 )
247 doDataRef = pexConfig.Field(
248 dtype=bool,
249 default=False,
250 doc="Return a simulated gen2 butler dataRef.",
251 )
252 doGenerateImage = pexConfig.Field(
253 dtype=bool,
254 default=False,
255 doc="Return the generated output image if True.",
256 )
257 doGenerateData = pexConfig.Field(
258 dtype=bool,
259 default=False,
260 doc="Return a non-image data structure if True.",
261 )
262 doGenerateAmpDict = pexConfig.Field(
263 dtype=bool,
264 default=False,
265 doc="Return a dict of exposure amplifiers instead of an afwImage.Exposure.",
266 )
267
268
269class IsrMock(pipeBase.Task):
270 """Class to generate consistent mock images for ISR testing.
271
272 ISR testing currently relies on one-off fake images that do not
273 accurately mimic the full set of detector effects. This class
274 uses the test camera/detector/amplifier structure defined in
275 `lsst.afw.cameraGeom.testUtils` to avoid making the test data
276 dependent on any of the actual obs package formats.
277 """
278 ConfigClass = IsrMockConfig
279 _DefaultName = "isrMock"
280
281 def __init__(self, **kwargs):
282 super().__init__(**kwargs)
283 self.rng = np.random.RandomState(self.config.rngSeed)
284 # The coefficients have units adu for IsrTask and have
285 # units electron for IsrTaskLSST.
286 self.crosstalkCoeffs = np.array([[0.0, 0.0, 0.0, 0.0, 0.0, -1e-3, 0.0, 0.0],
287 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
288 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
289 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
290 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
291 [1e-2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
292 [1e-2, 0.0, 0.0, 2.2e-2, 0.0, 0.0, 0.0, 0.0],
293 [1e-2, 5e-3, 5e-4, 3e-3, 4e-2, 5e-3, 5e-3, 0.0]])
294 if getDebugFrame(self._display, "mockCrosstalkCoeffs"):
295 self.crosstalkCoeffs = np.array([[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
296 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
297 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
298 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
299 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
300 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
301 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
302 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
303 # Generic gaussian BF kernel
304 self.bfKernel = np.array([[1., 4., 7., 4., 1.],
305 [4., 16., 26., 16., 4.],
306 [7., 26., 41., 26., 7.],
307 [4., 16., 26., 16., 4.],
308 [1., 4., 7., 4., 1.]]) / 273.0
309
310 def run(self):
311 """Generate a mock ISR product, and return it.
312
313 Returns
314 -------
315 image : `lsst.afw.image.Exposure`
316 Simulated ISR image with signals added.
317 dataProduct :
318 Simulated ISR data products.
319 None :
320 Returned if no valid configuration was found.
321
322 Raises
323 ------
324 RuntimeError
325 Raised if both doGenerateImage and doGenerateData are specified.
326 """
327 if self.config.doGenerateImage and self.config.doGenerateData:
328 raise RuntimeError("Only one of doGenerateImage and doGenerateData may be specified.")
329 elif self.config.doGenerateImage:
330 return self.makeImage()
331 elif self.config.doGenerateData:
332 return self.makeData()
333 else:
334 return None
335
336 def makeData(self):
337 """Generate simulated ISR data.
338
339 Currently, only the class defined crosstalk coefficient
340 matrix, brighter-fatter kernel, a constant unity transmission
341 curve, or a simple single-entry defect list can be generated.
342
343 Returns
344 -------
345 dataProduct :
346 Simulated ISR data product.
347 """
348 if sum(map(bool, [self.config.doBrighterFatter,
349 self.config.doDeferredCharge,
350 self.config.doDefects,
351 self.config.doTransmissionCurve,
352 self.config.doCrosstalkCoeffs,
353 self.config.doLinearizer])) != 1:
354 raise RuntimeError("Only one data product can be generated at a time.")
355 elif self.config.doBrighterFatter:
356 return self.makeBfKernel()
357 elif self.config.doDeferredCharge:
358 return self.makeDeferredChargeCalib()
359 elif self.config.doDefects:
360 return self.makeDefectList()
361 elif self.config.doTransmissionCurve:
362 return self.makeTransmissionCurve()
363 elif self.config.doCrosstalkCoeffs:
364 return self.crosstalkCoeffs
365 elif self.config.doLinearizer:
366 return self.makeLinearizer()
367 else:
368 return None
369
370 def makeBfKernel(self):
371 """Generate a simple Gaussian brighter-fatter kernel.
372
373 Returns
374 -------
375 kernel : `numpy.ndarray`
376 Simulated brighter-fatter kernel.
377 """
378 return self.bfKernel
379
381 """Generate a CTI calibration.
382 """
383
384 raise NotImplementedError("Mock deferred charge is not implemented for IsrMock.")
385
386 def makeDefectList(self):
387 """Generate a simple single-entry defect list.
388
389 Returns
390 -------
391 defectList : `lsst.meas.algorithms.Defects`
392 Simulated defect list
393 """
395 lsst.geom.Extent2I(40, 50))])
396
398 """Generate the simulated crosstalk coefficients.
399
400 Returns
401 -------
402 coeffs : `numpy.ndarray`
403 Simulated crosstalk coefficients.
404 """
405
406 return self.crosstalkCoeffs
407
409 """Generate a simulated flat transmission curve.
410
411 Returns
412 -------
413 transmission : `lsst.afw.image.TransmissionCurve`
414 Simulated transmission curve.
415 """
416
417 return afwImage.TransmissionCurve.makeIdentity()
418
419 def makeLinearity(self):
420 """Generate a linearity dataset.
421
422 Returns
423 -------
424 linearizer : `lsst.ip.isr.Linearizer`
425 """
426 raise NotImplementedError("Linearizer not implemented for isrMock.")
427
428 def makeImage(self):
429 """Generate a simulated ISR image.
430
431 Returns
432 -------
433 exposure : `lsst.afw.image.Exposure` or `dict`
434 Simulated ISR image data.
435
436 Notes
437 -----
438 This method currently constructs a "raw" data image by:
439
440 * Generating a simulated sky with noise
441 * Adding a single Gaussian "star"
442 * Adding the fringe signal
443 * Multiplying the frame by the simulated flat
444 * Adding dark current (and noise)
445 * Adding a bias offset (and noise)
446 * Adding an overscan gradient parallel to the pixel y-axis
447 * Simulating crosstalk by adding a scaled version of each
448 amplifier to each other amplifier.
449
450 The exposure with image data constructed this way is in one of
451 three formats.
452
453 * A single image, with overscan and prescan regions retained
454 * A single image, with overscan and prescan regions trimmed
455 * A `dict`, containing the amplifer data indexed by the
456 amplifier name.
457
458 The nonlinearity, CTE, and brighter fatter are currently not
459 implemented.
460
461 Note that this method generates an image in the reverse
462 direction as the ISR processing, as the output image here has
463 had a series of instrument effects added to an idealized
464 exposure.
465 """
466 exposure = self.getExposure()
467
468 for idx, amp in enumerate(exposure.getDetector()):
469 bbox = None
470 if self.config.isTrimmed is True:
471 bbox = amp.getBBox()
472 else:
473 bbox = amp.getRawDataBBox()
474
475 ampData = exposure.image[bbox]
476
477 if self.config.doAddSky is True:
478 self.amplifierAddNoise(ampData, self.config.skyLevel, np.sqrt(self.config.skyLevel))
479
480 if self.config.doAddSource is True:
481 for sourceAmp, sourceFlux, sourceX, sourceY in zip(self.config.sourceAmp,
482 self.config.sourceFlux,
483 self.config.sourceX,
484 self.config.sourceY):
485 if idx == sourceAmp:
486 self.amplifierAddSource(ampData, sourceFlux, sourceX, sourceY)
487
488 if self.config.doAddFringe is True:
489 self.amplifierAddFringe(amp, ampData, np.array(self.config.fringeScale),
490 x0=np.array(self.config.fringeX0),
491 y0=np.array(self.config.fringeY0))
492
493 if self.config.doAddFlat is True:
494 if ampData.getArray().sum() == 0.0:
495 self.amplifierAddNoise(ampData, 1.0, 0.0)
496 u0 = exposure.getDimensions().getX()
497 v0 = exposure.getDimensions().getY()
498 self.amplifierMultiplyFlat(amp, ampData, self.config.flatDrop, u0=u0, v0=v0)
499
500 if self.config.doAddDark is True:
501 self.amplifierAddNoise(ampData,
502 self.config.darkRate * self.config.darkTime / self.config.gain,
503 np.sqrt(self.config.darkRate
504 * self.config.darkTime / self.config.gain))
505
506 if self.config.doAddCrosstalk is True:
507 ctCalib = CrosstalkCalib()
508 for idxS, ampS in enumerate(exposure.getDetector()):
509 for idxT, ampT in enumerate(exposure.getDetector()):
510 ampDataT = exposure.image[ampT.getBBox()
511 if self.config.isTrimmed else ampT.getRawDataBBox()]
512 outAmp = ctCalib.extractAmp(exposure.getImage(), ampS, ampT,
513 isTrimmed=self.config.isTrimmed)
514 self.amplifierAddCT(outAmp, ampDataT, self.crosstalkCoeffs[idxS][idxT])
515
516 for amp in exposure.getDetector():
517 bbox = None
518 if self.config.isTrimmed is True:
519 bbox = amp.getBBox()
520 else:
521 bbox = amp.getRawDataBBox()
522
523 ampData = exposure.image[bbox]
524
525 if self.config.doAddBias is True:
526 self.amplifierAddNoise(ampData, self.config.biasLevel,
527 self.config.readNoise / self.config.gain)
528
529 if self.config.doAddOverscan is True:
530 oscanBBox = amp.getRawHorizontalOverscanBBox()
531 oscanData = exposure.image[oscanBBox]
532 self.amplifierAddNoise(oscanData, self.config.biasLevel,
533 self.config.readNoise / self.config.gain)
534
535 self.amplifierAddYGradient(ampData, -1.0 * self.config.overscanScale,
536 1.0 * self.config.overscanScale)
537 self.amplifierAddYGradient(oscanData, -1.0 * self.config.overscanScale,
538 1.0 * self.config.overscanScale)
539
540 if self.config.doGenerateAmpDict is True:
541 expDict = dict()
542 for amp in exposure.getDetector():
543 expDict[amp.getName()] = exposure
544 return expDict
545 else:
546 return exposure
547
548 # afw primatives to construct the image structure
549 def getCamera(self, isForAssembly=False):
550 """Construct a test camera object.
551
552 Parameters
553 -------
554 isForAssembly : `bool`
555 If True, construct a camera with "super raw"
556 orientation (all amplifiers have LL readout
557 corner but still contains the necessary flip
558 and offset info needed for assembly. This is
559 needed if isLsstLike is True. If False, return
560 a camera with bboxes flipped and offset to the
561 correct orientation given the readout corner.
562
563 Returns
564 -------
565 camera : `lsst.afw.cameraGeom.camera`
566 Test camera.
567 """
568 cameraWrapper = afwTestUtils.CameraWrapper(
569 plateScale=self.config.plateScale,
570 radialDistortion=self.config.radialDistortion,
571 isLsstLike=self.config.isLsstLike and isForAssembly,
572 )
573 camera = cameraWrapper.camera
574 return camera
575
576 def getExposure(self, isTrimmed=None):
577 """Construct a test exposure.
578
579 The test exposure has a simple WCS set, as well as a list of
580 unlikely header keywords that can be removed during ISR
581 processing to exercise that code.
582
583 Parameters
584 ----------
585 isTrimmed : `bool` or `None`, optional
586 Override the configuration isTrimmed?
587
588 Returns
589 -------
590 exposure : `lsst.afw.exposure.Exposure`
591 Construct exposure containing masked image of the
592 appropriate size.
593 """
594 if isTrimmed is None:
595 _isTrimmed = self.config.isTrimmed
596 else:
597 _isTrimmed = isTrimmed
598
599 camera = self.getCamera(isForAssembly=self.config.isLsstLike)
600 detector = camera[self.config.detectorIndex]
601 image = afwUtils.makeImageFromCcd(
602 detector,
603 isTrimmed=_isTrimmed,
604 showAmpGain=False,
605 rcMarkSize=0,
606 binSize=1,
607 imageFactory=afwImage.ImageF,
608 )
609
610 var = afwImage.ImageF(image.getDimensions())
611 mask = afwImage.Mask(image.getDimensions())
612 image.assign(0.0)
613
614 maskedImage = afwImage.makeMaskedImage(image, mask, var)
615 exposure = afwImage.makeExposure(maskedImage)
616 exposure.setDetector(detector)
617 exposure.setWcs(self.getWcs())
618
619 visitInfo = afwImage.VisitInfo(exposureTime=self.config.expTime, darkTime=self.config.darkTime)
620 exposure.getInfo().setVisitInfo(visitInfo)
621 # Set a dummy ID.
622 exposure.getInfo().setId(12345)
623
624 metadata = exposure.getMetadata()
625 metadata.add("SHEEP", 7.3, "number of sheep on farm")
626 metadata.add("MONKEYS", 155, "monkeys per tree")
627 metadata.add("VAMPIRES", 4, "How scary are vampires.")
628
629 ccd = exposure.getDetector()
630 newCcd = ccd.rebuild()
631 newCcd.clear()
632 readoutMap = {
633 'LL': ReadoutCorner.LL,
634 'LR': ReadoutCorner.LR,
635 'UR': ReadoutCorner.UR,
636 'UL': ReadoutCorner.UL,
637 }
638 for amp in ccd:
639 newAmp = amp.rebuild()
640 newAmp.setLinearityCoeffs((0., 1., 0., 0.))
641 newAmp.setLinearityType("Polynomial")
642 newAmp.setGain(self.config.gain)
643 newAmp.setSuspectLevel(25000.0)
644 newAmp.setSaturation(32000.0)
645 readoutCorner = amp.getReadoutCorner().name
646
647 # Apply flips to bbox where needed
648 imageBBox = amp.getRawDataBBox()
649 rawBbox = amp.getRawBBox()
650 parallelOscanBBox = amp.getRawParallelOverscanBBox()
651 serialOscanBBox = amp.getRawSerialOverscanBBox()
652 prescanBBox = amp.getRawPrescanBBox()
653
654 if self.config.isLsstLike:
655 # This follows cameraGeom.testUtils
656 xoffset, yoffset = amp.getRawXYOffset()
657 offext = lsst.geom.Extent2I(xoffset, yoffset)
658 flipx = bool(amp.getRawFlipX())
659 flipy = bool(amp.getRawFlipY())
660 if flipx:
661 xExt = rawBbox.getDimensions().getX()
662 rawBbox.flipLR(xExt)
663 imageBBox.flipLR(xExt)
664 parallelOscanBBox.flipLR(xExt)
665 serialOscanBBox.flipLR(xExt)
666 prescanBBox.flipLR(xExt)
667 if flipy:
668 yExt = rawBbox.getDimensions().getY()
669 rawBbox.flipTB(yExt)
670 imageBBox.flipTB(yExt)
671 parallelOscanBBox.flipTB(yExt)
672 serialOscanBBox.flipTB(yExt)
673 prescanBBox.flipTB(yExt)
674 if not flipx and not flipy:
675 readoutCorner = 'LL'
676 elif flipx and not flipy:
677 readoutCorner = 'LR'
678 elif flipx and flipy:
679 readoutCorner = 'UR'
680 elif not flipx and flipy:
681 readoutCorner = 'UL'
682 rawBbox.shift(offext)
683 imageBBox.shift(offext)
684 parallelOscanBBox.shift(offext)
685 serialOscanBBox.shift(offext)
686 prescanBBox.shift(offext)
687 newAmp.setReadoutCorner(readoutMap[readoutCorner])
688 newAmp.setRawBBox(rawBbox)
689 newAmp.setRawDataBBox(imageBBox)
690 newAmp.setRawParallelOverscanBBox(parallelOscanBBox)
691 newAmp.setRawSerialOverscanBBox(serialOscanBBox)
692 newAmp.setRawPrescanBBox(prescanBBox)
693 newAmp.setRawFlipX(False)
694 newAmp.setRawFlipY(False)
695 no_offset = lsst.geom.Extent2I(0, 0)
696 newAmp.setRawXYOffset(no_offset)
697
698 newCcd.append(newAmp)
699
700 exposure.setDetector(newCcd.finish())
701
702 exposure.image.array[:] = np.zeros(exposure.getImage().getDimensions()).transpose()
703 exposure.mask.array[:] = np.zeros(exposure.getMask().getDimensions()).transpose()
704 exposure.variance.array[:] = np.zeros(exposure.getVariance().getDimensions()).transpose()
705
706 return exposure
707
708 def getWcs(self):
709 """Construct a dummy WCS object.
710
711 Taken from the deprecated ip_isr/examples/exampleUtils.py.
712
713 This is not guaranteed, given the distortion and pixel scale
714 listed in the afwTestUtils camera definition.
715
716 Returns
717 -------
718 wcs : `lsst.afw.geom.SkyWcs`
719 Test WCS transform.
720 """
721 return afwGeom.makeSkyWcs(crpix=lsst.geom.Point2D(0.0, 100.0),
722 crval=lsst.geom.SpherePoint(45.0, 25.0, lsst.geom.degrees),
723 cdMatrix=afwGeom.makeCdMatrix(scale=1.0*lsst.geom.degrees))
724
725 def localCoordToExpCoord(self, ampData, x, y):
726 """Convert between a local amplifier coordinate and the full
727 exposure coordinate.
728
729 Parameters
730 ----------
731 ampData : `lsst.afw.image.ImageF`
732 Amplifier image to use for conversions.
733 x : `int`
734 X-coordinate of the point to transform.
735 y : `int`
736 Y-coordinate of the point to transform.
737
738 Returns
739 -------
740 u : `int`
741 Transformed x-coordinate.
742 v : `int`
743 Transformed y-coordinate.
744
745 Notes
746 -----
747 The output is transposed intentionally here, to match the
748 internal transpose between numpy and afw.image coordinates.
749 """
750 u = x + ampData.getBBox().getBeginX()
751 v = y + ampData.getBBox().getBeginY()
752
753 return (v, u)
754
755 # Simple data values.
756 def amplifierAddNoise(self, ampData, mean, sigma, rng=None):
757 """Add Gaussian noise to an amplifier's image data.
758
759 This method operates in the amplifier coordinate frame.
760
761 Parameters
762 ----------
763 ampData : `lsst.afw.image.ImageF`
764 Amplifier image to operate on.
765 mean : `float`
766 Mean value of the Gaussian noise.
767 sigma : `float`
768 Sigma of the Gaussian noise.
769 rng : `np.random.RandomState`, optional
770 Random state to use instead of self.rng.
771 """
772 if rng is not None:
773 _rng = rng
774 else:
775 _rng = self.rng
776
777 ampArr = ampData.array
778 ampArr[:] = ampArr[:] + _rng.normal(mean, sigma,
779 size=ampData.getDimensions()).transpose()
780
781 def amplifierAddYGradient(self, ampData, start, end):
782 """Add a y-axis linear gradient to an amplifier's image data.
783
784 This method operates in the amplifier coordinate frame.
785
786 Parameters
787 ----------
788 ampData : `lsst.afw.image.ImageF`
789 Amplifier image to operate on.
790 start : `float`
791 Start value of the gradient (at y=0).
792 end : `float`
793 End value of the gradient (at y=ymax).
794 """
795 nPixY = ampData.getDimensions().getY()
796 ampArr = ampData.array
797 ampArr[:] = ampArr[:] + (np.interp(range(nPixY), (0, nPixY - 1), (start, end)).reshape(nPixY, 1)
798 + np.zeros(ampData.getDimensions()).transpose())
799
800 def amplifierAddSource(self, ampData, scale, x0, y0):
801 """Add a single Gaussian source to an amplifier.
802
803 This method operates in the amplifier coordinate frame.
804
805 Parameters
806 ----------
807 ampData : `lsst.afw.image.ImageF`
808 Amplifier image to operate on.
809 scale : `float`
810 Peak flux of the source to add.
811 x0 : `float`
812 X-coordinate of the source peak.
813 y0 : `float`
814 Y-coordinate of the source peak.
815 """
816 for x in range(0, ampData.getDimensions().getX()):
817 for y in range(0, ampData.getDimensions().getY()):
818 ampData.array[y][x] = (ampData.array[y][x]
819 + scale * np.exp(-0.5 * ((x - x0)**2 + (y - y0)**2) / 3.0**2))
820
821 def amplifierAddCT(self, ampDataSource, ampDataTarget, scale):
822 """Add a scaled copy of an amplifier to another, simulating crosstalk.
823
824 This method operates in the amplifier coordinate frame.
825
826 Parameters
827 ----------
828 ampDataSource : `lsst.afw.image.ImageF`
829 Amplifier image to add scaled copy from.
830 ampDataTarget : `lsst.afw.image.ImageF`
831 Amplifier image to add scaled copy to.
832 scale : `float`
833 Flux scale of the copy to add to the target.
834
835 Notes
836 -----
837 This simulates simple crosstalk between amplifiers.
838 """
839 ampDataTarget.array[:] = (ampDataTarget.array[:]
840 + scale * ampDataSource.array[:])
841
842 # Functional form data values.
843 def amplifierAddFringe(self, amp, ampData, scale, x0=100, y0=0):
844 """Add a fringe-like ripple pattern to an amplifier's image data.
845
846 Parameters
847 ----------
848 amp : `~lsst.afw.ampInfo.AmpInfoRecord`
849 Amplifier to operate on. Needed for amp<->exp coordinate
850 transforms.
851 ampData : `lsst.afw.image.ImageF`
852 Amplifier image to operate on.
853 scale : `numpy.array` or `float`
854 Peak intensity scaling for the ripple.
855 x0 : `numpy.array` or `float`, optional
856 Fringe center
857 y0 : `numpy.array` or `float`, optional
858 Fringe center
859
860 Notes
861 -----
862 This uses an offset sinc function to generate a ripple
863 pattern. True fringes have much finer structure, but this
864 pattern should be visually identifiable. The (x, y)
865 coordinates are in the frame of the amplifier, and (u, v) in
866 the frame of the full trimmed image.
867 """
868 for x in range(0, ampData.getDimensions().getX()):
869 for y in range(0, ampData.getDimensions().getY()):
870 (u, v) = self.localCoordToExpCoord(amp, x, y)
871 ampData.getArray()[y][x] = np.sum((ampData.getArray()[y][x]
872 + scale * np.sinc(((u - x0) / 50)**2
873 + ((v - y0) / 50)**2)))
874
875 def amplifierMultiplyFlat(self, amp, ampData, fracDrop, u0=100.0, v0=100.0):
876 """Multiply an amplifier's image data by a flat-like pattern.
877
878 Parameters
879 ----------
880 amp : `lsst.afw.ampInfo.AmpInfoRecord`
881 Amplifier to operate on. Needed for amp<->exp coordinate
882 transforms.
883 ampData : `lsst.afw.image.ImageF`
884 Amplifier image to operate on.
885 fracDrop : `float`
886 Fractional drop from center to edge of detector along x-axis.
887 u0 : `float`
888 Peak location in detector coordinates.
889 v0 : `float`
890 Peak location in detector coordinates.
891
892 Notes
893 -----
894 This uses a 2-d Gaussian to simulate an illumination pattern
895 that falls off towards the edge of the detector. The (x, y)
896 coordinates are in the frame of the amplifier, and (u, v) in
897 the frame of the full trimmed image.
898 """
899 if fracDrop >= 1.0:
900 raise RuntimeError("Flat fractional drop cannot be greater than 1.0")
901
902 sigma = u0 / np.sqrt(-2.0 * np.log(fracDrop))
903
904 for x in range(0, ampData.getDimensions().getX()):
905 for y in range(0, ampData.getDimensions().getY()):
906 (u, v) = self.localCoordToExpCoord(amp, x, y)
907 f = np.exp(-0.5 * ((u - u0)**2 + (v - v0)**2) / sigma**2)
908 ampData.array[y][x] = (ampData.array[y][x] * f)
909
910
912 """Generate a raw exposure suitable for ISR.
913 """
914 def __init__(self, **kwargs):
915 super().__init__(**kwargs)
916 self.config.isTrimmed = False
917 self.config.doGenerateImage = True
918 self.config.doGenerateAmpDict = False
919 self.config.doAddOverscan = True
920 self.config.doAddSky = True
921 self.config.doAddSource = True
922 self.config.doAddCrosstalk = False
923 self.config.doAddBias = True
924 self.config.doAddDark = True
925
926
928 """Generate a trimmed raw exposure.
929 """
930 def __init__(self, **kwargs):
931 super().__init__(**kwargs)
932 self.config.isTrimmed = True
933 self.config.doAddOverscan = False
934
935
937 """Generate a trimmed raw exposure.
938 """
939 def __init__(self, **kwargs):
940 super().__init__(**kwargs)
941 self.config.isTrimmed = True
942 self.config.doGenerateImage = True
943 self.config.doAddOverscan = False
944 self.config.doAddSky = True
945 self.config.doAddSource = True
946 self.config.doAddCrosstalk = False
947
948 self.config.doAddBias = False
949 self.config.doAddDark = False
950 self.config.doAddFlat = False
951 self.config.doAddFringe = True
952
953 self.config.biasLevel = 0.0
954 self.config.readNoise = 10.0
955
956
958 """Generate a raw exposure dict suitable for ISR.
959 """
960 def __init__(self, **kwargs):
961 super().__init__(**kwargs)
962 self.config.doGenerateAmpDict = True
963
964
966 """Parent class for those that make master calibrations.
967 """
968 def __init__(self, **kwargs):
969 super().__init__(**kwargs)
970 self.config.isTrimmed = True
971 self.config.doGenerateImage = True
972 self.config.doAddOverscan = False
973 self.config.doAddSky = False
974 self.config.doAddSource = False
975 self.config.doAddCrosstalk = False
976
977 self.config.doAddBias = False
978 self.config.doAddDark = False
979 self.config.doAddFlat = False
980 self.config.doAddFringe = False
981
982
984 """Simulated master bias calibration.
985 """
986 def __init__(self, **kwargs):
987 super().__init__(**kwargs)
988 self.config.doAddBias = True
989 self.config.readNoise = 10.0
990
991
993 """Simulated master dark calibration.
994 """
995 def __init__(self, **kwargs):
996 super().__init__(**kwargs)
997 self.config.doAddDark = True
998 self.config.darkTime = 1.0
999
1000
1002 """Simulated master flat calibration.
1003 """
1004 def __init__(self, **kwargs):
1005 super().__init__(**kwargs)
1006 self.config.doAddFlat = True
1007
1008
1010 """Simulated master fringe calibration.
1011 """
1012 def __init__(self, **kwargs):
1013 super().__init__(**kwargs)
1014 self.config.doAddFringe = True
1015
1016
1018 """Simulated untrimmed master fringe calibration.
1019 """
1020 def __init__(self, **kwargs):
1021 super().__init__(**kwargs)
1022 self.config.isTrimmed = False
1023
1024
1026 """Simulated brighter-fatter kernel.
1027 """
1028 def __init__(self, **kwargs):
1029 super().__init__(**kwargs)
1030 self.config.doGenerateImage = False
1031 self.config.doGenerateData = True
1032 self.config.doBrighterFatter = True
1033 self.config.doDefects = False
1034 self.config.doCrosstalkCoeffs = False
1035 self.config.doTransmissionCurve = False
1036
1037
1039 """Simulated deferred charge calibration.
1040 """
1041 def __init__(self, **kwargs):
1042 super().__init__(**kwargs)
1043 self.config.doGenerateImage = False
1044 self.config.doGenerateData = True
1045 self.config.doDeferredCharge = True
1046 self.config.doDefects = False
1047 self.config.doCrosstalkCoeffs = False
1048 self.config.doTransmissionCurve = False
1049
1050
1052 """Simulated defect list.
1053 """
1054 def __init__(self, **kwargs):
1055 super().__init__(**kwargs)
1056 self.config.doGenerateImage = False
1057 self.config.doGenerateData = True
1058 self.config.doBrighterFatter = False
1059 self.config.doDefects = True
1060 self.config.doCrosstalkCoeffs = False
1061 self.config.doTransmissionCurve = False
1062
1063
1065 """Simulated crosstalk coefficient matrix.
1066 """
1067 def __init__(self, **kwargs):
1068 super().__init__(**kwargs)
1069 self.config.doGenerateImage = False
1070 self.config.doGenerateData = True
1071 self.config.doBrighterFatter = False
1072 self.config.doDefects = False
1073 self.config.doCrosstalkCoeffs = True
1074 self.config.doTransmissionCurve = False
1075
1076
1078 """Simulated transmission curve.
1079 """
1080 def __init__(self, **kwargs):
1081 super().__init__(**kwargs)
1082 self.config.doGenerateImage = False
1083 self.config.doGenerateData = True
1084 self.config.doBrighterFatter = False
1085 self.config.doDefects = False
1086 self.config.doCrosstalkCoeffs = False
1087 self.config.doTransmissionCurve = True
1088
1089
1090class MockDataContainer(object):
1091 """Container for holding ISR mock objects.
1092 """
1093 dataId = "isrMock Fake Data"
1094 darkval = 2. # electron/sec
1095 oscan = 250. # adu
1096 gradient = .10
1097 exptime = 15.0 # seconds
1098 darkexptime = 15.0 # seconds
1099
1100 def __init__(self, **kwargs):
1101 if 'config' in kwargs.keys():
1102 self.config = kwargs['config']
1103 else:
1104 self.config = None
1105
1106 def expectImage(self):
1107 if self.config is None:
1108 self.config = IsrMockConfig()
1109 self.config.doGenerateImage = True
1110 self.config.doGenerateData = False
1111
1112 def expectData(self):
1113 if self.config is None:
1114 self.config = IsrMockConfig()
1115 self.config.doGenerateImage = False
1116 self.config.doGenerateData = True
1117
1118 def get(self, dataType, **kwargs):
1119 """Return an appropriate data product.
1120
1121 Parameters
1122 ----------
1123 dataType : `str`
1124 Type of data product to return.
1125
1126 Returns
1127 -------
1128 mock : IsrMock.run() result
1129 The output product.
1130 """
1131 if "_filename" in dataType:
1132 self.expectData()
1133 return tempfile.mktemp(), "mock"
1134 elif 'transmission_' in dataType:
1135 self.expectData()
1136 return TransmissionMock(config=self.config).run()
1137 elif dataType == 'ccdExposureId':
1138 self.expectData()
1139 return 20090913
1140 elif dataType == 'camera':
1141 self.expectData()
1142 return IsrMock(config=self.config).getCamera()
1143 elif dataType == 'raw':
1144 self.expectImage()
1145 return RawMock(config=self.config).run()
1146 elif dataType == 'bias':
1147 self.expectImage()
1148 return BiasMock(config=self.config).run()
1149 elif dataType == 'dark':
1150 self.expectImage()
1151 return DarkMock(config=self.config).run()
1152 elif dataType == 'flat':
1153 self.expectImage()
1154 return FlatMock(config=self.config).run()
1155 elif dataType == 'fringe':
1156 self.expectImage()
1157 return FringeMock(config=self.config).run()
1158 elif dataType == 'defects':
1159 self.expectData()
1160 return DefectMock(config=self.config).run()
1161 elif dataType == 'bfKernel':
1162 self.expectData()
1163 return BfKernelMock(config=self.config).run()
1164 elif dataType == 'linearizer':
1165 return None
1166 elif dataType == 'crosstalkSources':
1167 return None
1168 else:
1169 raise RuntimeError("ISR DataRefMock cannot return %s.", dataType)
1170
1171
1173 """Container for mock fringe data.
1174 """
1175 dataId = "isrMock Fake Data"
1176 darkval = 2. # electron/sec
1177 oscan = 250. # adu
1178 gradient = .10
1179 exptime = 15 # seconds
1180 darkexptime = 40. # seconds
1181
1182 def __init__(self, **kwargs):
1183 if 'config' in kwargs.keys():
1184 self.config = kwargs['config']
1185 else:
1186 self.config = IsrMockConfig()
1187 self.config.isTrimmed = True
1188 self.config.doAddFringe = True
1189 self.config.readNoise = 10.0
1190
1191 def get(self, dataType, **kwargs):
1192 """Return an appropriate data product.
1193
1194 Parameters
1195 ----------
1196 dataType : `str`
1197 Type of data product to return.
1198
1199 Returns
1200 -------
1201 mock : IsrMock.run() result
1202 The output product.
1203 """
1204 if "_filename" in dataType:
1205 return tempfile.mktemp(), "mock"
1206 elif 'transmission_' in dataType:
1207 return TransmissionMock(config=self.config).run()
1208 elif dataType == 'ccdExposureId':
1209 return 20090913
1210 elif dataType == 'camera':
1211 return IsrMock(config=self.config).getCamera()
1212 elif dataType == 'raw':
1213 return CalibratedRawMock(config=self.config).run()
1214 elif dataType == 'bias':
1215 return BiasMock(config=self.config).run()
1216 elif dataType == 'dark':
1217 return DarkMock(config=self.config).run()
1218 elif dataType == 'flat':
1219 return FlatMock(config=self.config).run()
1220 elif dataType == 'fringe':
1221 fringes = []
1222 configCopy = copy.deepcopy(self.config)
1223 for scale, x, y in zip(self.config.fringeScale, self.config.fringeX0, self.config.fringeY0):
1224 configCopy.fringeScale = [1.0]
1225 configCopy.fringeX0 = [x]
1226 configCopy.fringeY0 = [y]
1227 fringes.append(FringeMock(config=configCopy).run())
1228 return fringes
1229 elif dataType == 'defects':
1230 return DefectMock(config=self.config).run()
1231 elif dataType == 'bfKernel':
1232 return BfKernelMock(config=self.config).run()
1233 elif dataType == 'linearizer':
1234 return None
1235 elif dataType == 'crosstalkSources':
1236 return None
1237 else:
1238 return None
Represent a 2-dimensional array of bitmask pixels.
Definition Mask.h:81
Information about a single exposure of an imaging camera.
Definition VisitInfo.h:68
An integer coordinate rectangle.
Definition Box.h:55
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
localCoordToExpCoord(self, ampData, x, y)
Definition isrMock.py:725
amplifierAddNoise(self, ampData, mean, sigma, rng=None)
Definition isrMock.py:756
amplifierAddSource(self, ampData, scale, x0, y0)
Definition isrMock.py:800
amplifierMultiplyFlat(self, amp, ampData, fracDrop, u0=100.0, v0=100.0)
Definition isrMock.py:875
amplifierAddYGradient(self, ampData, start, end)
Definition isrMock.py:781
getCamera(self, isForAssembly=False)
Definition isrMock.py:549
getExposure(self, isTrimmed=None)
Definition isrMock.py:576
amplifierAddCT(self, ampDataSource, ampDataTarget, scale)
Definition isrMock.py:821
amplifierAddFringe(self, amp, ampData, scale, x0=100, y0=0)
Definition isrMock.py:843
__init__(self, **kwargs)
Definition isrMock.py:281
get(self, dataType, **kwargs)
Definition isrMock.py:1118
get(self, dataType, **kwargs)
Definition isrMock.py:1191
__init__(self, **kwargs)
Definition isrMock.py:914
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition SkyWcs.cc:521
Eigen::Matrix2d makeCdMatrix(lsst::geom::Angle const &scale, lsst::geom::Angle const &orientation=0 *lsst::geom::degrees, bool flipX=false)
Make a WCS CD matrix.
Definition SkyWcs.cc:133
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename std::shared_ptr< Image< ImagePixelT > > image, typename std::shared_ptr< Mask< MaskPixelT > > mask=Mask< MaskPixelT >(), typename std::shared_ptr< Image< VariancePixelT > > variance=Image< VariancePixelT >())
A function to return a MaskedImage of the correct type (cf.
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
Definition Exposure.h:484