LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
mockCoadd.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010, 2011, 2012 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
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 <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 
23 import lsst.afw.image
24 import lsst.afw.geom
25 import lsst.pex.config
26 import lsst.afw.table
27 import lsst.pipe.base
28 from lsst.pipe.tasks.makeSkyMap import MakeSkyMapTask
29 from lsst.pipe.tasks.makeCoaddTempExp import MakeCoaddTempExpTask
30 from lsst.pipe.tasks.assembleCoadd import AssembleCoaddTask
31 from .mockObject import MockObjectTask
32 from .mockObservation import MockObservationTask
33 from .mockSelect import MockSelectImagesTask
34 
35 class MockCoaddConfig(lsst.pex.config.Config):
36  makeSkyMap = lsst.pex.config.ConfigurableField(
37  doc = "SkyMap builder subtask",
38  target = MakeSkyMapTask
39  )
40  mockObject = lsst.pex.config.ConfigurableField(
41  doc = "Subtask that generates and draws the objects/sources in the mock images",
42  target = MockObjectTask
43  )
44  mockObservation = lsst.pex.config.ConfigurableField(
45  doc = "Subtask that generates the Wcs, Psf, Calib, etc. of mock images",
46  target = MockObservationTask
47  )
48  coaddName = lsst.pex.config.Field(
49  doc = "Coadd name used as a prefix for other datasets",
50  dtype = str,
51  optional = False,
52  default = "deep"
53  )
54  nObservations = lsst.pex.config.Field(
55  doc = "Number of mock observations to generate.",
56  dtype = int,
57  optional = False,
58  default = 12
59  )
60  edgeBuffer = lsst.pex.config.Field(
61  doc = ("Number of pixels by which to grow object bounding boxes when determining whether they land "
62  " completely on a generated image"),
63  dtype = int,
64  optional = False,
65  default = 5
66  )
67 
68  def setupSkyMapPatches(self, nPatches=2, patchSize=400, pixelScale = 0.2*lsst.afw.geom.arcseconds):
69  """
70  Set the nested [discrete] skymap config parameters such that the full tract
71  has nPatches x nPatches patches of the given size and pixel scale.
72  """
73  self.makeSkyMap.skyMap['discrete'].patchInnerDimensions = [patchSize, patchSize]
74  self.makeSkyMap.skyMap['discrete'].pixelScale = pixelScale.asArcseconds()
75  # multiply by 0.5 because we want a half-width; subtract 0.49 to ensure that we get the right
76  # number after skyMap.TractInfo rounds up.
77  radius = (0.5 * nPatches - 0.49) * patchSize * pixelScale.asDegrees()
78  self.makeSkyMap.skyMap['discrete'].radiusList = [radius]
79 
80  def setDefaults(self):
81  self.makeSkyMap.skyMap.name = 'discrete'
82  self.makeSkyMap.skyMap['discrete'].raList = [90.0]
83  self.makeSkyMap.skyMap['discrete'].decList = [0.0]
84  self.makeSkyMap.skyMap['discrete'].patchBorder = 10
85  self.makeSkyMap.skyMap['discrete'].projection = "TAN"
86  self.makeSkyMap.skyMap['discrete'].tractOverlap = 0.0
87  self.setupSkyMapPatches()
88 
89 class MockCoaddTask(lsst.pipe.base.CmdLineTask):
90  """MockCoaddTask is a driver task for creating mock coadds. As opposed to more realistic
91  simulations, MockCoadd generates and uses extremely simple "toy" data that can be used to more
92  rigorously test the behavior of high-level task code because the expected results are
93  more easily predicted. In particular, calexps are generated directly from the truth catalog,
94  and contain only zero-noise stars that are created using the same Psf, Calib, and Wcs that will
95  be attached to the mock calexp.
96 
97  In addition to creating the mock calexps and truth catalogs, MockCoadd also contains driver
98  code to run the MakeSkyMap, MakeCoaddTempExp, and AssembleCoadd tasks on the mock calexps,
99  and code to directly create a mock coadd image using CoaddPsf, which can be compared to the
100  output of the regular coadd tasks to check that the coadd code and CoaddPsf are consistent.
101 
102  Note that aside from MakeSkyMapTask, the coadd tasks are *not* subtasks of MockCoaddTasks,
103  and their configs are not part of MockCoaddConfig; these are created locally within
104  MockCoaddTask methods when needed, as not all coadd task config options are appropriate
105  for the mock data generated by MockCoadd.
106  """
107 
108  ConfigClass = MockCoaddConfig
109 
110  _DefaultName = "MockCoadd"
111 
112  def __init__(self, **kwds):
113  """Construct a MockCoaddTask and the subtasks used for generating skymaps, objects,
114  and observations (i.e. calexp parameters).
115  """
116  lsst.pipe.base.CmdLineTask.__init__(self, **kwds)
117  self.makeSubtask("makeSkyMap")
118  self.makeSubtask("mockObject")
119  self.makeSubtask("mockObservation")
121  self.objectIdKey = self.schema.addField("objectId", type=long, doc="foreign key to truth catalog")
122  self.exposureIdKey = self.schema.addField("exposureId", type=long,
123  doc="foreign key to observation catalog")
124  self.centroidInBBoxKey = self.schema.addField(
125  "centroidInBBox", type="Flag",
126  doc="set if this source's center position is inside the generated image's bbox"
127  )
128  self.partialOverlapKey = self.schema.addField(
129  "partialOverlap", type="Flag",
130  doc="set if this source was not completely inside the generated image"
131  )
132 
133  def buildSkyMap(self, butler):
134  """Build the skymap for the mock dataset."""
135  return self.makeSkyMap.run(butler.dataRef(self.config.coaddName + "Coadd_skyMap")).skyMap
136 
137  def buildTruthCatalog(self, butler=None, skyMap=None, tract=0):
138  """Create and save (if butler is not None) a truth catalog containing all the mock objects.
139 
140  Must be run after buildSkyMap.
141 
142  Most of the work is delegated to the mockObject subtask.
143  """
144  if skyMap is None:
145  skyMap = butler.get(self.config.coaddName + "Coadd_skyMap")
146  catalog = self.mockObject.run(tractInfo=skyMap[tract])
147  if butler is not None:
148  butler.put(catalog, "truth", tract=tract)
149  return catalog
150 
151  def buildObservationCatalog(self, butler=None, skyMap=None, tract=0, camera=None):
152  """Create and save (if butler is not None) an ExposureCatalog of simulated observations,
153  containing the Psfs, Wcss, Calibs, etc. of the calexps to be simulated.
154 
155  Must be run after buildSkyMap.
156 
157  Most of the work is delegated to the mockObservation subtask.
158  """
159  if skyMap is None:
160  skyMap = butler.get(self.config.coaddName + "Coadd_skyMap")
161  if camera is None:
162  camera = butler.get("camera")
163  catalog = self.mockObservation.run(butler=butler,
164  n=self.config.nObservations, camera=camera,
165  tractInfo=skyMap[tract])
166  if butler is not None:
167  butler.put(catalog, "observations", tract=tract)
168  return catalog
169 
170  def buildInputImages(self, butler, obsCatalog=None, truthCatalog=None, tract=0):
171  """Use the truth catalog and observation catalog to create and save (if butler is not None)
172  mock calexps and an ExposureCatalog ('simsrc') that contains information about which objects
173  appear partially or fully in each exposure.
174 
175  Must be run after buildTruthCatalog and buildObservationCatalog.
176  """
177  if obsCatalog is None:
178  obsCatalog = butler.get("observations", tract=tract)
179  if truthCatalog is None:
180  truthCatalog = butler.get("truth", tract=tract)
181  ccdKey = obsCatalog.getSchema().find("ccd").key
182  visitKey = obsCatalog.getSchema().find("visit").key
183  simSrcCatalog = lsst.afw.table.SimpleCatalog(self.schema)
184  for obsRecord in obsCatalog:
185  ccd = obsRecord.getI(ccdKey)
186  visit = obsRecord.getI(visitKey)
187  self.log.info("Generating image for visit={visit}, ccd={ccd}".format(ccd=ccd, visit=visit))
188  exposure = lsst.afw.image.ExposureF(obsRecord.getBBox())
189  exposure.setCalib(obsRecord.getCalib())
190  exposure.setWcs(obsRecord.getWcs())
191  exposure.setPsf(obsRecord.getPsf())
192  for truthRecord in truthCatalog:
193  status = self.mockObject.drawSource(truthRecord, exposure, buffer=self.config.edgeBuffer)
194  if status:
195  simSrcRecord = simSrcCatalog.addNew()
196  simSrcRecord.setCoord(truthRecord.getCoord())
197  simSrcRecord.setL(self.objectIdKey, truthRecord.getId())
198  simSrcRecord.setL(self.exposureIdKey, obsRecord.getId())
199  simSrcRecord.setFlag(self.centroidInBBoxKey, obsRecord.contains(truthRecord.getCoord()))
200  simSrcRecord.setFlag(self.partialOverlapKey, status == 1)
201  self.log.info(" added object {id}".format(id=truthRecord.getId()))
202  exposure.getMaskedImage().getVariance().set(1.0)
203  if butler is not None:
204  butler.put(exposure, "calexp", ccd=ccd, visit=visit)
205  if butler is not None:
206  butler.put(simSrcCatalog, "simsrc", tract=tract)
207  return simSrcCatalog
208 
209  def buildAllInputs(self, butler):
210  """Convenience function that calls buildSkyMap, buildObservationCatalog, buildTruthCatalog,
211  and buildInputImages.
212  """
213  skyMap = self.buildSkyMap(butler)
214  observations = self.buildObservationCatalog(butler, skyMap=skyMap)
215  truth = self.buildTruthCatalog(butler, skyMap=skyMap)
216  simSrcCatalog = self.buildInputImages(butler, obsCatalog=observations, truthCatalog=truth)
217 
218  def makeCoaddTask(self, cls):
219  """Helper function to create a Coadd task with configuration appropriate for the simulations.
220 
221  MockCoaddTask does not include MakeCoaddTempExpTask or AssembleCoaddTask as subtasks, because
222  we want explicit control over their configs, rather than leaving this up to the user.
223  However, we have to install our own SelectImages task for both of these, so it made sense
224  to have a single method that would create one of these two tasks, set the config values we
225  want, and install the custom SelectImagesTask.
226  """
227  config = cls.ConfigClass()
228  config.coaddName = self.config.coaddName
229  config.select.retarget(MockSelectImagesTask)
230  if cls == MakeCoaddTempExpTask:
231  config.bgSubtracted = True
232  config.doPsfMatch = False
233  elif cls == AssembleCoaddTask:
234  config.doMatchBackgrounds = False
235  return cls(config)
236 
237  def iterPatchRefs(self, butler, tractInfo):
238  """Generator that iterates over the patches in a tract, yielding dataRefs.
239  """
240  nPatchX, nPatchY = tractInfo.getNumPatches()
241  for iPatchX in range(nPatchX):
242  for iPatchY in range(nPatchY):
243  patchRef = butler.dataRef(self.config.coaddName + "Coadd",
244  tract=tractInfo.getId(), patch="%d,%d" % (iPatchX,iPatchY),
245  filter='r')
246  yield patchRef
247 
248  def buildCoadd(self, butler, skyMap=None, tract=0):
249  """Run the coadd tasks (MakeCoaddTempExp and AssembleCoadd) on the mock data.
250 
251  Must be run after buildInputImages.
252  """
253  if skyMap is None:
254  skyMap = butler.get(self.config.coaddName + "Coadd_skyMap")
255  tractInfo = skyMap[tract]
256  makeCoaddTempExpTask = self.makeCoaddTask(MakeCoaddTempExpTask)
257  assembleCoaddTask = self.makeCoaddTask(AssembleCoaddTask)
258  for patchRef in self.iterPatchRefs(butler, tractInfo):
259  makeCoaddTempExpTask.run(patchRef)
260  for patchRef in self.iterPatchRefs(butler, tractInfo):
261  assembleCoaddTask.run(patchRef)
262 
263  def buildMockCoadd(self, butler, truthCatalog=None, skyMap=None, tract=0):
264  """Directly create a simulation of the coadd, using the CoaddPsf of the coadd exposure
265  and the truth catalog.
266 
267  Must be run after buildCoadd.
268  """
269  if truthCatalog is None:
270  truthCatalog = butler.get("truth", tract=tract)
271  if skyMap is None:
272  skyMap = butler.get(self.config.coaddName + "Coadd_skyMap")
273  tractInfo = skyMap[tract]
274  tractWcs = tractInfo.getWcs()
275  for patchRef in self.iterPatchRefs(butler, tractInfo):
276  exposure = patchRef.get(self.config.coaddName + "Coadd")
277  exposure.getMaskedImage().getImage().set(0.0)
279  exposure.getInfo().getCoaddInputs().ccds, exposure.getWcs()
280  )
281  exposure.setPsf(coaddPsf)
282  for truthRecord in truthCatalog:
283  self.mockObject.drawSource(truthRecord, exposure, buffer=0)
284  patchRef.put(exposure, self.config.coaddName + "Coadd_mock")
285 
286 def run(root):
287  """Convenience function to create and run MockCoaddTask with default settings.
288  """
289  from .simpleMapper import makeDataRepo
290  butler = makeDataRepo(root=root)
291  task = MockCoaddTask()
292  task.buildAllInputs(butler)
293  task.buildCoadd(butler)
294  task.buildMockCoadd(butler)
static Schema makeMinimalSchema()
Return a minimal schema for Simple tables and records.
Definition: Simple.h:121
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
CoaddPsf is the Psf derived to be used for non-PSF-matched Coadd images.
Definition: CoaddPsf.h:45