Coverage for python/lsst/pipe/tasks/makeDiscreteSkyMap.py: 94%

1# This file is part of pipe_tasks. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22__all__ = ["MakeDiscreteSkyMapConfig", "MakeDiscreteSkyMapTask"] 

23 

24import lsst.sphgeom 

25 

26import lsst.geom as geom 

27import lsst.pex.config as pexConfig 

28import lsst.pipe.base as pipeBase 

29from lsst.skymap import DiscreteSkyMap, BaseSkyMap 

30from lsst.utils.timer import timeMethod 

31 

32 

33class MakeDiscreteSkyMapConfig(pexConfig.Config): 

34 """Config for MakeDiscreteSkyMapTask. 

35 """ 

36 

37 coaddName = pexConfig.Field( 

38 doc="coadd name, e.g. deep, goodSeeing, chiSquared", 

39 dtype=str, 

40 default="deep", 

41 ) 

42 skyMap = pexConfig.ConfigField( 

43 dtype=BaseSkyMap.ConfigClass, 

44 doc="SkyMap configuration parameters, excluding position and radius" 

45 ) 

46 borderSize = pexConfig.Field( 

47 doc="additional border added to the bounding box of the calexps, in degrees", 

48 dtype=float, 

49 default=0.0 

50 ) 

51 doAppend = pexConfig.Field( 

52 doc="append another tract to an existing DiscreteSkyMap on disk, if present?", 

53 dtype=bool, 

54 default=False 

55 ) 

56 doWrite = pexConfig.Field( 

57 doc="persist the skyMap?", 

58 dtype=bool, 

59 default=True, 

60 ) 

61 

62 def setDefaults(self): 

63 self.skyMap.tractOverlap = 0.0 

64 

65 

66class MakeDiscreteSkyMapTask(pipeBase.Task): 

67 """Make a DiscreteSkyMap in a repository, using the bounding box of a set of calexps. 

68 

69 The command-line and run signatures and config are sufficiently different from MakeSkyMapTask 

70 that we don't inherit from it, but it is a replacement, so we use the same config/metadata names. 

71 """ 

72 

73 ConfigClass = MakeDiscreteSkyMapConfig 

74 _DefaultName = "makeDiscreteSkyMap" 

75 

76 @timeMethod 

77 def run(self, wcs_bbox_tuple_list, oldSkyMap=None): 

78 """Make a SkyMap from the bounds of the given set of calexp metadata. 

79 

80 Parameters 

81 ---------- 

82 wcs_bbox_tuple_list : `iterable` 

83 A list of tuples with each element expected to be a (Wcs, Box2I) pair. 

84 oldSkyMap : `lsst.skymap.DiscreteSkyMap`, optional 

85 The SkyMap to extend if appending. 

86 

87 Returns 

88 ------- 

89 skyMap : `lsst.pipe.base.Struct` 

90 Sky map returned as a struct with attributes: 

91 

92 ``skyMap`` 

93 The returned SkyMap (`lsst.skyMap.SkyMap`). 

94 """ 

95 self.log.info("Extracting bounding boxes of %d images", len(wcs_bbox_tuple_list)) 

96 points = [] 

97 for wcs, boxI in wcs_bbox_tuple_list: 

98 boxD = geom.Box2D(boxI) 

99 points.extend(wcs.pixelToSky(corner).getVector() for corner in boxD.getCorners()) 

100 if len(points) == 0: 100 ↛ 101line 100 didn't jump to line 101 because the condition on line 100 was never true

101 raise RuntimeError("No data found from which to compute convex hull") 

102 self.log.info("Computing spherical convex hull") 

103 polygon = lsst.sphgeom.ConvexPolygon.convexHull(points) 

104 if polygon is None: 104 ↛ 105line 104 didn't jump to line 105 because the condition on line 104 was never true

105 raise RuntimeError( 

106 "Failed to compute convex hull of the vertices of all calexp bounding boxes; " 

107 "they may not be hemispherical." 

108 ) 

109 circle = polygon.getBoundingCircle() 

110 

111 skyMapConfig = DiscreteSkyMap.ConfigClass() 

112 if oldSkyMap: 

113 skyMapConfig.raList.extend(oldSkyMap.config.raList) 

114 skyMapConfig.decList.extend(oldSkyMap.config.decList) 

115 skyMapConfig.radiusList.extend(oldSkyMap.config.radiusList) 

116 configIntersection = {k: getattr(self.config.skyMap, k) 

117 for k in self.config.skyMap.toDict() 

118 if k in skyMapConfig} 

119 skyMapConfig.update(**configIntersection) 

120 circleCenter = lsst.sphgeom.LonLat(circle.getCenter()) 

121 skyMapConfig.raList.append(circleCenter[0].asDegrees()) 

122 skyMapConfig.decList.append(circleCenter[1].asDegrees()) 

123 circleRadiusDeg = circle.getOpeningAngle().asDegrees() 

124 skyMapConfig.radiusList.append(circleRadiusDeg + self.config.borderSize) 

125 skyMap = DiscreteSkyMap(skyMapConfig) 

126 

127 for tractInfo in skyMap: 

128 wcs = tractInfo.getWcs() 

129 posBox = geom.Box2D(tractInfo.getBBox()) 

130 pixelPosList = ( 

131 posBox.getMin(), 

132 geom.Point2D(posBox.getMaxX(), posBox.getMinY()), 

133 posBox.getMax(), 

134 geom.Point2D(posBox.getMinX(), posBox.getMaxY()), 

135 ) 

136 skyPosList = [wcs.pixelToSky(pos).getPosition(geom.degrees) for pos in pixelPosList] 

137 posStrList = ["(%0.3f, %0.3f)" % tuple(skyPos) for skyPos in skyPosList] 

138 self.log.info("tract %s has corners %s (RA, Dec deg) and %s x %s patches", 

139 tractInfo.getId(), ", ".join(posStrList), 

140 tractInfo.getNumPatches()[0], tractInfo.getNumPatches()[1]) 

141 return pipeBase.Struct( 

142 skyMap=skyMap 

143 )