LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
makeCoaddTempExp.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 import numpy
23 
24 import lsst.pex.config as pexConfig
25 import lsst.daf.persistence as dafPersist
26 import lsst.afw.image as afwImage
27 import lsst.coadd.utils as coaddUtils
28 import lsst.pipe.base as pipeBase
30 import lsst.log as log
31 import lsst.utils as utils
32 from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig
33 from .coaddBase import CoaddBaseTask, makeSkyInfo
34 from .warpAndPsfMatch import WarpAndPsfMatchTask
35 from .coaddHelpers import groupPatchExposures, getGroupDataRef
36 
37 __all__ = ["MakeCoaddTempExpTask", "MakeWarpTask", "MakeWarpConfig"]
38 
39 
40 class MissingExposureError(Exception):
41  """Raised when data cannot be retrieved for an exposure.
42  When processing patches, sometimes one exposure is missing; this lets us
43  distinguish bewteen that case, and other errors.
44  """
45  pass
46 
47 
48 class MakeCoaddTempExpConfig(CoaddBaseTask.ConfigClass):
49  """Config for MakeCoaddTempExpTask
50  """
51  warpAndPsfMatch = pexConfig.ConfigurableField(
52  target=WarpAndPsfMatchTask,
53  doc="Task to warp and PSF-match calexp",
54  )
55  doWrite = pexConfig.Field(
56  doc="persist <coaddName>Coadd_<warpType>Warp",
57  dtype=bool,
58  default=True,
59  )
60  bgSubtracted = pexConfig.Field(
61  doc="Work with a background subtracted calexp?",
62  dtype=bool,
63  default=True,
64  )
65  coaddPsf = pexConfig.ConfigField(
66  doc="Configuration for CoaddPsf",
67  dtype=CoaddPsfConfig,
68  )
69  makeDirect = pexConfig.Field(
70  doc="Make direct Warp/Coadds",
71  dtype=bool,
72  default=True,
73  )
74  makePsfMatched = pexConfig.Field(
75  doc="Make Psf-Matched Warp/Coadd?",
76  dtype=bool,
77  default=False,
78  )
79 
80  doWriteEmptyWarps = pexConfig.Field(
81  dtype=bool,
82  default=False,
83  doc="Write out warps even if they are empty"
84  )
85 
86  hasFakes = pexConfig.Field(
87  doc="Should be set to True if fake sources have been inserted into the input data.",
88  dtype=bool,
89  default=False,
90  )
91  doApplySkyCorr = pexConfig.Field(dtype=bool, default=False, doc="Apply sky correction?")
92 
93  def validate(self):
94  CoaddBaseTask.ConfigClass.validate(self)
95  if not self.makePsfMatched and not self.makeDirect:
96  raise RuntimeError("At least one of config.makePsfMatched and config.makeDirect must be True")
97  if self.doPsfMatch:
98  # Backwards compatibility.
99  log.warn("Config doPsfMatch deprecated. Setting makePsfMatched=True and makeDirect=False")
100  self.makePsfMatched = True
101  self.makeDirect = False
102 
103  def setDefaults(self):
104  CoaddBaseTask.ConfigClass.setDefaults(self)
105  self.warpAndPsfMatch.psfMatch.kernel.active.kernelSize = self.matchingKernelSize
106 
107 
113 
114 
116  r"""!Warp and optionally PSF-Match calexps onto an a common projection.
117 
118  @anchor MakeCoaddTempExpTask_
119 
120  @section pipe_tasks_makeCoaddTempExp_Contents Contents
121 
122  - @ref pipe_tasks_makeCoaddTempExp_Purpose
123  - @ref pipe_tasks_makeCoaddTempExp_Initialize
124  - @ref pipe_tasks_makeCoaddTempExp_IO
125  - @ref pipe_tasks_makeCoaddTempExp_Config
126  - @ref pipe_tasks_makeCoaddTempExp_Debug
127  - @ref pipe_tasks_makeCoaddTempExp_Example
128 
129  @section pipe_tasks_makeCoaddTempExp_Purpose Description
130 
131  Warp and optionally PSF-Match calexps onto a common projection, by
132  performing the following operations:
133  - Group calexps by visit/run
134  - For each visit, generate a Warp by calling method @ref makeTempExp.
135  makeTempExp loops over the visit's calexps calling @ref WarpAndPsfMatch
136  on each visit
137 
138  The result is a `directWarp` (and/or optionally a `psfMatchedWarp`).
139 
140  @section pipe_tasks_makeCoaddTempExp_Initialize Task Initialization
141 
142  @copydoc \_\_init\_\_
143 
144  This task has one special keyword argument: passing reuse=True will cause
145  the task to skip the creation of warps that are already present in the
146  output repositories.
147 
148  @section pipe_tasks_makeCoaddTempExp_IO Invoking the Task
149 
150  This task is primarily designed to be run from the command line.
151 
152  The main method is `runDataRef`, which takes a single butler data reference for the patch(es)
153  to process.
154 
155  @copydoc run
156 
157  WarpType identifies the types of convolutions applied to Warps (previously CoaddTempExps).
158  Only two types are available: direct (for regular Warps/Coadds) and psfMatched
159  (for Warps/Coadds with homogenized PSFs). We expect to add a third type, likelihood,
160  for generating likelihood Coadds with Warps that have been correlated with their own PSF.
161 
162  @section pipe_tasks_makeCoaddTempExp_Config Configuration parameters
163 
164  See @ref MakeCoaddTempExpConfig and parameters inherited from
165  @link lsst.pipe.tasks.coaddBase.CoaddBaseConfig CoaddBaseConfig @endlink
166 
167  @subsection pipe_tasks_MakeCoaddTempExp_psfMatching Guide to PSF-Matching Configs
168 
169  To make `psfMatchedWarps`, select `config.makePsfMatched=True`. The subtask
170  @link lsst.ip.diffim.modelPsfMatch.ModelPsfMatchTask ModelPsfMatchTask @endlink
171  is responsible for the PSF-Matching, and its config is accessed via `config.warpAndPsfMatch.psfMatch`.
172  The optimal configuration depends on aspects of dataset: the pixel scale, average PSF FWHM and
173  dimensions of the PSF kernel. These configs include the requested model PSF, the matching kernel size,
174  padding of the science PSF thumbnail and spatial sampling frequency of the PSF.
175 
176  *Config Guidelines*: The user must specify the size of the model PSF to which to match by setting
177  `config.modelPsf.defaultFwhm` in units of pixels. The appropriate values depends on science case.
178  In general, for a set of input images, this config should equal the FWHM of the visit
179  with the worst seeing. The smallest it should be set to is the median FWHM. The defaults
180  of the other config options offer a reasonable starting point.
181  The following list presents the most common problems that arise from a misconfigured
182  @link lsst.ip.diffim.modelPsfMatch.ModelPsfMatchTask ModelPsfMatchTask @endlink
183  and corresponding solutions. All assume the default Alard-Lupton kernel, with configs accessed via
184  ```config.warpAndPsfMatch.psfMatch.kernel['AL']```. Each item in the list is formatted as:
185  Problem: Explanation. *Solution*
186 
187  *Troublshooting PSF-Matching Configuration:*
188  - Matched PSFs look boxy: The matching kernel is too small. _Increase the matching kernel size.
189  For example:_
190 
191  config.warpAndPsfMatch.psfMatch.kernel['AL'].kernelSize=27 # default 21
192 
193  Note that increasing the kernel size also increases runtime.
194  - Matched PSFs look ugly (dipoles, quadropoles, donuts): unable to find good solution
195  for matching kernel. _Provide the matcher with more data by either increasing
196  the spatial sampling by decreasing the spatial cell size,_
197 
198  config.warpAndPsfMatch.psfMatch.kernel['AL'].sizeCellX = 64 # default 128
199  config.warpAndPsfMatch.psfMatch.kernel['AL'].sizeCellY = 64 # default 128
200 
201  _or increasing the padding around the Science PSF, for example:_
202 
203  config.warpAndPsfMatch.psfMatch.autoPadPsfTo=1.6 # default 1.4
204 
205  Increasing `autoPadPsfTo` increases the minimum ratio of input PSF dimensions to the
206  matching kernel dimensions, thus increasing the number of pixels available to fit
207  after convolving the PSF with the matching kernel.
208  Optionally, for debugging the effects of padding, the level of padding may be manually
209  controlled by setting turning off the automatic padding and setting the number
210  of pixels by which to pad the PSF:
211 
212  config.warpAndPsfMatch.psfMatch.doAutoPadPsf = False # default True
213  config.warpAndPsfMatch.psfMatch.padPsfBy = 6 # pixels. default 0
214 
215  - Deconvolution: Matching a large PSF to a smaller PSF produces
216  a telltale noise pattern which looks like ripples or a brain.
217  _Increase the size of the requested model PSF. For example:_
218 
219  config.modelPsf.defaultFwhm = 11 # Gaussian sigma in units of pixels.
220 
221  - High frequency (sometimes checkered) noise: The matching basis functions are too small.
222  _Increase the width of the Gaussian basis functions. For example:_
223 
224  config.warpAndPsfMatch.psfMatch.kernel['AL'].alardSigGauss=[1.5, 3.0, 6.0]
225  # from default [0.7, 1.5, 3.0]
226 
227 
228  @section pipe_tasks_makeCoaddTempExp_Debug Debug variables
229 
230  MakeCoaddTempExpTask has no debug output, but its subtasks do.
231 
232  @section pipe_tasks_makeCoaddTempExp_Example A complete example of using MakeCoaddTempExpTask
233 
234  This example uses the package ci_hsc to show how MakeCoaddTempExp fits
235  into the larger Data Release Processing.
236  Set up by running:
237 
238  setup ci_hsc
239  cd $CI_HSC_DIR
240  # if not built already:
241  python $(which scons) # this will take a while
242 
243  The following assumes that `processCcd.py` and `makeSkyMap.py` have previously been run
244  (e.g. by building `ci_hsc` above) to generate a repository of calexps and an
245  output respository with the desired SkyMap. The command,
246 
247  makeCoaddTempExp.py $CI_HSC_DIR/DATA --rerun ci_hsc \
248  --id patch=5,4 tract=0 filter=HSC-I \
249  --selectId visit=903988 ccd=16 --selectId visit=903988 ccd=17 \
250  --selectId visit=903988 ccd=23 --selectId visit=903988 ccd=24 \
251  --config doApplyUberCal=False makePsfMatched=True modelPsf.defaultFwhm=11
252 
253  writes a direct and PSF-Matched Warp to
254  - `$CI_HSC_DIR/DATA/rerun/ci_hsc/deepCoadd/HSC-I/0/5,4/warp-HSC-I-0-5,4-903988.fits` and
255  - `$CI_HSC_DIR/DATA/rerun/ci_hsc/deepCoadd/HSC-I/0/5,4/psfMatchedWarp-HSC-I-0-5,4-903988.fits`
256  respectively.
257 
258  @note PSF-Matching in this particular dataset would benefit from adding
259  `--configfile ./matchingConfig.py` to
260  the command line arguments where `matchingConfig.py` is defined by:
261 
262  echo "
263  config.warpAndPsfMatch.psfMatch.kernel['AL'].kernelSize=27
264  config.warpAndPsfMatch.psfMatch.kernel['AL'].alardSigGauss=[1.5, 3.0, 6.0]" > matchingConfig.py
265 
266 
267  Add the option `--help` to see more options.
268  """
269  ConfigClass = MakeCoaddTempExpConfig
270  _DefaultName = "makeCoaddTempExp"
271 
272  def __init__(self, reuse=False, **kwargs):
273  CoaddBaseTask.__init__(self, **kwargs)
274  self.reuse = reuse
275  self.makeSubtask("warpAndPsfMatch")
276  if self.config.hasFakes:
277  self.calexpType = "fakes_calexp"
278  else:
279  self.calexpType = "calexp"
280 
281  @pipeBase.timeMethod
282  def runDataRef(self, patchRef, selectDataList=[]):
283  """!Produce <coaddName>Coadd_<warpType>Warp images by warping and optionally PSF-matching.
284 
285  @param[in] patchRef: data reference for sky map patch. Must include keys "tract", "patch",
286  plus the camera-specific filter key (e.g. "filter" or "band")
287  @return: dataRefList: a list of data references for the new <coaddName>Coadd_directWarps
288  if direct or both warp types are requested and <coaddName>Coadd_psfMatchedWarps if only psfMatched
289  warps are requested.
290 
291  @warning: this task assumes that all exposures in a warp (coaddTempExp) have the same filter.
292 
293  @warning: this task sets the PhotoCalib of the coaddTempExp to the PhotoCalib of the first calexp
294  with any good pixels in the patch. For a mosaic camera the resulting PhotoCalib should be ignored
295  (assembleCoadd should determine zeropoint scaling without referring to it).
296  """
297  skyInfo = self.getSkyInfo(patchRef)
298 
299  # DataRefs to return are of type *_directWarp unless only *_psfMatchedWarp requested
300  if self.config.makePsfMatched and not self.config.makeDirect:
301  primaryWarpDataset = self.getTempExpDatasetName("psfMatched")
302  else:
303  primaryWarpDataset = self.getTempExpDatasetName("direct")
304 
305  calExpRefList = self.selectExposures(patchRef, skyInfo, selectDataList=selectDataList)
306 
307  if len(calExpRefList) == 0:
308  self.log.warn("No exposures to coadd for patch %s", patchRef.dataId)
309  return None
310  self.log.info("Selected %d calexps for patch %s", len(calExpRefList), patchRef.dataId)
311  calExpRefList = [calExpRef for calExpRef in calExpRefList if calExpRef.datasetExists(self.calexpType)]
312  self.log.info("Processing %d existing calexps for patch %s", len(calExpRefList), patchRef.dataId)
313 
314  groupData = groupPatchExposures(patchRef, calExpRefList, self.getCoaddDatasetName(),
315  primaryWarpDataset)
316  self.log.info("Processing %d warp exposures for patch %s", len(groupData.groups), patchRef.dataId)
317 
318  dataRefList = []
319  for i, (tempExpTuple, calexpRefList) in enumerate(groupData.groups.items()):
320  tempExpRef = getGroupDataRef(patchRef.getButler(), primaryWarpDataset,
321  tempExpTuple, groupData.keys)
322  if self.reuse and tempExpRef.datasetExists(datasetType=primaryWarpDataset, write=True):
323  self.log.info("Skipping makeCoaddTempExp for %s; output already exists.", tempExpRef.dataId)
324  dataRefList.append(tempExpRef)
325  continue
326  self.log.info("Processing Warp %d/%d: id=%s", i, len(groupData.groups), tempExpRef.dataId)
327 
328  # TODO: mappers should define a way to go from the "grouping keys" to a numeric ID (#2776).
329  # For now, we try to get a long integer "visit" key, and if we can't, we just use the index
330  # of the visit in the list.
331  try:
332  visitId = int(tempExpRef.dataId["visit"])
333  except (KeyError, ValueError):
334  visitId = i
335 
336  calExpList = []
337  ccdIdList = []
338  dataIdList = []
339 
340  for calExpInd, calExpRef in enumerate(calexpRefList):
341  self.log.info("Reading calexp %s of %s for Warp id=%s", calExpInd+1, len(calexpRefList),
342  calExpRef.dataId)
343  try:
344  ccdId = calExpRef.get("ccdExposureId", immediate=True)
345  except Exception:
346  ccdId = calExpInd
347  try:
348  # We augment the dataRef here with the tract, which is harmless for loading things
349  # like calexps that don't need the tract, and necessary for meas_mosaic outputs,
350  # which do.
351  calExpRef = calExpRef.butlerSubset.butler.dataRef(self.calexpType,
352  dataId=calExpRef.dataId,
353  tract=skyInfo.tractInfo.getId())
354  calExp = self.getCalibratedExposure(calExpRef, bgSubtracted=self.config.bgSubtracted)
355  except Exception as e:
356  self.log.warn("Calexp %s not found; skipping it: %s", calExpRef.dataId, e)
357  continue
358 
359  if self.config.doApplySkyCorr:
360  self.applySkyCorr(calExpRef, calExp)
361 
362  calExpList.append(calExp)
363  ccdIdList.append(ccdId)
364  dataIdList.append(calExpRef.dataId)
365 
366  exps = self.run(calExpList, ccdIdList, skyInfo, visitId, dataIdList).exposures
367 
368  if any(exps.values()):
369  dataRefList.append(tempExpRef)
370  else:
371  self.log.warn("Warp %s could not be created", tempExpRef.dataId)
372 
373  if self.config.doWrite:
374  for (warpType, exposure) in exps.items(): # compatible w/ Py3
375  if exposure is not None:
376  self.log.info("Persisting %s" % self.getTempExpDatasetName(warpType))
377  tempExpRef.put(exposure, self.getTempExpDatasetName(warpType))
378 
379  return dataRefList
380 
381  def run(self, calExpList, ccdIdList, skyInfo, visitId=0, dataIdList=None):
382  """Create a Warp from inputs
383 
384  We iterate over the multiple calexps in a single exposure to construct
385  the warp (previously called a coaddTempExp) of that exposure to the
386  supplied tract/patch.
387 
388  Pixels that receive no pixels are set to NAN; this is not correct
389  (violates LSST algorithms group policy), but will be fixed up by
390  interpolating after the coaddition.
391 
392  @param calexpRefList: List of data references for calexps that (may)
393  overlap the patch of interest
394  @param skyInfo: Struct from CoaddBaseTask.getSkyInfo() with geometric
395  information about the patch
396  @param visitId: integer identifier for visit, for the table that will
397  produce the CoaddPsf
398  @return a pipeBase Struct containing:
399  - exposures: a dictionary containing the warps requested:
400  "direct": direct warp if config.makeDirect
401  "psfMatched": PSF-matched warp if config.makePsfMatched
402  """
403  warpTypeList = self.getWarpTypeList()
404 
405  totGoodPix = {warpType: 0 for warpType in warpTypeList}
406  didSetMetadata = {warpType: False for warpType in warpTypeList}
407  coaddTempExps = {warpType: self._prepareEmptyExposure(skyInfo) for warpType in warpTypeList}
408  inputRecorder = {warpType: self.inputRecorder.makeCoaddTempExpRecorder(visitId, len(calExpList))
409  for warpType in warpTypeList}
410 
411  modelPsf = self.config.modelPsf.apply() if self.config.makePsfMatched else None
412  if dataIdList is None:
413  dataIdList = ccdIdList
414 
415  for calExpInd, (calExp, ccdId, dataId) in enumerate(zip(calExpList, ccdIdList, dataIdList)):
416  self.log.info("Processing calexp %d of %d for this Warp: id=%s",
417  calExpInd+1, len(calExpList), dataId)
418 
419  try:
420  warpedAndMatched = self.warpAndPsfMatch.run(calExp, modelPsf=modelPsf,
421  wcs=skyInfo.wcs, maxBBox=skyInfo.bbox,
422  makeDirect=self.config.makeDirect,
423  makePsfMatched=self.config.makePsfMatched)
424  except Exception as e:
425  self.log.warn("WarpAndPsfMatch failed for calexp %s; skipping it: %s", dataId, e)
426  continue
427  try:
428  numGoodPix = {warpType: 0 for warpType in warpTypeList}
429  for warpType in warpTypeList:
430  exposure = warpedAndMatched.getDict()[warpType]
431  if exposure is None:
432  continue
433  coaddTempExp = coaddTempExps[warpType]
434  if didSetMetadata[warpType]:
435  mimg = exposure.getMaskedImage()
436  mimg *= (coaddTempExp.getPhotoCalib().getInstFluxAtZeroMagnitude() /
437  exposure.getPhotoCalib().getInstFluxAtZeroMagnitude())
438  del mimg
439  numGoodPix[warpType] = coaddUtils.copyGoodPixels(
440  coaddTempExp.getMaskedImage(), exposure.getMaskedImage(), self.getBadPixelMask())
441  totGoodPix[warpType] += numGoodPix[warpType]
442  self.log.debug("Calexp %s has %d good pixels in this patch (%.1f%%) for %s",
443  dataId, numGoodPix[warpType],
444  100.0*numGoodPix[warpType]/skyInfo.bbox.getArea(), warpType)
445  if numGoodPix[warpType] > 0 and not didSetMetadata[warpType]:
446  coaddTempExp.setPhotoCalib(exposure.getPhotoCalib())
447  coaddTempExp.setFilter(exposure.getFilter())
448  coaddTempExp.getInfo().setVisitInfo(exposure.getInfo().getVisitInfo())
449  # PSF replaced with CoaddPsf after loop if and only if creating direct warp
450  coaddTempExp.setPsf(exposure.getPsf())
451  didSetMetadata[warpType] = True
452 
453  # Need inputRecorder for CoaddApCorrMap for both direct and PSF-matched
454  inputRecorder[warpType].addCalExp(calExp, ccdId, numGoodPix[warpType])
455 
456  except Exception as e:
457  self.log.warn("Error processing calexp %s; skipping it: %s", dataId, e)
458  continue
459 
460  for warpType in warpTypeList:
461  self.log.info("%sWarp has %d good pixels (%.1f%%)",
462  warpType, totGoodPix[warpType], 100.0*totGoodPix[warpType]/skyInfo.bbox.getArea())
463 
464  if totGoodPix[warpType] > 0 and didSetMetadata[warpType]:
465  inputRecorder[warpType].finish(coaddTempExps[warpType], totGoodPix[warpType])
466  if warpType == "direct":
467  coaddTempExps[warpType].setPsf(
468  CoaddPsf(inputRecorder[warpType].coaddInputs.ccds, skyInfo.wcs,
469  self.config.coaddPsf.makeControl()))
470  else:
471  if not self.config.doWriteEmptyWarps:
472  # No good pixels. Exposure still empty
473  coaddTempExps[warpType] = None
474 
475  result = pipeBase.Struct(exposures=coaddTempExps)
476  return result
477 
478  def getCalibratedExposure(self, dataRef, bgSubtracted):
479  """Return one calibrated Exposure, possibly with an updated SkyWcs.
480 
481  @param[in] dataRef a sensor-level data reference
482  @param[in] bgSubtracted return calexp with background subtracted? If False get the
483  calexp's background background model and add it to the calexp.
484  @return calibrated exposure
485 
486  @raises MissingExposureError If data for the exposure is not available.
487 
488  If config.doApplyUberCal, the exposure will be photometrically
489  calibrated via the `jointcal_photoCalib` dataset and have its SkyWcs
490  updated to the `jointcal_wcs`, otherwise it will be calibrated via the
491  Exposure's own PhotoCalib and have the original SkyWcs.
492  """
493  try:
494  exposure = dataRef.get(self.calexpType, immediate=True)
495  except dafPersist.NoResults as e:
496  raise MissingExposureError('Exposure not found: %s ' % str(e)) from e
497 
498  if not bgSubtracted:
499  background = dataRef.get("calexpBackground", immediate=True)
500  mi = exposure.getMaskedImage()
501  mi += background.getImage()
502  del mi
503 
504  if self.config.doApplyUberCal:
505  if self.config.useMeasMosaic:
506  from lsst.meas.mosaic import applyMosaicResultsExposure
507  # NOTE: this changes exposure in-place, updating its Calib and Wcs.
508  # Save the calibration error, as it gets overwritten with zero.
509  calibrationErr = exposure.getPhotoCalib().getCalibrationErr()
510  try:
511  applyMosaicResultsExposure(dataRef, calexp=exposure)
512  except dafPersist.NoResults as e:
513  raise MissingExposureError('Mosaic calibration not found: %s ' % str(e)) from e
514  photoCalib = afwImage.PhotoCalib(exposure.getPhotoCalib().getCalibrationMean(),
515  calibrationErr,
516  exposure.getBBox())
517  else:
518  photoCalib = dataRef.get("jointcal_photoCalib")
519  skyWcs = dataRef.get("jointcal_wcs")
520  exposure.setWcs(skyWcs)
521  else:
522  photoCalib = exposure.getPhotoCalib()
523 
524  exposure.maskedImage = photoCalib.calibrateImage(exposure.maskedImage,
525  includeScaleUncertainty=self.config.includeCalibVar)
526  exposure.maskedImage /= photoCalib.getCalibrationMean()
527  exposure.setPhotoCalib(photoCalib)
528  # TODO: The images will have a calibration of 1.0 everywhere once RFC-545 is implemented.
529  # exposure.setCalib(afwImage.Calib(1.0))
530  return exposure
531 
532  @staticmethod
533  def _prepareEmptyExposure(skyInfo):
534  """Produce an empty exposure for a given patch"""
535  exp = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs)
536  exp.getMaskedImage().set(numpy.nan, afwImage.Mask
537  .getPlaneBitMask("NO_DATA"), numpy.inf)
538  return exp
539 
540  def getWarpTypeList(self):
541  """Return list of requested warp types per the config.
542  """
543  warpTypeList = []
544  if self.config.makeDirect:
545  warpTypeList.append("direct")
546  if self.config.makePsfMatched:
547  warpTypeList.append("psfMatched")
548  return warpTypeList
549 
550  def applySkyCorr(self, dataRef, calexp):
551  """Apply correction to the sky background level
552 
553  Sky corrections can be generated with the 'skyCorrection.py'
554  executable in pipe_drivers. Because the sky model used by that
555  code extends over the entire focal plane, this can produce
556  better sky subtraction.
557 
558  The calexp is updated in-place.
559 
560  Parameters
561  ----------
562  dataRef : `lsst.daf.persistence.ButlerDataRef`
563  Data reference for calexp.
564  calexp : `lsst.afw.image.Exposure` or `lsst.afw.image.MaskedImage`
565  Calibrated exposure.
566  """
567  bg = dataRef.get("skyCorr")
568  if isinstance(calexp, afwImage.Exposure):
569  calexp = calexp.getMaskedImage()
570  calexp -= bg.getImage()
571 
572 
573 class MakeWarpConnections(pipeBase.PipelineTaskConnections,
574  dimensions=("tract", "patch", "skymap", "instrument", "visit"),
575  defaultTemplates={"coaddName": "deep"}):
576  calExpList = cT.Input(
577  doc="Input exposures to be resampled and optionally PSF-matched onto a SkyMap projection/patch",
578  name="calexp",
579  storageClass="ExposureF",
580  dimensions=("instrument", "visit", "detector"),
581  multiple=True,
582  )
583  backgroundList = cT.Input(
584  doc="Input backgrounds to be added back into the calexp if bgSubtracted=False",
585  name="calexpBackground",
586  storageClass="Background",
587  dimensions=("instrument", "visit", "detector"),
588  multiple=True,
589  )
590  skyCorrList = cT.Input(
591  doc="Input Sky Correction to be subtracted from the calexp if doApplySkyCorr=True",
592  name="skyCorr",
593  storageClass="Background",
594  dimensions=("instrument", "visit", "detector"),
595  multiple=True,
596  )
597  skyMap = cT.Input(
598  doc="Input definition of geometry/bbox and projection/wcs for warped exposures",
599  name="{coaddName}Coadd_skyMap",
600  storageClass="SkyMap",
601  dimensions=("skymap",),
602  )
603  direct = cT.Output(
604  doc=("Output direct warped exposure (previously called CoaddTempExp), produced by resampling ",
605  "calexps onto the skyMap patch geometry."),
606  name="{coaddName}Coadd_directWarp",
607  storageClass="ExposureF",
608  dimensions=("tract", "patch", "skymap", "visit", "instrument"),
609  )
610  psfMatched = cT.Output(
611  doc=("Output PSF-Matched warped exposure (previously called CoaddTempExp), produced by resampling ",
612  "calexps onto the skyMap patch geometry and PSF-matching to a model PSF."),
613  name="{coaddName}Coadd_psfMatchedWarp",
614  storageClass="ExposureF",
615  dimensions=("tract", "patch", "skymap", "visit", "instrument"),
616  )
617 
618  def __init__(self, *, config=None):
619  super().__init__(config=config)
620  if config.bgSubtracted:
621  self.inputs.remove("backgroundList")
622  if not config.doApplySkyCorr:
623  self.inputs.remove("skyCorrList")
624  if not config.makeDirect:
625  self.outputs.remove("direct")
626  if not config.makePsfMatched:
627  self.outputs.remove("psfMatched")
628 
629 
630 class MakeWarpConfig(pipeBase.PipelineTaskConfig, MakeCoaddTempExpConfig,
631  pipelineConnections=MakeWarpConnections):
632 
633  def validate(self):
634  super().validate()
635  # TODO: Remove this constraint after DM-17062
636  if self.doApplyUberCal:
637  raise RuntimeError("Gen3 MakeWarpTask cannot apply meas_mosaic or jointcal results."
638  "Please set doApplyUbercal=False.")
639 
640 
641 class MakeWarpTask(MakeCoaddTempExpTask, pipeBase.PipelineTask):
642  """Warp and optionally PSF-Match calexps onto an a common projection
643 
644  First Draft of a Gen3 compatible MakeWarpTask which
645  currently does not handle doApplyUberCal=True.
646  """
647  ConfigClass = MakeWarpConfig
648  _DefaultName = "makeWarp"
649 
650  @utils.inheritDoc(pipeBase.PipelineTask)
651  def runQuantum(self, butlerQC, inputRefs, outputRefs):
652  """
653  Notes
654  ----
655  Construct warps for requested warp type for single epoch
656 
657  PipelineTask (Gen3) entry point to warp and optionally PSF-match
658  calexps. This method is analogous to `runDataRef`.
659  """
660  # Read in all inputs.
661  inputs = butlerQC.get(inputRefs)
662 
663  # Construct skyInfo expected by `run`. We remove the SkyMap itself
664  # from the dictionary so we can pass it as kwargs later.
665  skyMap = inputs.pop("skyMap")
666  quantumDataId = butlerQC.quantum.dataId
667  skyInfo = makeSkyInfo(skyMap, tractId=quantumDataId['tract'], patchId=quantumDataId['patch'])
668 
669  # Construct list of input DataIds expected by `run`
670  dataIdList = [ref.dataId for ref in inputRefs.calExpList]
671 
672  # Construct list of packed integer IDs expected by `run`
673  ccdIdList = [dataId.pack("visit_detector") for dataId in dataIdList]
674 
675  # Extract integer visitId requested by `run`
676  visits = [dataId['visit'] for dataId in dataIdList]
677  assert(all(visits[0] == visit for visit in visits))
678  visitId = visits[0]
679 
680  self.prepareCalibratedExposures(**inputs)
681  results = self.run(**inputs, visitId=visitId, ccdIdList=ccdIdList, dataIdList=dataIdList,
682  skyInfo=skyInfo)
683  if self.config.makeDirect:
684  butlerQC.put(results.exposures["direct"], outputRefs.direct)
685  if self.config.makePsfMatched:
686  butlerQC.put(results.exposures["psfMatched"], outputRefs.psfMatched)
687 
688  def prepareCalibratedExposures(self, calExpList, backgroundList=None, skyCorrList=None):
689  """Calibrate and add backgrounds to input calExpList in place
690 
691  TODO DM-17062: apply jointcal/meas_mosaic here
692 
693  Parameters
694  ----------
695  calExpList : `list` of `lsst.afw.image.Exposure`
696  Sequence of calexps to be modified in place
697  backgroundList : `list` of `lsst.afw.math.backgroundList`
698  Sequence of backgrounds to be added back in if bgSubtracted=False
699  skyCorrList : `list` of `lsst.afw.math.backgroundList`
700  Sequence of background corrections to be subtracted if doApplySkyCorr=True
701  """
702  backgroundList = len(calExpList)*[None] if backgroundList is None else backgroundList
703  skyCorrList = len(calExpList)*[None] if skyCorrList is None else skyCorrList
704  for calexp, background, skyCorr in zip(calExpList, backgroundList, skyCorrList):
705  mi = calexp.maskedImage
706  if not self.config.bgSubtracted:
707  mi += background.getImage()
708  if self.config.doApplySkyCorr:
709  mi -= skyCorr.getImage()
def getCoaddDatasetName(self, warpType="direct")
Definition: coaddBase.py:150
def getGroupDataRef(butler, datasetType, groupTuple, keys)
Definition: coaddHelpers.py:99
Base class for coaddition.
Definition: coaddBase.py:101
The photometric calibration of an exposure.
Definition: PhotoCalib.h:114
A class to contain the data, WCS, and other information needed to describe an image of the sky...
Definition: Exposure.h:72
def makeSkyInfo(skyMap, tractId, patchId)
Definition: coaddBase.py:249
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations...
daf::base::PropertySet * set
Definition: fits.cc:902
Warp and optionally PSF-Match calexps onto an a common projection.
bool any(CoordinateExpr< N > const &expr) noexcept
Return true if any elements are true.
Definition: Log.h:706
def getSkyInfo(self, patchRef)
Use getSkyinfo to return the skyMap, tract and patch information, wcs and the outer bbox of the patch...
Definition: coaddBase.py:134
def getTempExpDatasetName(self, warpType="direct")
Definition: coaddBase.py:164
bool all(CoordinateExpr< N > const &expr) noexcept
Return true if all elements are true.
def getBadPixelMask(self)
Convenience method to provide the bitmask from the mask plane names.
Definition: coaddBase.py:199
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
CoaddPsf is the Psf derived to be used for non-PSF-matched Coadd images.
Definition: CoaddPsf.h:58
def __init__(self, minimum, dataRange, Q)
def getCalibratedExposure(self, dataRef, bgSubtracted)
def selectExposures(self, patchRef, skyInfo=None, selectDataList=[])
Select exposures to coadd.
Definition: coaddBase.py:114
def runDataRef(self, patchRef, selectDataList=[])
Produce <coaddName>Coadd_<warpType>Warp images by warping and optionally PSF-matching.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
def run(self, calExpList, ccdIdList, skyInfo, visitId=0, dataIdList=None)
int copyGoodPixels(lsst::afw::image::MaskedImage< ImagePixelT, lsst::afw::image::MaskPixel, lsst::afw::image::VariancePixel > &destImage, lsst::afw::image::MaskedImage< ImagePixelT, lsst::afw::image::MaskPixel, lsst::afw::image::VariancePixel > const &srcImage, lsst::afw::image::MaskPixel const badPixelMask)
copy good pixels from one masked image to another
def groupPatchExposures(patchDataRef, calexpDataRefList, coaddDatasetType="deepCoadd", tempExpDatasetType="deepCoadd_directWarp")
Definition: coaddHelpers.py:60