LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
baseSkyMap.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 """
23 @todo
24 - Consider tweaking pixel scale so the average scale is as specified, rather than the scale at the center
25 """
26 import lsst.pex.config as pexConfig
27 import lsst.afw.geom as afwGeom
28 from . import detail
29 
30 __all__ = ["BaseSkyMap"]
31 
32 class BaseSkyMapConfig(pexConfig.Config):
33  patchInnerDimensions = pexConfig.ListField(
34  doc = "dimensions of inner region of patches (x,y pixels)",
35  dtype = int,
36  length = 2,
37  default = (4000, 4000),
38  )
39  patchBorder = pexConfig.Field(
40  doc = "border between patch inner and outer bbox (pixels)",
41  dtype = int,
42  default = 100,
43  )
44  tractOverlap = pexConfig.Field(
45  doc = "minimum overlap between adjacent sky tracts, on the sky (deg)",
46  dtype = float,
47  default = 1.0,
48  )
49  pixelScale = pexConfig.Field(
50  doc = "nominal pixel scale (arcsec/pixel)",
51  dtype = float,
52  default = 0.333
53  )
54  projection = pexConfig.Field(
55  doc = """one of the FITS WCS projection codes, such as:
56  - STG: stereographic projection
57  - MOL: Molleweide's projection
58  - TAN: tangent-plane projection
59  """,
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 
70 class BaseSkyMap(object):
71  """A collection of overlapping Tracts that map part or all of the sky.
72 
73  See TractInfo for more information.
74 
75  BaseSkyMap is an abstract base class. Subclasses must do the following:
76  @li define __init__ and have it construct the TractInfo objects and put them in _tractInfoList
77  @li define __getstate__ and __setstate__ to allow pickling (the butler saves sky maps using pickle);
78  see DodecaSkyMap for an example of how to do this. (Most of that code could be moved
79  into this base class, but that would make it harder to handle older versions of pickle data.)
80  """
81  ConfigClass = BaseSkyMapConfig
82  def __init__(self, config=None):
83  """Construct a BaseSkyMap
84 
85  @param[in] config: an instance of self.ConfigClass; if None the default config is used
86  """
87  if config is None:
88  config = self.ConfigClass()
89  config.freeze() # just to be sure, e.g. for pickling
90  self.config = config
91  self._tractInfoList = []
93  pixelScale = afwGeom.Angle(self.config.pixelScale, afwGeom.arcseconds),
94  projection = self.config.projection,
95  rotation = afwGeom.Angle(self.config.rotation, afwGeom.degrees),
96  )
97 
98  def findTract(self, coord):
99  """Find the tract whose center is nearest the specified coord.
100 
101  @param[in] coord: sky coordinate (afwCoord.Coord)
102  @return TractInfo of tract whose center is nearest the specified coord
103 
104  @warning:
105  - if tracts do not cover the whole sky then the returned tract may not include the coord
106 
107  @note
108  - This routine will be more efficient if coord is ICRS.
109  - If coord is equidistant between multiple sky tract centers then one is arbitrarily chosen.
110  - The default implementation is not very efficient; subclasses may wish to override.
111  """
112  icrsCoord = coord.toIcrs()
113  distTractInfoList = []
114  for tractInfo in self:
115  angSep = icrsCoord.angularSeparation(tractInfo.getCtrCoord()).asDegrees()
116  distTractInfoList.append((angSep, tractInfo))
117  distTractInfoList.sort()
118  return distTractInfoList[0][1]
119 
120  def findTractPatchList(self, coordList):
121  """Find tracts and patches that overlap a region
122 
123  @param[in] coordList: list of sky coordinates (afwCoord.Coord)
124  @return list of (TractInfo, list of PatchInfo) for tracts and patches that contain,
125  or may contain, the specified region. The list will be empty if there is no overlap.
126 
127  @warning this uses a naive algorithm that may find some tracts and patches that do not overlap
128  the region (especially if the region is not a rectangle aligned along patch x,y).
129  """
130  retList = []
131  for tractInfo in self:
132  patchList = tractInfo.findPatchList(coordList)
133  if patchList:
134  retList.append((tractInfo, patchList))
135  return retList
136 
137  def __getitem__(self, ind):
138  return self._tractInfoList[ind]
139 
140  def __iter__(self):
141  return iter(self._tractInfoList)
142 
143  def __len__(self):
144  return len(self._tractInfoList)
int iter