LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
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 struct
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 
37 
38 class BaseSkyMapConfig(pexConfig.Config):
39  patchInnerDimensions = pexConfig.ListField(
40  doc="dimensions of inner region of patches (x,y pixels)",
41  dtype=int,
42  length=2,
43  default=(4000, 4000),
44  )
45  patchBorder = pexConfig.Field(
46  doc="border between patch inner and outer bbox (pixels)",
47  dtype=int,
48  default=100,
49  )
50  tractOverlap = pexConfig.Field(
51  doc="minimum overlap between adjacent sky tracts, on the sky (deg)",
52  dtype=float,
53  default=1.0,
54  )
55  pixelScale = pexConfig.Field(
56  doc="nominal pixel scale (arcsec/pixel)",
57  dtype=float,
58  default=0.333
59  )
60  projection = pexConfig.Field(
61  doc="one of the FITS WCS projection codes, such as:"
62  "- STG: stereographic projection"
63  "- MOL: Molleweide's projection"
64  "- TAN: tangent-plane projection",
65  dtype=str,
66  default="STG",
67  )
68  rotation = pexConfig.Field(
69  doc="Rotation for WCS (deg)",
70  dtype=float,
71  default=0,
72  )
73 
74 
75 class BaseSkyMap:
76  """A collection of overlapping Tracts that map part or all of the sky.
77 
78  See TractInfo for more information.
79 
80  Parameters
81  ----------
82  config : `BaseSkyMapConfig` or None (optional)
83  The configuration for this SkyMap; if None use the default config.
84 
85  Notes
86  -----
87  BaseSkyMap is an abstract base class. Subclasses must do the following:
88  define ``__init__`` and have it construct the TractInfo objects and put
89  them in ``__tractInfoList__`` define ``__getstate__`` and ``__setstate__``
90  to allow pickling (the butler saves sky maps using pickle);
91  see DodecaSkyMap for an example of how to do this. (Most of that code could
92  be moved into this base class, but that would make it harder to handle
93  older versions of pickle data.) define updateSha1 to add any
94  subclass-specific state to the hash.
95 
96  All SkyMap subclasses must be conceptually immutable; they must always
97  refer to the same set of mathematical tracts and patches even if the in-
98  memory representation of those objects changes.
99  """
100 
101  ConfigClass = BaseSkyMapConfig
102 
103  def __init__(self, config=None):
104  if config is None:
105  config = self.ConfigClassConfigClass()
106  config.freeze() # just to be sure, e.g. for pickling
107  self.configconfig = config
108  self._tractInfoList_tractInfoList = []
109  self._wcsFactory_wcsFactory = detail.WcsFactory(
110  pixelScale=Angle(self.configconfig.pixelScale, arcseconds),
111  projection=self.configconfig.projection,
112  rotation=Angle(self.configconfig.rotation, degrees),
113  )
114  self._sha1_sha1 = None
115 
116  def findTract(self, coord):
117  """Find the tract whose center is nearest the specified coord.
118 
119  Parameters
120  ----------
121  coord : `lsst.geom.SpherePoint`
122  ICRS sky coordinate to search for.
123 
124  Returns
125  -------
126  result : `TractInfo`
127  TractInfo of tract whose center is nearest the specified coord.
128 
129  Notes
130  -----
131  - If coord is equidistant between multiple sky tract centers then one
132  is arbitrarily chosen.
133 
134  - The default implementation is not very efficient; subclasses may wish
135  to override.
136 
137  **Warning:**
138  If tracts do not cover the whole sky then the returned tract may not
139  include the coord.
140  """
141  distTractInfoList = []
142  for i, tractInfo in enumerate(self):
143  angSep = coord.separation(tractInfo.getCtrCoord()).asDegrees()
144  # include index in order to disambiguate identical angSep values
145  distTractInfoList.append((angSep, i, tractInfo))
146  distTractInfoList.sort()
147  return distTractInfoList[0][2]
148 
149  def findTractPatchList(self, coordList):
150  """Find tracts and patches that overlap a region.
151 
152  Parameters
153  ----------
154  coordList : `list` of `lsst.geom.SpherePoint`
155  List of ICRS sky coordinates to search for.
156 
157  Returns
158  -------
159  reList : `list` of (`TractInfo`, `list` of `PatchInfo`)
160  For tracts and patches that contain, or may contain, the specified
161  region. The list will be empty if there is no overlap.
162 
163  Notes
164  -----
165  **warning:**
166  This uses a naive algorithm that may find some tracts and patches
167  that do not overlap the region (especially if the region is not a
168  rectangle aligned along patch x, y).
169  """
170  retList = []
171  for tractInfo in self:
172  patchList = tractInfo.findPatchList(coordList)
173  if patchList:
174  retList.append((tractInfo, patchList))
175  return retList
176 
177  def findClosestTractPatchList(self, coordList):
178  """Find closest tract and patches that overlap coordinates.
179 
180  Parameters
181  ----------
182  coordList : `lsst.geom.SpherePoint`
183  List of ICRS sky coordinates to search for.
184 
185  Returns
186  -------
187  retList : `list`
188  list of (TractInfo, list of PatchInfo) for tracts and patches
189  that contain, or may contain, the specified region.
190  The list will be empty if there is no overlap.
191  """
192  retList = []
193  for coord in coordList:
194  tractInfo = self.findTractfindTract(coord)
195  patchList = tractInfo.findPatchList(coordList)
196  if patchList and not (tractInfo, patchList) in retList:
197  retList.append((tractInfo, patchList))
198  return retList
199 
200  def __getitem__(self, ind):
201  return self._tractInfoList_tractInfoList[ind]
202 
203  def __iter__(self):
204  return iter(self._tractInfoList_tractInfoList)
205 
206  def __len__(self):
207  return len(self._tractInfoList_tractInfoList)
208 
209  def __hash__(self):
210  return hash(self.getSha1getSha1())
211 
212  def __eq__(self, other):
213  try:
214  return self.getSha1getSha1() == other.getSha1()
215  except AttributeError:
216  return NotImplemented
217 
218  def __ne__(self, other):
219  return not (self == other)
220 
221  def logSkyMapInfo(self, log):
222  """Write information about a sky map to supplied log
223 
224  Parameters
225  ----------
226  log : `lsst.log.Log`
227  Log object that information about skymap will be written
228  """
229  log.info("sky map has %s tracts" % (len(self),))
230  for tractInfo in self:
231  wcs = tractInfo.getWcs()
232  posBox = geom.Box2D(tractInfo.getBBox())
233  pixelPosList = (
234  posBox.getMin(),
235  geom.Point2D(posBox.getMaxX(), posBox.getMinY()),
236  posBox.getMax(),
237  geom.Point2D(posBox.getMinX(), posBox.getMaxY()),
238  )
239  skyPosList = [wcs.pixelToSky(pos).getPosition(geom.degrees) for pos in pixelPosList]
240  posStrList = ["(%0.3f, %0.3f)" % tuple(skyPos) for skyPos in skyPosList]
241  log.info("tract %s has corners %s (RA, Dec deg) and %s x %s patches" %
242  (tractInfo.getId(), ", ".join(posStrList),
243  tractInfo.getNumPatches()[0], tractInfo.getNumPatches()[1]))
244 
245  def getSha1(self):
246  """Return a SHA1 hash that uniquely identifies this SkyMap instance.
247 
248  Returns
249  -------
250  sha1 : `bytes`
251  A 20-byte hash that uniquely identifies this SkyMap instance.
252 
253  Notes
254  -----
255  Subclasses should almost always override ``updateSha1`` instead of
256  this function to add subclass-specific state to the hash.
257  """
258  if self._sha1_sha1 is None:
259  sha1 = hashlib.sha1()
260  sha1.update(type(self).__name__.encode('utf-8'))
261  configPacked = struct.pack(
262  "<iiidd3sd",
263  self.configconfig.patchInnerDimensions[0],
264  self.configconfig.patchInnerDimensions[1],
265  self.configconfig.patchBorder,
266  self.configconfig.tractOverlap,
267  self.configconfig.pixelScale,
268  self.configconfig.projection.encode('ascii'),
269  self.configconfig.rotation
270  )
271  sha1.update(configPacked)
272  self.updateSha1updateSha1(sha1)
273  self._sha1_sha1 = sha1.digest()
274  return self._sha1_sha1
275 
276  def updateSha1(self, sha1):
277  """Add subclass-specific state or configuration options to the SHA1.
278 
279  Parameters
280  ----------
281  sha1 : `hashlib.sha1`
282  A hashlib object on which `update` can be called to add
283  additional state to the hash.
284 
285  Notes
286  -----
287  This method is conceptually "protected" : it should be reimplemented by
288  all subclasses, but called only by the base class implementation of
289  `getSha1` .
290  """
291  raise NotImplementedError()
292 
293  SKYMAP_RUN_COLLECTION_NAME = "skymaps"
294 
295  SKYMAP_DATASET_TYPE_NAME = "skyMap"
296 
297  def register(self, name, butler):
298  """Add skymap, tract, and patch Dimension entries to the given Gen3
299  Butler.
300 
301  Parameters
302  ----------
303  name : `str`
304  The name of the skymap.
305  butler : `lsst.daf.butler.Butler`
306  The butler to add to.
307 
308  Raises
309  ------
310  lsst.daf.butler.registry.ConflictingDefinitionError
311  Raised if a different skymap exists with the same name.
312 
313  Notes
314  -----
315  Registering the same skymap multiple times (with the exact same
316  definition) is safe, but inefficient; most of the work of computing
317  the rows to be inserted must be done first in order to check for
318  consistency between the new skymap and any existing one.
319 
320  Re-registering a skymap with different tract and/or patch definitions
321  but the same summary information may not be detected as a conflict but
322  will never result in updating the skymap; there is intentionally no
323  way to modify a registered skymap (aside from manual administrative
324  operations on the database), as it is hard to guarantee that this can
325  be done without affecting reproducibility.
326  """
327  nxMax = 0
328  nyMax = 0
329  tractRecords = []
330  patchRecords = []
331  for tractInfo in self:
332  nx, ny = tractInfo.getNumPatches()
333  nxMax = max(nxMax, nx)
334  nyMax = max(nyMax, ny)
335  region = tractInfo.getOuterSkyPolygon()
336  centroid = SpherePoint(region.getCentroid())
337  tractRecords.append({
338  "skymap": name,
339  "tract": tractInfo.getId(),
340  "region": region,
341  "ra": centroid.getRa().asDegrees(),
342  "dec": centroid.getDec().asDegrees(),
343  })
344  for patchInfo in tractInfo:
345  cellX, cellY = patchInfo.getIndex()
346  patchRecords.append({
347  "skymap": name,
348  "tract": tractInfo.getId(),
349  "patch": tractInfo.getSequentialPatchIndex(patchInfo),
350  "cell_x": cellX,
351  "cell_y": cellY,
352  "region": patchInfo.getOuterSkyPolygon(tractInfo.getWcs()),
353  })
354  skyMapRecord = {
355  "skymap": name,
356  "hash": self.getSha1getSha1(),
357  "tract_max": len(self),
358  "patch_nx_max": nxMax,
359  "patch_ny_max": nyMax,
360  }
361  butler.registry.registerRun(self.SKYMAP_RUN_COLLECTION_NAMESKYMAP_RUN_COLLECTION_NAME)
362  # Kind of crazy that we've got three different capitalizations of
363  # "skymap" here, but that's what the various conventions (or at least
364  # precedents) dictate.
365  from lsst.daf.butler import DatasetType
366  from lsst.daf.butler.registry import ConflictingDefinitionError
367  datasetType = DatasetType(
368  name=self.SKYMAP_DATASET_TYPE_NAMESKYMAP_DATASET_TYPE_NAME,
369  dimensions=["skymap"],
370  storageClass="SkyMap",
371  universe=butler.registry.dimensions
372  )
373  butler.registry.registerDatasetType(datasetType)
374  with butler.transaction():
375  try:
376  inserted = butler.registry.syncDimensionData("skymap", skyMapRecord)
377  except ConflictingDefinitionError as err:
378  raise ConflictingDefinitionError(
379  f"SkyMap with hash {self.getSha1().hex()} is already registered with a different name."
380  ) from err
381  if inserted:
382  butler.registry.insertDimensionData("tract", *tractRecords)
383  butler.registry.insertDimensionData("patch", *patchRecords)
384  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:297
def findClosestTractPatchList(self, coordList)
Definition: baseSkyMap.py:177
def findTractPatchList(self, coordList)
Definition: baseSkyMap.py:149
def __init__(self, config=None)
Definition: baseSkyMap.py:103