LSST Applications g1cfbe01172+01aa18f939,g20cdd03214+31e6b93548,g28da252d5a+ea8665a95b,g2bbee38e9b+9ec6cc348d,g2bc492864f+9ec6cc348d,g347aa1857d+9ec6cc348d,g3a166c0a6a+9ec6cc348d,g4322eb9e3a+65eff1e020,g461a3dce89+b86e4b8053,g50ff169b8f+f991eae79d,g52b1c1532d+b86e4b8053,g607f77f49a+31e6b93548,g78056777b3+8ae2798781,g858d7b2824+31e6b93548,g8cd86fa7b1+4851e61ca4,g9ddcbc5298+f24b38b85a,ga1e77700b3+3309dba821,gae0086650b+b86e4b8053,gb0e22166c9+6076c0b52b,gbb886bcc26+dccb771098,gbd462c55f0+dc07f8e65d,gc0c51c7ec2+31e6b93548,gc120e1dc64+a417ce3171,gc28159a63d+9ec6cc348d,gc2a6998b7e+f95f64aeae,gcdd4ae20e8+507450c4cd,gcf0d15dbbd+507450c4cd,gd1535ee943+bcf88ba65f,gd598c5cd71+66126f91fb,gdaeeff99f8+006e14e809,gdbce86181e+39d5515b1a,ge3d4d395c2+b12d4d6a95,ge79ae78c31+9ec6cc348d,gf048a9a2f4+d9c36e6b63,gfbcc870c63+ea41c4420b,w.2024.27
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"""
23todo: Consider tweaking pixel scale so the average scale is as specified,
24rather than the scale at the center.
25"""
26
27__all__ = ["BaseSkyMapConfig", "BaseSkyMap"]
28
29import hashlib
30import numpy as np
31
32import lsst.geom as geom
33import lsst.pex.config as pexConfig
34from lsst.geom import Angle, arcseconds, degrees
35from . import detail
36from .tractBuilder import tractBuilderRegistry
37
38
39class 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
77 ``config.tractBuilder["legacy"].patchInnerDimensions``.
78
79 Returns
80 -------
81 innerDimensions : `list` [`int`, `int`]
82 """
83 return self.tractBuilder["legacy"].patchInnerDimensions
84
85 def setPatchInnerDimensions(self, value):
86 """Set the patch inner dimensions, for backwards compatibility.
87
88 This value is only used with the ``legacy`` tract builder,
89 and is ignored otherwise. In general, the config should be
90 accessed directly with
91 ``config.tractBuilder["legacy"].patchInnerDimensions``.
92
93 Parameters
94 ----------
95 value : `list` [`int`, `int`]
96 """
97 self.tractBuilder["legacy"].patchInnerDimensions = value
98
99 patchInnerDimensions = property(getPatchInnerDimensions, setPatchInnerDimensions)
100
101 def getPatchBorder(self):
102 """Get the patch border, for backwards compatibility.
103
104 This value is only used with the ``legacy`` tract builder,
105 and is ignored otherwise. In general, the config should be
106 accessed directly with config.tractBuilder["legacy"].patchBorder.
107
108 Returns
109 -------
110 border: `int`
111 """
112 return self.tractBuilder["legacy"].patchBorder
113
114 def setPatchBorder(self, value):
115 """Set the patch border, for backwards compatibility.
116
117 This value is only used with the ``legacy`` tract builder,
118 and is ignored otherwise. In general, the config should be
119 accessed directly with config.tractBuilder["legacy"].patchBorder.
120
121 Parameters
122 ----------
123 border: `int`
124 """
125 self.tractBuilder["legacy"].patchBorder = value
126
127 patchBorder = property(getPatchBorder, setPatchBorder)
128
129
131 """A collection of overlapping Tracts that map part or all of the sky.
132
133 See TractInfo for more information.
134
135 Parameters
136 ----------
137 config : `BaseSkyMapConfig` or `None` (optional)
138 The configuration for this SkyMap; if None use the default config.
139
140 Notes
141 -----
142 BaseSkyMap is an abstract base class. Subclasses must do the following:
143 define ``__init__`` and have it construct the TractInfo objects and put
144 them in ``__tractInfoList__`` define ``__getstate__`` and ``__setstate__``
145 to allow pickling (the butler saves sky maps using pickle);
146 see DodecaSkyMap for an example of how to do this. (Most of that code could
147 be moved into this base class, but that would make it harder to handle
148 older versions of pickle data.) define updateSha1 to add any
149 subclass-specific state to the hash.
150
151 All SkyMap subclasses must be conceptually immutable; they must always
152 refer to the same set of mathematical tracts and patches even if the in-
153 memory representation of those objects changes.
154 """
155
156 ConfigClass = BaseSkyMapConfig
157
158 def __init__(self, config=None):
159 if config is None:
160 config = self.ConfigClass()
161 config.freeze() # just to be sure, e.g. for pickling
162 self.config = config
165 pixelScale=Angle(self.config.pixelScale, arcseconds),
166 projection=self.config.projection,
167 rotation=Angle(self.config.rotation, degrees),
168 )
169 self._sha1 = None
170 self._tractBuilder = config.tractBuilder.apply()
171
172 def findTract(self, coord):
173 """Find the tract whose center is nearest the specified coord.
174
175 Parameters
176 ----------
177 coord : `lsst.geom.SpherePoint`
178 ICRS sky coordinate to search for.
179
180 Returns
181 -------
182 result : `TractInfo`
183 TractInfo of tract whose center is nearest the specified coord.
184
185 Notes
186 -----
187 - If coord is equidistant between multiple sky tract centers then one
188 is arbitrarily chosen.
189
190 - The default implementation is not very efficient; subclasses may wish
191 to override.
192
193 .. warning::
194
195 If tracts do not cover the whole sky then the returned tract may not
196 include the coord.
197 """
198 distTractInfoList = []
199 for i, tractInfo in enumerate(self):
200 angSep = coord.separation(tractInfo.getCtrCoord()).asDegrees()
201 # include index in order to disambiguate identical angSep values
202 distTractInfoList.append((angSep, i, tractInfo))
203 distTractInfoList.sort()
204 return distTractInfoList[0][2]
205
206 def findTractIdArray(self, ra, dec, degrees=False):
207 """Find array of tract IDs with vectorized operations (where
208 supported).
209
210 If a given sky map does not support vectorized operations, then a loop
211 over findTract will be called.
212
213 Parameters
214 ----------
215 ra : `numpy.ndarray`
216 Array of Right Ascension. Units are radians unless
217 degrees=True.
218 dec : `numpy.ndarray`
219 Array of Declination. Units are radians unless
220 degrees=True.
221 degrees : `bool`, optional
222 Input ra, dec arrays are degrees if `True`.
223
224 Returns
225 -------
226 tractId : `numpy.ndarray`
227 Array of tract IDs
228
229 Notes
230 -----
231 - If coord is equidistant between multiple sky tract centers then one
232 is arbitrarily chosen.
233
234 .. warning::
235
236 If tracts do not cover the whole sky then the returned tract may not
237 include the given ra/dec.
238 """
239 units = geom.degrees if degrees else geom.radians
240 coords = [geom.SpherePoint(r*units, d*units) for r, d in zip(np.atleast_1d(ra),
241 np.atleast_1d(dec))]
242
243 tractId = np.array([self.findTract(coord).getId() for coord in coords])
244
245 return tractId
246
247 def findTractPatchList(self, coordList):
248 """Find tracts and patches that overlap a region.
249
250 Parameters
251 ----------
252 coordList : `list` of `lsst.geom.SpherePoint`
253 List of ICRS sky coordinates to search for.
254
255 Returns
256 -------
257 reList : `list` of (`TractInfo`, `list` of `PatchInfo`)
258 For tracts and patches that contain, or may contain, the specified
259 region. The list will be empty if there is no overlap.
260
261 Notes
262 -----
263 .. warning::
264
265 This uses a naive algorithm that may find some tracts and patches
266 that do not overlap the region (especially if the region is not a
267 rectangle aligned along patch x, y).
268 """
269 retList = []
270 for tractInfo in self:
271 patchList = tractInfo.findPatchList(coordList)
272 if patchList:
273 retList.append((tractInfo, patchList))
274 return retList
275
276 def findClosestTractPatchList(self, coordList):
277 """Find closest tract and patches that overlap coordinates.
278
279 Parameters
280 ----------
281 coordList : `lsst.geom.SpherePoint`
282 List of ICRS sky coordinates to search for.
283
284 Returns
285 -------
286 retList : `list`
287 list of (TractInfo, list of PatchInfo) for tracts and patches
288 that contain, or may contain, the specified region.
289 The list will be empty if there is no overlap.
290 """
291 retList = []
292 for coord in coordList:
293 tractInfo = self.findTract(coord)
294 patchList = tractInfo.findPatchList(coordList)
295 if patchList and not (tractInfo, patchList) in retList:
296 retList.append((tractInfo, patchList))
297 return retList
298
299 def __getitem__(self, ind):
300 return self._tractInfoList[ind]
301
302 def __iter__(self):
303 return iter(self._tractInfoList)
304
305 def __len__(self):
306 return len(self._tractInfoList)
307
308 def __hash__(self):
309 return hash(self.getSha1())
310
311 def __eq__(self, other):
312 try:
313 return self.getSha1() == other.getSha1()
314 except AttributeError:
315 return NotImplemented
316
317 def __ne__(self, other):
318 return not (self == other)
319
320 def logSkyMapInfo(self, log):
321 """Write information about a sky map to supplied log
322
323 Parameters
324 ----------
325 log : `logging.Logger`
326 Log object that information about skymap will be written.
327 """
328 log.info("sky map has %s tracts" % (len(self),))
329 for tractInfo in self:
330 wcs = tractInfo.getWcs()
331 posBox = geom.Box2D(tractInfo.getBBox())
332 pixelPosList = (
333 posBox.getMin(),
334 geom.Point2D(posBox.getMaxX(), posBox.getMinY()),
335 posBox.getMax(),
336 geom.Point2D(posBox.getMinX(), posBox.getMaxY()),
337 )
338 skyPosList = [wcs.pixelToSky(pos).getPosition(geom.degrees) for pos in pixelPosList]
339 posStrList = ["(%0.3f, %0.3f)" % tuple(skyPos) for skyPos in skyPosList]
340 log.info("tract %s has corners %s (RA, Dec deg) and %s x %s patches" %
341 (tractInfo.getId(), ", ".join(posStrList),
342 tractInfo.getNumPatches()[0], tractInfo.getNumPatches()[1]))
343
344 def getSha1(self):
345 """Return a SHA1 hash that uniquely identifies this SkyMap instance.
346
347 Returns
348 -------
349 sha1 : `bytes`
350 A 20-byte hash that uniquely identifies this SkyMap instance.
351
352 Notes
353 -----
354 Subclasses should almost always override ``updateSha1`` instead of
355 this function to add subclass-specific state to the hash.
356 """
357 if self._sha1 is None:
358 sha1 = hashlib.sha1()
359 sha1.update(type(self).__name__.encode('utf-8'))
360 configPacked = self._tractBuilder.getPackedConfig(self.config)
361 sha1.update(configPacked)
362 self.updateSha1(sha1)
363 self._sha1 = sha1.digest()
364 return self._sha1
365
366 def updateSha1(self, sha1):
367 """Add subclass-specific state or configuration options to the SHA1.
368
369 Parameters
370 ----------
371 sha1 : `hashlib.sha1`
372 A hashlib object on which `update` can be called to add
373 additional state to the hash.
374
375 Notes
376 -----
377 This method is conceptually "protected" : it should be reimplemented by
378 all subclasses, but called only by the base class implementation of
379 `getSha1` .
380 """
381 raise NotImplementedError()
382
383 SKYMAP_RUN_COLLECTION_NAME = "skymaps"
384
385 SKYMAP_DATASET_TYPE_NAME = "skyMap"
386
387 def register(self, name, butler):
388 """Add skymap, tract, and patch Dimension entries to the given Gen3
389 Butler.
390
391 Parameters
392 ----------
393 name : `str`
394 The name of the skymap.
395 butler : `lsst.daf.butler.Butler`
396 The butler to add to.
397
398 Raises
399 ------
400 lsst.daf.butler.registry.ConflictingDefinitionError
401 Raised if a different skymap exists with the same name.
402
403 Notes
404 -----
405 Registering the same skymap multiple times (with the exact same
406 definition) is safe, but inefficient; most of the work of computing
407 the rows to be inserted must be done first in order to check for
408 consistency between the new skymap and any existing one.
409
410 Re-registering a skymap with different tract and/or patch definitions
411 but the same summary information may not be detected as a conflict but
412 will never result in updating the skymap; there is intentionally no
413 way to modify a registered skymap (aside from manual administrative
414 operations on the database), as it is hard to guarantee that this can
415 be done without affecting reproducibility.
416 """
417 numPatches = [tractInfo.getNumPatches() for tractInfo in self]
418 nxMax = max(nn[0] for nn in numPatches)
419 nyMax = max(nn[1] for nn in numPatches)
420
421 skyMapRecord = {
422 "skymap": name,
423 "hash": self.getSha1(),
424 "tract_max": len(self),
425 "patch_nx_max": nxMax,
426 "patch_ny_max": nyMax,
427 }
429 # Kind of crazy that we've got three different capitalizations of
430 # "skymap" here, but that's what the various conventions (or at least
431 # precedents) dictate.
432 from lsst.daf.butler import DatasetType
433 from lsst.daf.butler.registry import ConflictingDefinitionError
434 datasetType = DatasetType(
435 name=self.SKYMAP_DATASET_TYPE_NAME,
436 dimensions=["skymap"],
437 storageClass="SkyMap",
438 universe=butler.dimensions
439 )
440 butler.registry.registerDatasetType(datasetType)
441 with butler.transaction():
442 try:
443 inserted = butler.registry.syncDimensionData("skymap", skyMapRecord)
444 except ConflictingDefinitionError as err:
445 raise ConflictingDefinitionError(
446 f"SkyMap with hash {self.getSha1().hex()} is already registered with a different name."
447 ) from err
448 if inserted:
449 for tractInfo in self:
450 tractId = tractInfo.getId()
451 tractRegion = tractInfo.getOuterSkyPolygon()
452 tractWcs = tractInfo.getWcs()
453 tractRecord = dict(
454 skymap=name,
455 tract=tractId,
456 region=tractRegion,
457 )
458 butler.registry.insertDimensionData("tract", tractRecord)
459
460 patchRecords = []
461 for patchInfo in tractInfo:
462 xx, yy = patchInfo.getIndex()
463 patchRecords.append(
464 dict(
465 skymap=name,
466 tract=tractId,
467 patch=tractInfo.getSequentialPatchIndex(patchInfo),
468 cell_x=xx,
469 cell_y=yy,
470 region=patchInfo.getOuterSkyPolygon(tractWcs),
471 )
472 )
473 butler.registry.insertDimensionData("patch", *patchRecords)
474
475 butler.put(self, datasetType, {"skymap": name}, run=self.SKYMAP_RUN_COLLECTION_NAMESKYMAP_RUN_COLLECTION_NAME)
int max
A class representing an angle.
Definition Angle.h:128
A floating-point coordinate rectangle geometry.
Definition Box.h:413
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
findTractIdArray(self, ra, dec, degrees=False)
findClosestTractPatchList(self, coordList)
findTractPatchList(self, coordList)