23 todo: Consider tweaking pixel scale so the average scale is as specified,
24 rather than the scale at the center.
27 __all__ = [
"BaseSkyMapConfig",
"BaseSkyMap"]
34 from lsst.geom import SpherePoint, Angle, arcseconds, degrees
36 from .tractBuilder
import tractBuilderRegistry
40 tractBuilder = tractBuilderRegistry.makeField(
41 doc=
"Algorithm for creating patches within the tract.",
45 tractOverlap = pexConfig.Field(
46 doc=
"minimum overlap between adjacent sky tracts, on the sky (deg)",
50 pixelScale = pexConfig.Field(
51 doc=
"nominal pixel scale (arcsec/pixel)",
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",
63 rotation = pexConfig.Field(
64 doc=
"Rotation for WCS (deg)",
72 """Get the patch inner dimensions, for backwards compatibility.
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.
80 innerDimensions : `list` [`int`, `int`]
82 return self.
tractBuildertractBuilder[
"legacy"].patchInnerDimensions
85 """Set the patch inner dimensions, for backwards compatibility.
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.
93 value : `list` [`int`, `int`]
95 self.
tractBuildertractBuilder[
"legacy"].patchInnerDimensions = value
97 patchInnerDimensions = property(getPatchInnerDimensions, setPatchInnerDimensions)
100 """Get the patch border, for backwards compatibility.
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.
110 return self.
tractBuildertractBuilder[
"legacy"].patchBorder
113 """Set the patch border, for backwards compatibility.
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.
123 self.
tractBuildertractBuilder[
"legacy"].patchBorder = value
125 patchBorder = property(getPatchBorder, setPatchBorder)
129 """A collection of overlapping Tracts that map part or all of the sky.
131 See TractInfo for more information.
135 config : `BaseSkyMapConfig` or None (optional)
136 The configuration for this SkyMap; if None use the default config.
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.
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.
154 ConfigClass = BaseSkyMapConfig
163 pixelScale=
Angle(self.
configconfig.pixelScale, arcseconds),
164 projection=self.
configconfig.projection,
165 rotation=
Angle(self.
configconfig.rotation, degrees),
167 self.
_sha1_sha1 =
None
168 self.
_tractBuilder_tractBuilder = config.tractBuilder.apply()
171 """Find the tract whose center is nearest the specified coord.
175 coord : `lsst.geom.SpherePoint`
176 ICRS sky coordinate to search for.
181 TractInfo of tract whose center is nearest the specified coord.
185 - If coord is equidistant between multiple sky tract centers then one
186 is arbitrarily chosen.
188 - The default implementation is not very efficient; subclasses may wish
192 If tracts do not cover the whole sky then the returned tract may not
195 distTractInfoList = []
196 for i, tractInfo
in enumerate(self):
197 angSep = coord.separation(tractInfo.getCtrCoord()).asDegrees()
199 distTractInfoList.append((angSep, i, tractInfo))
200 distTractInfoList.sort()
201 return distTractInfoList[0][2]
204 """Find array of tract IDs with vectorized operations (where supported).
206 If a given sky map does not support vectorized operations, then a loop
207 over findTract will be called.
212 Array of Right Ascension. Units are radians unless
215 Array of Declination. Units are radians unless
217 degrees : `bool`, optional
218 Input ra, dec arrays are degrees if True.
222 tractId : `np.ndarray`
227 - If coord is equidistant between multiple sky tract centers then one
228 is arbitrarily chosen.
231 If tracts do not cover the whole sky then the returned tract may not
232 include the given ra/dec.
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),
238 tractId = np.array([self.
findTractfindTract(coord).getId()
for coord
in coords])
243 """Find tracts and patches that overlap a region.
247 coordList : `list` of `lsst.geom.SpherePoint`
248 List of ICRS sky coordinates to search for.
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.
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).
264 for tractInfo
in self:
265 patchList = tractInfo.findPatchList(coordList)
267 retList.append((tractInfo, patchList))
271 """Find closest tract and patches that overlap coordinates.
275 coordList : `lsst.geom.SpherePoint`
276 List of ICRS sky coordinates to search for.
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.
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))
303 return hash(self.
getSha1getSha1())
307 return self.
getSha1getSha1() == other.getSha1()
308 except AttributeError:
309 return NotImplemented
312 return not (self == other)
315 """Write information about a sky map to supplied log
320 Log object that information about skymap will be written
322 log.info(
"sky map has %s tracts" % (len(self),))
323 for tractInfo
in self:
324 wcs = tractInfo.getWcs()
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]))
339 """Return a SHA1 hash that uniquely identifies this SkyMap instance.
344 A 20-byte hash that uniquely identifies this SkyMap instance.
348 Subclasses should almost always override ``updateSha1`` instead of
349 this function to add subclass-specific state to the hash.
351 if self.
_sha1_sha1
is None:
352 sha1 = hashlib.sha1()
353 sha1.update(
type(self).__name__.encode(
'utf-8'))
355 sha1.update(configPacked)
357 self.
_sha1_sha1 = sha1.digest()
358 return self.
_sha1_sha1
361 """Add subclass-specific state or configuration options to the SHA1.
365 sha1 : `hashlib.sha1`
366 A hashlib object on which `update` can be called to add
367 additional state to the hash.
371 This method is conceptually "protected" : it should be reimplemented by
372 all subclasses, but called only by the base class implementation of
375 raise NotImplementedError()
377 SKYMAP_RUN_COLLECTION_NAME =
"skymaps"
379 SKYMAP_DATASET_TYPE_NAME =
"skyMap"
382 """Add skymap, tract, and patch Dimension entries to the given Gen3
388 The name of the skymap.
389 butler : `lsst.daf.butler.Butler`
390 The butler to add to.
394 lsst.daf.butler.registry.ConflictingDefinitionError
395 Raised if a different skymap exists with the same name.
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.
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.
415 for tractInfo
in self:
416 nx, ny = tractInfo.getNumPatches()
417 nxMax =
max(nxMax, nx)
418 nyMax =
max(nyMax, ny)
419 region = tractInfo.getOuterSkyPolygon()
421 tractRecords.append({
423 "tract": tractInfo.getId(),
425 "ra": centroid.getRa().asDegrees(),
426 "dec": centroid.getDec().asDegrees(),
428 for patchInfo
in tractInfo:
429 cellX, cellY = patchInfo.getIndex()
430 patchRecords.append({
432 "tract": tractInfo.getId(),
433 "patch": tractInfo.getSequentialPatchIndex(patchInfo),
436 "region": patchInfo.getOuterSkyPolygon(tractInfo.getWcs()),
441 "tract_max": len(self),
442 "patch_nx_max": nxMax,
443 "patch_ny_max": nyMax,
449 from lsst.daf.butler
import DatasetType
450 from lsst.daf.butler.registry
import ConflictingDefinitionError
451 datasetType = DatasetType(
453 dimensions=[
"skymap"],
454 storageClass=
"SkyMap",
455 universe=butler.registry.dimensions
457 butler.registry.registerDatasetType(datasetType)
458 with butler.transaction():
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."
466 butler.registry.insertDimensionData(
"tract", *tractRecords)
467 butler.registry.insertDimensionData(
"patch", *patchRecords)
A class representing an angle.
A floating-point coordinate rectangle geometry.
Point in an unspecified spherical coordinate system.
def getPatchInnerDimensions(self)
def setPatchBorder(self, value)
def setPatchInnerDimensions(self, value)
def findTract(self, coord)
def register(self, name, butler)
def updateSha1(self, sha1)
def findClosestTractPatchList(self, coordList)
string SKYMAP_RUN_COLLECTION_NAME
string SKYMAP_DATASET_TYPE_NAME
def logSkyMapInfo(self, log)
def findTractPatchList(self, coordList)
def findTractIdArray(self, ra, dec, degrees=False)
def __init__(self, config=None)
def __getitem__(self, ind)