LSSTApplications  20.0.0
LSSTDataManagementBasePackage
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.ConfigClass()
106  config.freeze() # just to be sure, e.g. for pickling
107  self.config = config
108  self._tractInfoList = []
110  pixelScale=Angle(self.config.pixelScale, arcseconds),
111  projection=self.config.projection,
112  rotation=Angle(self.config.rotation, degrees),
113  )
114  self._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.findTract(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[ind]
202 
203  def __iter__(self):
204  return iter(self._tractInfoList)
205 
206  def __len__(self):
207  return len(self._tractInfoList)
208 
209  def __hash__(self):
210  return hash(self.getSha1())
211 
212  def __eq__(self, other):
213  try:
214  return self.getSha1() == 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 is None:
259  sha1 = hashlib.sha1()
260  sha1.update(type(self).__name__.encode('utf-8'))
261  configPacked = struct.pack(
262  "<iiidd3sd",
263  self.config.patchInnerDimensions[0],
264  self.config.patchInnerDimensions[1],
265  self.config.patchBorder,
266  self.config.tractOverlap,
267  self.config.pixelScale,
268  self.config.projection.encode('ascii'),
269  self.config.rotation
270  )
271  sha1.update(configPacked)
272  self.updateSha1(sha1)
273  self._sha1 = sha1.digest()
274  return self._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  def register(self, name, registry):
294  """Add SkyMap, Tract, and Patch Dimension entries to the given Gen3
295  Butler Registry.
296 
297  Parameters
298  ----------
299  name : `str`
300  The name of the skymap.
301  registry : `lsst.daf.butler.Registry`
302  The registry to add to.
303  """
304  nxMax = 0
305  nyMax = 0
306  records = {
307  "skymap": [],
308  "tract": [],
309  "patch": [],
310  }
311  for tractInfo in self:
312  nx, ny = tractInfo.getNumPatches()
313  nxMax = max(nxMax, nx)
314  nyMax = max(nyMax, ny)
315  region = tractInfo.getOuterSkyPolygon()
316  centroid = SpherePoint(region.getCentroid())
317  records["tract"].append({
318  "skymap": name,
319  "tract": tractInfo.getId(),
320  "region": region,
321  "ra": centroid.getRa().asDegrees(),
322  "dec": centroid.getDec().asDegrees(),
323  })
324  for patchInfo in tractInfo:
325  cellX, cellY = patchInfo.getIndex()
326  records["patch"].append({
327  "skymap": name,
328  "tract": tractInfo.getId(),
329  "patch": tractInfo.getSequentialPatchIndex(patchInfo),
330  "cell_x": cellX,
331  "cell_y": cellY,
332  "region": patchInfo.getOuterSkyPolygon(tractInfo.getWcs()),
333  })
334  records["skymap"].append({
335  "skymap": name,
336  "hash": self.getSha1(),
337  "tract_max": len(self),
338  "patch_nx_max": nxMax,
339  "patch_ny_max": nyMax,
340  })
341  with registry.transaction():
342  for dimension, recordsForDimension in records.items():
343  registry.insertDimensionData(dimension, *recordsForDimension)
lsst.skymap.baseSkyMap.BaseSkyMap.logSkyMapInfo
def logSkyMapInfo(self, log)
Definition: baseSkyMap.py:221
lsst.skymap.baseSkyMap.BaseSkyMap.__len__
def __len__(self)
Definition: baseSkyMap.py:206
lsst.skymap.baseSkyMap.BaseSkyMap.__ne__
def __ne__(self, other)
Definition: baseSkyMap.py:218
lsst.skymap.baseSkyMap.BaseSkyMap.config
config
Definition: baseSkyMap.py:107
lsst.skymap.baseSkyMap.BaseSkyMap._sha1
_sha1
Definition: baseSkyMap.py:114
lsst.skymap.baseSkyMap.BaseSkyMap.getSha1
def getSha1(self)
Definition: baseSkyMap.py:245
ast::append
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
lsst.skymap.baseSkyMap.BaseSkyMap._wcsFactory
_wcsFactory
Definition: baseSkyMap.py:109
lsst.skymap.baseSkyMap.BaseSkyMap.register
def register(self, name, registry)
Definition: baseSkyMap.py:293
lsst.skymap.baseSkyMap.BaseSkyMap._tractInfoList
_tractInfoList
Definition: baseSkyMap.py:108
lsst.skymap.baseSkyMap.BaseSkyMap.findClosestTractPatchList
def findClosestTractPatchList(self, coordList)
Definition: baseSkyMap.py:177
lsst.skymap.baseSkyMap.BaseSkyMapConfig
Definition: baseSkyMap.py:38
lsst.skymap.baseSkyMap.BaseSkyMap.__eq__
def __eq__(self, other)
Definition: baseSkyMap.py:212
lsst.skymap.baseSkyMap.BaseSkyMap.__init__
def __init__(self, config=None)
Definition: baseSkyMap.py:103
max
int max
Definition: BoundedField.cc:104
lsst.skymap.baseSkyMap.BaseSkyMap
Definition: baseSkyMap.py:75
lsst.skymap.baseSkyMap.BaseSkyMap.ConfigClass
ConfigClass
Definition: baseSkyMap.py:101
lsst::geom
Definition: geomOperators.dox:4
type
table::Key< int > type
Definition: Detector.cc:163
lsst.skymap.baseSkyMap.BaseSkyMap.__iter__
def __iter__(self)
Definition: baseSkyMap.py:203
lsst::geom::Point< double, 2 >
lsst.skymap.baseSkyMap.BaseSkyMap.findTractPatchList
def findTractPatchList(self, coordList)
Definition: baseSkyMap.py:149
lsst.skymap.baseSkyMap.BaseSkyMap.findTract
def findTract(self, coord)
Definition: baseSkyMap.py:116
lsst::geom::Angle
A class representing an angle.
Definition: Angle.h:127
lsst::geom::SpherePoint
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
lsst::geom::Box2D
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
lsst.skymap.detail.wcsFactory.WcsFactory
Definition: wcsFactory.py:29
lsst.skymap.baseSkyMap.BaseSkyMap.__hash__
def __hash__(self)
Definition: baseSkyMap.py:209
astshim.fitsChanContinued.iter
def iter(self)
Definition: fitsChanContinued.py:88
lsst.skymap.baseSkyMap.BaseSkyMap.__getitem__
def __getitem__(self, ind)
Definition: baseSkyMap.py:200
lsst.skymap.baseSkyMap.BaseSkyMap.updateSha1
def updateSha1(self, sha1)
Definition: baseSkyMap.py:276