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
baseSkyMap.py
Go to the documentation of this file.
1 # LSST Data Management System
2 # Copyright 2008, 2009, 2010 LSST Corporation.
3 #
4 # This product includes software developed by the
5 # LSST Project (http://www.lsst.org/).
6 #
7 # This program is free software: you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation, either version 3 of the License, or
10 # (at your option) any later version.
11 #
12 # This program is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
16 #
17 # You should have received a copy of the LSST License Statement and
18 # the GNU General Public License along with this program. If not,
19 # see <http://www.lsstcorp.org/LegalNotices/>.
20 #
21 
22 """
23 todo: Consider tweaking pixel scale so the average scale is as specified,
24 rather than the scale at the center.
25 """
26 
27 __all__ = ["BaseSkyMapConfig", "BaseSkyMap"]
28 
29 import hashlib
30 import numpy as np
31 
32 import lsst.geom as geom
33 import lsst.pex.config as pexConfig
34 from lsst.geom import SpherePoint, Angle, arcseconds, degrees
35 from . import detail
36 from .tractBuilder import tractBuilderRegistry
37 
38 
39 class BaseSkyMapConfig(pexConfig.Config):
40  tractBuilder = tractBuilderRegistry.makeField(
41  doc="Algorithm for creating patches within the tract.",
42  default="legacy"
43  )
44 
45  tractOverlap = pexConfig.Field(
46  doc="minimum overlap between adjacent sky tracts, on the sky (deg)",
47  dtype=float,
48  default=1.0,
49  )
50  pixelScale = pexConfig.Field(
51  doc="nominal pixel scale (arcsec/pixel)",
52  dtype=float,
53  default=0.333
54  )
55  projection = pexConfig.Field(
56  doc="one of the FITS WCS projection codes, such as:"
57  "- STG: stereographic projection"
58  "- MOL: Molleweide's projection"
59  "- TAN: tangent-plane projection",
60  dtype=str,
61  default="STG",
62  )
63  rotation = pexConfig.Field(
64  doc="Rotation for WCS (deg)",
65  dtype=float,
66  default=0,
67  )
68 
69  # Backwards compatibility
70  # We can't use the @property decorator because it makes pexConfig sad.
72  """Get the patch inner dimensions, for backwards compatibility.
73 
74  This value is only used with the ``legacy`` tract builder,
75  and is ignored otherwise. In general, the config should be
76  accessed directly with config.tractBuilder["legacy"].patchInnerDimensions.
77 
78  Returns
79  -------
80  innerDimensions : `list` [`int`, `int`]
81  """
82  return self.tractBuildertractBuilder["legacy"].patchInnerDimensions
83 
84  def setPatchInnerDimensions(self, value):
85  """Set the patch inner dimensions, for backwards compatibility.
86 
87  This value is only used with the ``legacy`` tract builder,
88  and is ignored otherwise. In general, the config should be
89  accessed directly with config.tractBuilder["legacy"].patchInnerDimensions.
90 
91  Parameters
92  ----------
93  value : `list` [`int`, `int`]
94  """
95  self.tractBuildertractBuilder["legacy"].patchInnerDimensions = value
96 
97  patchInnerDimensions = property(getPatchInnerDimensions, setPatchInnerDimensions)
98 
99  def getPatchBorder(self):
100  """Get the patch border, for backwards compatibility.
101 
102  This value is only used with the ``legacy`` tract builder,
103  and is ignored otherwise. In general, the config should be
104  accessed directly with config.tractBuilder["legacy"].patchBorder.
105 
106  Returns
107  -------
108  border: `int`
109  """
110  return self.tractBuildertractBuilder["legacy"].patchBorder
111 
112  def setPatchBorder(self, value):
113  """Set the patch border, for backwards compatibility.
114 
115  This value is only used with the ``legacy`` tract builder,
116  and is ignored otherwise. In general, the config should be
117  accessed directly with config.tractBuilder["legacy"].patchBorder.
118 
119  Parameters
120  -------
121  border: `int`
122  """
123  self.tractBuildertractBuilder["legacy"].patchBorder = value
124 
125  patchBorder = property(getPatchBorder, setPatchBorder)
126 
127 
129  """A collection of overlapping Tracts that map part or all of the sky.
130 
131  See TractInfo for more information.
132 
133  Parameters
134  ----------
135  config : `BaseSkyMapConfig` or None (optional)
136  The configuration for this SkyMap; if None use the default config.
137 
138  Notes
139  -----
140  BaseSkyMap is an abstract base class. Subclasses must do the following:
141  define ``__init__`` and have it construct the TractInfo objects and put
142  them in ``__tractInfoList__`` define ``__getstate__`` and ``__setstate__``
143  to allow pickling (the butler saves sky maps using pickle);
144  see DodecaSkyMap for an example of how to do this. (Most of that code could
145  be moved into this base class, but that would make it harder to handle
146  older versions of pickle data.) define updateSha1 to add any
147  subclass-specific state to the hash.
148 
149  All SkyMap subclasses must be conceptually immutable; they must always
150  refer to the same set of mathematical tracts and patches even if the in-
151  memory representation of those objects changes.
152  """
153 
154  ConfigClass = BaseSkyMapConfig
155 
156  def __init__(self, config=None):
157  if config is None:
158  config = self.ConfigClassConfigClass()
159  config.freeze() # just to be sure, e.g. for pickling
160  self.configconfig = config
161  self._tractInfoList_tractInfoList = []
162  self._wcsFactory_wcsFactory = detail.WcsFactory(
163  pixelScale=Angle(self.configconfig.pixelScale, arcseconds),
164  projection=self.configconfig.projection,
165  rotation=Angle(self.configconfig.rotation, degrees),
166  )
167  self._sha1_sha1 = None
168  self._tractBuilder_tractBuilder = config.tractBuilder.apply()
169 
170  def findTract(self, coord):
171  """Find the tract whose center is nearest the specified coord.
172 
173  Parameters
174  ----------
175  coord : `lsst.geom.SpherePoint`
176  ICRS sky coordinate to search for.
177 
178  Returns
179  -------
180  result : `TractInfo`
181  TractInfo of tract whose center is nearest the specified coord.
182 
183  Notes
184  -----
185  - If coord is equidistant between multiple sky tract centers then one
186  is arbitrarily chosen.
187 
188  - The default implementation is not very efficient; subclasses may wish
189  to override.
190 
191  **Warning:**
192  If tracts do not cover the whole sky then the returned tract may not
193  include the coord.
194  """
195  distTractInfoList = []
196  for i, tractInfo in enumerate(self):
197  angSep = coord.separation(tractInfo.getCtrCoord()).asDegrees()
198  # include index in order to disambiguate identical angSep values
199  distTractInfoList.append((angSep, i, tractInfo))
200  distTractInfoList.sort()
201  return distTractInfoList[0][2]
202 
203  def findTractIdArray(self, ra, dec, degrees=False):
204  """Find array of tract IDs with vectorized operations (where supported).
205 
206  If a given sky map does not support vectorized operations, then a loop
207  over findTract will be called.
208 
209  Parameters
210  ----------
211  ra : `np.ndarray`
212  Array of Right Ascension. Units are radians unless
213  degrees=True.
214  dec : `np.ndarray`
215  Array of Declination. Units are radians unless
216  degrees=True.
217  degrees : `bool`, optional
218  Input ra, dec arrays are degrees if True.
219 
220  Returns
221  -------
222  tractId : `np.ndarray`
223  Array of tract IDs
224 
225  Notes
226  -----
227  - If coord is equidistant between multiple sky tract centers then one
228  is arbitrarily chosen.
229 
230  **Warning:**
231  If tracts do not cover the whole sky then the returned tract may not
232  include the given ra/dec.
233  """
234  units = geom.degrees if degrees else geom.radians
235  coords = [geom.SpherePoint(r*units, d*units) for r, d in zip(np.atleast_1d(ra),
236  np.atleast_1d(dec))]
237 
238  tractId = np.array([self.findTractfindTract(coord).getId() for coord in coords])
239 
240  return tractId
241 
242  def findTractPatchList(self, coordList):
243  """Find tracts and patches that overlap a region.
244 
245  Parameters
246  ----------
247  coordList : `list` of `lsst.geom.SpherePoint`
248  List of ICRS sky coordinates to search for.
249 
250  Returns
251  -------
252  reList : `list` of (`TractInfo`, `list` of `PatchInfo`)
253  For tracts and patches that contain, or may contain, the specified
254  region. The list will be empty if there is no overlap.
255 
256  Notes
257  -----
258  **warning:**
259  This uses a naive algorithm that may find some tracts and patches
260  that do not overlap the region (especially if the region is not a
261  rectangle aligned along patch x, y).
262  """
263  retList = []
264  for tractInfo in self:
265  patchList = tractInfo.findPatchList(coordList)
266  if patchList:
267  retList.append((tractInfo, patchList))
268  return retList
269 
270  def findClosestTractPatchList(self, coordList):
271  """Find closest tract and patches that overlap coordinates.
272 
273  Parameters
274  ----------
275  coordList : `lsst.geom.SpherePoint`
276  List of ICRS sky coordinates to search for.
277 
278  Returns
279  -------
280  retList : `list`
281  list of (TractInfo, list of PatchInfo) for tracts and patches
282  that contain, or may contain, the specified region.
283  The list will be empty if there is no overlap.
284  """
285  retList = []
286  for coord in coordList:
287  tractInfo = self.findTractfindTract(coord)
288  patchList = tractInfo.findPatchList(coordList)
289  if patchList and not (tractInfo, patchList) in retList:
290  retList.append((tractInfo, patchList))
291  return retList
292 
293  def __getitem__(self, ind):
294  return self._tractInfoList_tractInfoList[ind]
295 
296  def __iter__(self):
297  return iter(self._tractInfoList_tractInfoList)
298 
299  def __len__(self):
300  return len(self._tractInfoList_tractInfoList)
301 
302  def __hash__(self):
303  return hash(self.getSha1getSha1())
304 
305  def __eq__(self, other):
306  try:
307  return self.getSha1getSha1() == other.getSha1()
308  except AttributeError:
309  return NotImplemented
310 
311  def __ne__(self, other):
312  return not (self == other)
313 
314  def logSkyMapInfo(self, log):
315  """Write information about a sky map to supplied log
316 
317  Parameters
318  ----------
319  log : `lsst.log.Log`
320  Log object that information about skymap will be written
321  """
322  log.info("sky map has %s tracts" % (len(self),))
323  for tractInfo in self:
324  wcs = tractInfo.getWcs()
325  posBox = geom.Box2D(tractInfo.getBBox())
326  pixelPosList = (
327  posBox.getMin(),
328  geom.Point2D(posBox.getMaxX(), posBox.getMinY()),
329  posBox.getMax(),
330  geom.Point2D(posBox.getMinX(), posBox.getMaxY()),
331  )
332  skyPosList = [wcs.pixelToSky(pos).getPosition(geom.degrees) for pos in pixelPosList]
333  posStrList = ["(%0.3f, %0.3f)" % tuple(skyPos) for skyPos in skyPosList]
334  log.info("tract %s has corners %s (RA, Dec deg) and %s x %s patches" %
335  (tractInfo.getId(), ", ".join(posStrList),
336  tractInfo.getNumPatches()[0], tractInfo.getNumPatches()[1]))
337 
338  def getSha1(self):
339  """Return a SHA1 hash that uniquely identifies this SkyMap instance.
340 
341  Returns
342  -------
343  sha1 : `bytes`
344  A 20-byte hash that uniquely identifies this SkyMap instance.
345 
346  Notes
347  -----
348  Subclasses should almost always override ``updateSha1`` instead of
349  this function to add subclass-specific state to the hash.
350  """
351  if self._sha1_sha1 is None:
352  sha1 = hashlib.sha1()
353  sha1.update(type(self).__name__.encode('utf-8'))
354  configPacked = self._tractBuilder_tractBuilder.getPackedConfig(self.configconfig)
355  sha1.update(configPacked)
356  self.updateSha1updateSha1(sha1)
357  self._sha1_sha1 = sha1.digest()
358  return self._sha1_sha1
359 
360  def updateSha1(self, sha1):
361  """Add subclass-specific state or configuration options to the SHA1.
362 
363  Parameters
364  ----------
365  sha1 : `hashlib.sha1`
366  A hashlib object on which `update` can be called to add
367  additional state to the hash.
368 
369  Notes
370  -----
371  This method is conceptually "protected" : it should be reimplemented by
372  all subclasses, but called only by the base class implementation of
373  `getSha1` .
374  """
375  raise NotImplementedError()
376 
377  SKYMAP_RUN_COLLECTION_NAME = "skymaps"
378 
379  SKYMAP_DATASET_TYPE_NAME = "skyMap"
380 
381  def register(self, name, butler):
382  """Add skymap, tract, and patch Dimension entries to the given Gen3
383  Butler.
384 
385  Parameters
386  ----------
387  name : `str`
388  The name of the skymap.
389  butler : `lsst.daf.butler.Butler`
390  The butler to add to.
391 
392  Raises
393  ------
394  lsst.daf.butler.registry.ConflictingDefinitionError
395  Raised if a different skymap exists with the same name.
396 
397  Notes
398  -----
399  Registering the same skymap multiple times (with the exact same
400  definition) is safe, but inefficient; most of the work of computing
401  the rows to be inserted must be done first in order to check for
402  consistency between the new skymap and any existing one.
403 
404  Re-registering a skymap with different tract and/or patch definitions
405  but the same summary information may not be detected as a conflict but
406  will never result in updating the skymap; there is intentionally no
407  way to modify a registered skymap (aside from manual administrative
408  operations on the database), as it is hard to guarantee that this can
409  be done without affecting reproducibility.
410  """
411  nxMax = 0
412  nyMax = 0
413  tractRecords = []
414  patchRecords = []
415  for tractInfo in self:
416  nx, ny = tractInfo.getNumPatches()
417  nxMax = max(nxMax, nx)
418  nyMax = max(nyMax, ny)
419  region = tractInfo.getOuterSkyPolygon()
420  centroid = SpherePoint(region.getCentroid())
421  tractRecords.append({
422  "skymap": name,
423  "tract": tractInfo.getId(),
424  "region": region,
425  "ra": centroid.getRa().asDegrees(),
426  "dec": centroid.getDec().asDegrees(),
427  })
428  for patchInfo in tractInfo:
429  cellX, cellY = patchInfo.getIndex()
430  patchRecords.append({
431  "skymap": name,
432  "tract": tractInfo.getId(),
433  "patch": tractInfo.getSequentialPatchIndex(patchInfo),
434  "cell_x": cellX,
435  "cell_y": cellY,
436  "region": patchInfo.getOuterSkyPolygon(tractInfo.getWcs()),
437  })
438  skyMapRecord = {
439  "skymap": name,
440  "hash": self.getSha1getSha1(),
441  "tract_max": len(self),
442  "patch_nx_max": nxMax,
443  "patch_ny_max": nyMax,
444  }
445  butler.registry.registerRun(self.SKYMAP_RUN_COLLECTION_NAMESKYMAP_RUN_COLLECTION_NAME)
446  # Kind of crazy that we've got three different capitalizations of
447  # "skymap" here, but that's what the various conventions (or at least
448  # precedents) dictate.
449  from lsst.daf.butler import DatasetType
450  from lsst.daf.butler.registry import ConflictingDefinitionError
451  datasetType = DatasetType(
452  name=self.SKYMAP_DATASET_TYPE_NAMESKYMAP_DATASET_TYPE_NAME,
453  dimensions=["skymap"],
454  storageClass="SkyMap",
455  universe=butler.registry.dimensions
456  )
457  butler.registry.registerDatasetType(datasetType)
458  with butler.transaction():
459  try:
460  inserted = butler.registry.syncDimensionData("skymap", skyMapRecord)
461  except ConflictingDefinitionError as err:
462  raise ConflictingDefinitionError(
463  f"SkyMap with hash {self.getSha1().hex()} is already registered with a different name."
464  ) from err
465  if inserted:
466  butler.registry.insertDimensionData("tract", *tractRecords)
467  butler.registry.insertDimensionData("patch", *patchRecords)
468  butler.put(self, datasetType, {"skymap": name}, run=self.SKYMAP_RUN_COLLECTION_NAMESKYMAP_RUN_COLLECTION_NAME)
int max
table::Key< int > type
Definition: Detector.cc:163
A class representing an angle.
Definition: Angle.h:127
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
def register(self, name, butler)
Definition: baseSkyMap.py:381
def findClosestTractPatchList(self, coordList)
Definition: baseSkyMap.py:270
def findTractPatchList(self, coordList)
Definition: baseSkyMap.py:242
def findTractIdArray(self, ra, dec, degrees=False)
Definition: baseSkyMap.py:203
def __init__(self, config=None)
Definition: baseSkyMap.py:156