LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
makeDiscreteSkyMap.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008-2015 AURA/LSST.
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 <https://www.lsstcorp.org/LegalNotices/>.
21 #
22 import sys
23 import traceback
24 import lsst.sphgeom
25 
26 import lsst.geom as geom
27 import lsst.pex.config as pexConfig
28 import lsst.pipe.base as pipeBase
29 from lsst.skymap import DiscreteSkyMap, BaseSkyMap
30 from lsst.pipe.base import ArgumentParser
31 from lsst.utils.timer import timeMethod
32 
33 
34 class MakeDiscreteSkyMapConfig(pexConfig.Config):
35  """Config for MakeDiscreteSkyMapTask
36  """
37  coaddName = pexConfig.Field(
38  doc="coadd name, e.g. deep, goodSeeing, chiSquared",
39  dtype=str,
40  default="deep",
41  )
42  skyMap = pexConfig.ConfigField(
43  dtype=BaseSkyMap.ConfigClass,
44  doc="SkyMap configuration parameters, excluding position and radius"
45  )
46  borderSize = pexConfig.Field(
47  doc="additional border added to the bounding box of the calexps, in degrees",
48  dtype=float,
49  default=0.0
50  )
51  doAppend = pexConfig.Field(
52  doc="append another tract to an existing DiscreteSkyMap on disk, if present?",
53  dtype=bool,
54  default=False
55  )
56  doWrite = pexConfig.Field(
57  doc="persist the skyMap?",
58  dtype=bool,
59  default=True,
60  )
61 
62  def setDefaults(self):
63  self.skyMapskyMap.tractOverlap = 0.0
64 
65 
66 class MakeDiscreteSkyMapRunner(pipeBase.TaskRunner):
67  """Run a task with all dataRefs at once, rather than one dataRef at a time.
68 
69  Call the run method of the task using two positional arguments:
70  - butler: data butler
71  - dataRefList: list of all dataRefs,
72  """
73  @staticmethod
74  def getTargetList(parsedCmd):
75  return [(parsedCmd.butler, parsedCmd.id.refList)]
76 
77  def __call__(self, args):
78  """
79  @param args Arguments for Task.run()
80 
81  @return:
82  - None if self.doReturnResults false
83  - A pipe_base Struct containing these fields if self.doReturnResults true:
84  - dataRef: the provided data reference
85  - metadata: task metadata after execution of run
86  - result: result returned by task run, or None if the task fails
87  """
88  butler, dataRefList = args
89  task = self.TaskClass(config=self.config, log=self.log)
90  result = None # in case the task fails
91  exitStatus = 0 # exit status for shell
92  if self.doRaise:
93  result = task.runDataRef(butler, dataRefList)
94  else:
95  try:
96  result = task.runDataRef(butler, dataRefList)
97  except Exception as e:
98  task.log.fatal("Failed: %s", e)
99  exitStatus = 1
100  if not isinstance(e, pipeBase.TaskError):
101  traceback.print_exc(file=sys.stderr)
102  for dataRef in dataRefList:
103  task.writeMetadata(dataRef)
104 
105  if self.doReturnResults:
106  return pipeBase.Struct(
107  dataRefList=dataRefList,
108  metadata=task.metadata,
109  result=result,
110  exitStatus=exitStatus,
111  )
112  else:
113  return pipeBase.Struct(
114  exitStatus=exitStatus,
115  )
116 
117 
118 class MakeDiscreteSkyMapTask(pipeBase.CmdLineTask):
119  """!Make a DiscreteSkyMap in a repository, using the bounding box of a set of calexps.
120 
121  The command-line and run signatures and config are sufficiently different from MakeSkyMapTask
122  that we don't inherit from it, but it is a replacement, so we use the same config/metadata names.
123  """
124  ConfigClass = MakeDiscreteSkyMapConfig
125  _DefaultName = "makeDiscreteSkyMap"
126  RunnerClass = MakeDiscreteSkyMapRunner
127 
128  def __init__(self, **kwargs):
129  super().__init__(**kwargs)
130 
131  def runDataRef(self, butler, dataRefList):
132  """Make a skymap from the bounds of the given set of calexps using the butler.
133 
134  Parameters
135  ----------
136  butler : `lsst.daf.persistence.Butler`
137  Gen2 data butler used to save the SkyMap
138  dataRefList : iterable
139  A list of Gen2 data refs of calexps used to determin the size and pointing of the SkyMap
140  Returns
141  -------
142  struct : `lsst.pipe.base.Struct`
143  The returned struct has one attribute, ``skyMap``, which holds the returned SkyMap
144  """
145  wcs_bbox_tuple_list = []
146  oldSkyMap = None
147  datasetName = self.config.coaddName + "Coadd_skyMap"
148  for dataRef in dataRefList:
149  if not dataRef.datasetExists("calexp"):
150  self.log.warning("CalExp for %s does not exist: ignoring", dataRef.dataId)
151  continue
152  wcs_bbox_tuple_list.append((dataRef.get("calexp_wcs", immediate=True),
153  dataRef.get("calexp_bbox", immediate=True)))
154  if self.config.doAppend and butler.datasetExists(datasetName):
155  oldSkyMap = butler.get(datasetName, immediate=True)
156  if not isinstance(oldSkyMap.config, DiscreteSkyMap.ConfigClass):
157  raise TypeError("Cannot append to existing non-discrete skymap")
158  compareLog = []
159  if not self.config.skyMap.compare(oldSkyMap.config, output=compareLog.append):
160  raise ValueError("Cannot append to existing skymap - configurations differ:", *compareLog)
161  result = self.runrun(wcs_bbox_tuple_list, oldSkyMap)
162  if self.config.doWrite:
163  butler.put(result.skyMap, datasetName)
164  return result
165 
166  @timeMethod
167  def run(self, wcs_bbox_tuple_list, oldSkyMap=None):
168  """Make a SkyMap from the bounds of the given set of calexp metadata.
169 
170  Parameters
171  ----------
172  wcs_bbox_tuple_list : iterable
173  A list of tuples with each element expected to be a (Wcs, Box2I) pair
174  oldSkyMap : `lsst.skymap.DiscreteSkyMap`, option
175  The SkyMap to extend if appending
176  Returns
177  -------
178  struct : `lsst.pipe.base.Struct
179  The returned struct has one attribute, ``skyMap``, which holds the returned SkyMap
180  """
181  self.log.info("Extracting bounding boxes of %d images", len(wcs_bbox_tuple_list))
182  points = []
183  for wcs, boxI in wcs_bbox_tuple_list:
184  boxD = geom.Box2D(boxI)
185  points.extend(wcs.pixelToSky(corner).getVector() for corner in boxD.getCorners())
186  if len(points) == 0:
187  raise RuntimeError("No data found from which to compute convex hull")
188  self.log.info("Computing spherical convex hull")
189  polygon = lsst.sphgeom.ConvexPolygon.convexHull(points)
190  if polygon is None:
191  raise RuntimeError(
192  "Failed to compute convex hull of the vertices of all calexp bounding boxes; "
193  "they may not be hemispherical."
194  )
195  circle = polygon.getBoundingCircle()
196 
197  skyMapConfig = DiscreteSkyMap.ConfigClass()
198  if oldSkyMap:
199  skyMapConfig.raList.extend(oldSkyMap.config.raList)
200  skyMapConfig.decList.extend(oldSkyMap.config.decList)
201  skyMapConfig.radiusList.extend(oldSkyMap.config.radiusList)
202  configIntersection = {k: getattr(self.config.skyMap, k)
203  for k in self.config.skyMap.toDict()
204  if k in skyMapConfig}
205  skyMapConfig.update(**configIntersection)
206  circleCenter = lsst.sphgeom.LonLat(circle.getCenter())
207  skyMapConfig.raList.append(circleCenter[0].asDegrees())
208  skyMapConfig.decList.append(circleCenter[1].asDegrees())
209  circleRadiusDeg = circle.getOpeningAngle().asDegrees()
210  skyMapConfig.radiusList.append(circleRadiusDeg + self.config.borderSize)
211  skyMap = DiscreteSkyMap(skyMapConfig)
212 
213  for tractInfo in skyMap:
214  wcs = tractInfo.getWcs()
215  posBox = geom.Box2D(tractInfo.getBBox())
216  pixelPosList = (
217  posBox.getMin(),
218  geom.Point2D(posBox.getMaxX(), posBox.getMinY()),
219  posBox.getMax(),
220  geom.Point2D(posBox.getMinX(), posBox.getMaxY()),
221  )
222  skyPosList = [wcs.pixelToSky(pos).getPosition(geom.degrees) for pos in pixelPosList]
223  posStrList = ["(%0.3f, %0.3f)" % tuple(skyPos) for skyPos in skyPosList]
224  self.log.info("tract %s has corners %s (RA, Dec deg) and %s x %s patches",
225  tractInfo.getId(), ", ".join(posStrList),
226  tractInfo.getNumPatches()[0], tractInfo.getNumPatches()[1])
227  return pipeBase.Struct(
228  skyMap=skyMap
229  )
230 
231  def _getConfigName(self):
232  """Return None to disable saving config
233 
234  There's only one SkyMap per repository, so the config is redundant, and checking it means we can't
235  easily overwrite or append to an existing repository.
236  """
237  return None
238 
239  def _getMetadataName(self):
240  """Return None to disable saving metadata
241 
242  The metadata is not interesting, and by not saving it we can eliminate a dataset type.
243  """
244  return None
245 
246  @classmethod
247  def _makeArgumentParser(cls):
248  parser = ArgumentParser(name=cls._DefaultName_DefaultName)
249  parser.add_id_argument(name="--id", datasetType="calexp", help="data ID, e.g. --id visit=123 ccd=1,2")
250  return parser
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
Make a DiscreteSkyMap in a repository, using the bounding box of a set of calexps.
def run(self, wcs_bbox_tuple_list, oldSkyMap=None)
static ConvexPolygon convexHull(std::vector< UnitVector3d > const &points)
convexHull returns the convex hull of the given set of points if it exists and throws an exception ot...
Definition: ConvexPolygon.h:65
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition: LonLat.h:48