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
makeDiscreteSkyMap.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 #
3 # LSST Data Management System
4 # Copyright 2008-2015 AURA/LSST.
5 #
6 # This product includes software developed by the
7 # LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20 # the GNU General Public License along with this program. If not,
21 # see <https://www.lsstcorp.org/LegalNotices/>.
22 #
23 import sys
24 import traceback
25 import lsst.geom
26 
27 import lsst.afw.image as afwImage
28 import lsst.afw.geom as afwGeom
29 import lsst.pex.config as pexConfig
30 import lsst.pipe.base as pipeBase
31 from lsst.skymap import DiscreteSkyMap, BaseSkyMap
32 
33 class MakeDiscreteSkyMapConfig(pexConfig.Config):
34  """Config for MakeDiscreteSkyMapTask
35  """
36  coaddName = pexConfig.Field(
37  doc = "coadd name, e.g. deep, goodSeeing, chiSquared",
38  dtype = str,
39  default = "deep",
40  )
41  skyMap = pexConfig.ConfigField(
42  dtype = BaseSkyMap.ConfigClass,
43  doc = "SkyMap configuration parameters, excluding position and radius"
44  )
45  borderSize = pexConfig.Field(
46  doc = "additional border added to the bounding box of the calexps, in degrees",
47  dtype = float,
48  default = 0.0
49  )
50  doAppend = pexConfig.Field(
51  doc = "append another tract to an existing DiscreteSkyMap on disk, if present?",
52  dtype = bool,
53  default = False
54  )
55  doWrite = pexConfig.Field(
56  doc = "persist the skyMap?",
57  dtype = bool,
58  default = True,
59  )
60 
61  def setDefaults(self):
62  self.skyMap.tractOverlap = 0.0
63 
64 class MakeDiscreteSkyMapRunner(pipeBase.TaskRunner):
65  """Run a task with all dataRefs at once, rather than one dataRef at a time.
66 
67  Call the run method of the task using two positional arguments:
68  - butler: data butler
69  - dataRefList: list of all dataRefs,
70  """
71  @staticmethod
72  def getTargetList(parsedCmd):
73  return [(parsedCmd.butler, parsedCmd.id.refList)]
74 
75  def __call__(self, args):
76  """
77  @param args Arguments for Task.run()
78 
79  @return:
80  - None if self.doReturnResults false
81  - A pipe_base Struct containing these fields if self.doReturnResults true:
82  - dataRef: the provided data reference
83  - metadata: task metadata after execution of run
84  - result: result returned by task run, or None if the task fails
85  """
86  butler, dataRefList = args
87  task = self.TaskClass(config=self.config, log=self.log)
88  result = None # in case the task fails
89  if self.doRaise:
90  result = task.run(butler, dataRefList)
91  else:
92  try:
93  result = task.run(butler, dataRefList)
94  except Exception, e:
95  task.log.fatal("Failed: %s" % e)
96  if not isinstance(e, pipeBase.TaskError):
97  traceback.print_exc(file=sys.stderr)
98  for dataRef in dataRefList:
99  task.writeMetadata(dataRef)
100 
101  if self.doReturnResults:
102  return pipeBase.Struct(
103  dataRefList = dataRefList,
104  metadata = task.metadata,
105  result = result,
106  )
107 
108 class MakeDiscreteSkyMapTask(pipeBase.CmdLineTask):
109  """!Make a DiscreteSkyMap in a repository, using the bounding box of a set of calexps.
110 
111  The command-line and run signatures and config are sufficiently different from MakeSkyMapTask
112  that we don't inherit from it, but it is a replacement, so we use the same config/metadata names.
113  """
114  ConfigClass = MakeDiscreteSkyMapConfig
115  _DefaultName = "makeDiscreteSkyMap"
116  RunnerClass = MakeDiscreteSkyMapRunner
117 
118  def __init__(self, **kwargs):
119  pipeBase.CmdLineTask.__init__(self, **kwargs)
120 
121  @pipeBase.timeMethod
122  def run(self, butler, dataRefList):
123  """!Make a skymap from the bounds of the given set of calexps.
124 
125  @param[in] butler data butler used to save the SkyMap
126  @param[in] dataRefList dataRefs of calexps used to determine the size and pointing of the SkyMap
127  @return a pipeBase Struct containing:
128  - skyMap: the constructed SkyMap
129  """
130  self.log.info("Extracting bounding boxes of %d images" % len(dataRefList))
131  points = []
132  for dataRef in dataRefList:
133  if not dataRef.datasetExists("calexp"):
134  self.log.warn("CalExp for %s does not exist: ignoring" % (dataRef.dataId,))
135  continue
136  md = dataRef.get("calexp_md", immediate=True)
137  wcs = afwImage.makeWcs(md)
138  # nb: don't need to worry about xy0 because Exposure saves Wcs with CRPIX shifted by (-x0, -y0).
139  boxI = afwGeom.Box2I(afwGeom.Point2I(0,0), afwGeom.Extent2I(md.get("NAXIS1"), md.get("NAXIS2")))
140  boxD = afwGeom.Box2D(boxI)
141  points.extend(tuple(wcs.pixelToSky(corner).getVector()) for corner in boxD.getCorners())
142  if len(points) == 0:
143  raise RuntimeError("No data found from which to compute convex hull")
144  self.log.info("Computing spherical convex hull")
145  polygon = lsst.geom.convexHull(points)
146  if polygon is None:
147  raise RuntimeError(
148  "Failed to compute convex hull of the vertices of all calexp bounding boxes; "
149  "they may not be hemispherical."
150  )
151  circle = polygon.getBoundingCircle()
152 
153  datasetName = self.config.coaddName + "Coadd_skyMap"
154 
155  skyMapConfig = DiscreteSkyMap.ConfigClass()
156  if self.config.doAppend and butler.datasetExists(datasetName):
157  oldSkyMap = butler.get(datasetName, immediate=True)
158  if not isinstance(oldSkyMap.config, DiscreteSkyMap.ConfigClass):
159  raise TypeError("Cannot append to existing non-discrete skymap")
160  compareLog = []
161  if not self.config.skyMap.compare(oldSkyMap.config, output=compareLog.append):
162  raise ValueError("Cannot append to existing skymap - configurations differ:", *compareLog)
163  skyMapConfig.raList.extend(oldSkyMap.config.raList)
164  skyMapConfig.decList.extend(oldSkyMap.config.decList)
165  skyMapConfig.radiusList.extend(oldSkyMap.config.radiusList)
166  skyMapConfig.update(**self.config.skyMap.toDict())
167  skyMapConfig.raList.append(circle.center[0])
168  skyMapConfig.decList.append(circle.center[1])
169  skyMapConfig.radiusList.append(circle.radius + self.config.borderSize)
170  skyMap = DiscreteSkyMap(skyMapConfig)
171 
172  for tractInfo in skyMap:
173  wcs = tractInfo.getWcs()
174  posBox = afwGeom.Box2D(tractInfo.getBBox())
175  pixelPosList = (
176  posBox.getMin(),
177  afwGeom.Point2D(posBox.getMaxX(), posBox.getMinY()),
178  posBox.getMax(),
179  afwGeom.Point2D(posBox.getMinX(), posBox.getMaxY()),
180  )
181  skyPosList = [wcs.pixelToSky(pos).getPosition(afwGeom.degrees) for pos in pixelPosList]
182  posStrList = ["(%0.3f, %0.3f)" % tuple(skyPos) for skyPos in skyPosList]
183  self.log.info("tract %s has corners %s (RA, Dec deg) and %s x %s patches" % \
184  (tractInfo.getId(), ", ".join(posStrList), \
185  tractInfo.getNumPatches()[0], tractInfo.getNumPatches()[1]))
186  if self.config.doWrite:
187  butler.put(skyMap, datasetName)
188  return pipeBase.Struct(
189  skyMap = skyMap
190  )
191 
192  def _getConfigName(self):
193  """Return None to disable saving config
194 
195  There's only one SkyMap per repository, so the config is redundant, and checking it means we can't
196  easily overwrite or append to an existing repository.
197  """
198  return None
199 
200  def _getMetadataName(self):
201  """Return None to disable saving metadata
202 
203  The metadata is not interesting, and by not saving it we can eliminate a dataset type.
204  """
205  return None
206 
An integer coordinate rectangle.
Definition: Box.h:53
Wcs::Ptr makeWcs(coord::Coord const &crval, geom::Point2D const &crpix, double CD11, double CD12, double CD21, double CD22)
Create a Wcs object from crval, crpix, CD, using CD elements (useful from python) ...
Definition: makeWcs.cc:141
def run
Make a skymap from the bounds of the given set of calexps.
A floating-point coordinate rectangle geometry.
Definition: Box.h:271
Make a DiscreteSkyMap in a repository, using the bounding box of a set of calexps.